.. Hydroflow documentation file, created by
   Pierre-Yves Gousenbourger on Wed Aug 2, 2023.

.. _solvers:

Available solvers
=================

All models implemented in Hydroflow are solved using finite volumes. With this numerical method, the domain is divided into computational cells in which the variables are constant (for a first-order scheme). If the mesh is unstructured and made of triangles, all cells will have two (at the boundary) or three neighbouring cells (inside the domain). These adjacent cells share common interfaces or intercell edges. At every iteration, the governing variables are updated for each cell with the following:

.. math::
   
   \bf{U}^{n+1}_{i} = \bf{U}^{n}_{i}-\frac{\Delta t}{\Omega _{i}}\sum_{j=1}^{nb}\bf{T}^{-1}_{j}\bf{F}^{*}_{j}(\bar{\bf{U}}_{j})\bf{L}_{j}+\bf{S}_{i}\Delta t


where :math:`\bf{U}` is the vector of the governing variables, :math:`n` is the temporal exponent, :math:`i` is the spatial index, :math:`\Delta t` is the timestep, :math:`\Omega _{i}` is the area of cell :math:`i`, :math:`nb` is the number of cells of the domain, :math:`\bf{T}` is a rotation matrix, :math:`F^{*}_{j}` is the flux at the intercell edge :math:`j`, :math:`L_{j}` is the length of the :math:`j`-th edge, and :math:`\bf{S}` are the source terms.

With this numerical method, discontinuities at the interfaces emerge that form local Riemann problems. Solving these Riemann problems results in fluxes accrosse the cell interfaces. The computation of the flux :math:`\bf{F}` through the interface will be done according to a chosen solver. Different types of solvers can be implemented for each model. Nevertheless, their implementation can vary from one model to another, depending on the physics of the considered models. Every solver implementation is hence specific to a model. Yet, for brevity reasons, we will here keep to the basics and use the Saint-Venant model as an example.

HLLC
----

The HLL solver was developed by Harten, Lax & van Leer (1983). It was then extended to the HLLC variation by Toro *et al.* (1994). This latter solver uses the characteristics of the model to compute the fluxes. More particularly, it uses wave speed estimates, based on the eigenvalues of the Jacobian matrix :math:`\bf{J} = \frac{\partial \bf{F(U)}}{\partial \bf{U}}` of the model, as weights to balance the fluxes from the left and right cells. Without getting into the details, the Saint-Venant model in two dimensions, considered as an augmented 1D problem at the local scale of an interface, has three characteristics:

.. math::
   
   \lambda_1 = u_n-c


.. math::
   
   \lambda_2 = c


.. math::
   
   \lambda_3 = u_n+c


with :math:`u_n` being the velocity of water in the direction normal to the edge and :math:`c=\sqrt{gh}` the celerity or the speed of propagation of the information in water where :math:`g` is the gravitational acceleration and :math:`h` the water depth.

These characteristics are ranked in increasing order of speed. Considering a flow going in the positive direction of a cartesian reference frame centred around the centre of the edge (:math:`u_n` is positive), they can all be positive if the flow is supercritical (i.e. :math:`u_n > c`). If the flow is subcritical (i.e. :math:`u_n < c`), :math:`\lambda_2` and :math:`\lambda_3` will be positive while :math:`\lambda_1` will be negative.

The wave speed estimates are then defined as:

.. math::
   
   \lambda^{-} = \max(\lambda_{1,L},\lambda_{1,R},0)


.. math::
   
   \lambda^{+} = \max(\lambda_{3,L},\lambda_{3,R},0)


Finally, the fluxes at the interfaces are computed in this way:

.. math::
   
   q^{*}_n = \frac{\lambda^{+}q_{n,L}-\lambda^{-}q_{n,R}+\lambda^{+}\lambda^{-}(z_{w,R}-z_{w,L})}{\lambda^{+}-\lambda^{-}}


.. math::
   
   \sigma^{*}_n = \frac{\lambda^{+}\sigma_{n,L}-\lambda^{-}\sigma_{n,R}+\lambda^{+}\lambda^{-}(q_{n,R}-q_{n,L})}{\lambda^{+}-\lambda^{-}}


.. math::
   
   \mu^{*}_{nt} =
   \begin{cases}
      v_{t,L}q^{*}_{n} & \text{if $u_{n} \geq 0$}\\
      v_{t,R}q^{*}_{n} & \text{if $u_{n} < 0$}
   \end{cases}


with :math:`n` and :math:`t` being respectively the indices referring to the directions normal and tangential to the edge, :math:`\sigma_n` and :math:`\mu_{nt}` the normal and tangential momentum fluxes to the edge respectively, and :math:`v_{t}` the velocity in the direction tangential to the edge.

The different formulation of the solver for the computation of the tangential momentum flux is due to the consideration of the contact discontinuity in the HLL solver, extending it as a result to the HLLC approach. 

Variations
^^^^^^^^^^

- The SVE model adds the Exner equation to those of the SV model. Extending the HLLC solver to the fully coupled SVE model leads to an overly diffusive behaviour of the model. Treating the fluxes associated with the shallow water equations differently from that of the Exner equation, in a Weakly Coupled way, showed better results. More information about these two variations can be found in Meurice & Soares-Frazão (2020).


Roe
---

In construction


Recapitulative table
--------------------

====== ================ ==========================
Solver Concerned models Variation
====== ================ ==========================
HLLC   SV/2P2L          / 
HLLC   SVE              Coupled or Weakly Coupled
Roe 
====== ================ ==========================


References
----------

Harten, A., Lax, P. D., & Leer, B. V. (1983). On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. *SIAM review*, 25(1), 35-61.

Meurice, R., & Soares-Frazão, S. (2020). A 2D HLL-based weakly coupled model for transient flows on mobile beds. *Journal of Hydroinformatics*, 22(5), 1351-1369.

Toro, E. F., Spruce, M., & Speares, W. (1994). Restoration of the contact surface in the HLL-Riemann solver. *Shock waves*, 4(1), 25-34.