.. _models:

******************************************************************************
Models
******************************************************************************

Due to computatonnal issues associated with the free surface Navier-Stokes equations, 
geophysical flows simulations are often carried out with reduced complexity models.
With vertically averaged models such as the Saint-Venant system, the water height being
an unknown of the system, there is no need to deal with moving meshes.
Freshkiss3D's numerical schemes are based on multilayer Saint-Venant systems
with mass-exchanges (i.e. layer-averaged Euler or layer-averaged Navier-Stokes systems)
of reduced complexity with respect to the Navier-Stokes equations but enriched with 
respect to the classic Shallow Water system.
Layer-averaged models implemented in Freshkiss3D include topography, viscosity and
friction source terms as well as variable density and tracer convection-diffusion.

.. figure:: _static/images/models.png
    :width: 850px
    :align: center

In this section the Navier-Stokes with free surface system is presented in the 
introduction and then the focus is put on layer-averaged models implemented
in Freshkiss3D. Given the complexity of the numerical methods involved, this
documentation only scratches the surface of Freshkiss3D's solver and
it is highly recommended to refer to the :ref:`references` section.


I - System of equations for free surface flows
================================================================================

The fluid domain along :math:`z` axis is delimited by free suface :math:`\eta(t,x,y)`
and topography :math:`z_b(x,y)` and water depth is given by :math:`h(t,x,y)=\eta(t,x,y)-z_b(x,y)`.

.. figure:: _static/images/free_surf.png
    :width: 400px
    :align: center

    Fluid domain in (x,z) plane

I.1 - Incompressible Navier-Stokes equations with free surface
--------------------------------------------------------------------------------

We consider the classical incompressible Navier-Stokes system:

.. math::
    \begin{eqnarray}
    \nabla. {\bf U}&=&0 \label{eq:nsi1}\\
    \rho \left( \partial_t \bf U + \bf U. \nabla \bf U \right) &=&
    \nabla. \boldsymbol{\Sigma} +\rho {\bf g}
    \end{eqnarray}

where the stress tensor :math:`\boldsymbol{\Sigma}` is given by:

.. math::
    \boldsymbol{\Sigma} = -p \boldsymbol{I} + \mu [\nabla {\bf U}+(\nabla {\bf U})^T]

with :math:`{\bf U}(t,x,y,z)=(u,v,w)^T` the velocity,
:math:`p` the fluid pressure,
:math:`{\bf g}= ( 0,0,-g)^T` the gravity forces, :math:`\mu` the dynamic viscosity
and :math:`\rho` the density.
The symbol :math:`\nabla` denotes :math:`\nabla = \left(\partial_x, \partial_y, \partial_z \right)^{\top}`.
Under the hydrostatic assumption, which consists in neglecting the vertical acceleration of the fluid,
and under the Boussinesq assumption, which consists in neglecting the density variations in 
the acceleration terms, the system of equations becomes:

.. important:: **(NS) The Navier-Stokes system with free surface**:

    .. math::
        \begin{align}
        \nabla. {\bf U} &= 0 \label{eq:incomp}\\
        \partial_t {\bf u} + \nabla_{x,y}. ( {\bf u}\otimes {\bf u})+
        \partial_z( {\bf u}w ) &= - \dfrac{1}{\rho_0} \nabla_{x,y} p + \nu \Delta {\bf u} \label{eq:mvthm}\\
        \partial_z p &= - \rho g \label{eq:verticalihm} 
        \end{align}

with :math:`{\bf u}(t,x,y,z)=(u,v)^T` the horizontal velocity,
:math:`\nu` the kinematic viscosity defined by :math:`\mu/\rho_0` and
:math:`\nabla_{x,y}` the horizontal projection of :math:`\nabla`.


I.2 - Tracer 
--------------------------------------------------------------------------------

The governing convection-diffusion equation for a tracer :math:`T` inside the fluid domain is given by:

.. important:: **(T) Tracer convection-diffusion equation**:

    .. math::
        \partial_t T + \nabla_{x,y}. ( {\bf u} T)+
        \partial_z (wT) = \nu_T \Delta T

with :math:`\nu_T` is defined by :math:`D_T/\rho_0` with :math:`D_T` the diffusion coefficient.

.. note::
    In Freshkiss3D
    
    * All default parameters are in the *conf_dict.py* file in the *problem* directory and could be changed (see details in :ref:`api`, ``problem``)
    * Tracers are optional and can be added to the simulation with the ``fk.Tracer`` class.
    * By default the density :math:`\rho` (see **Navier-Stokes system**) is set to :math:`\rho_0` but can be set as a function of the tracer with the ``variable_density`` key of ``numerical_parameters`` in the ``fk.Problem`` class. 
    * The default density state laws included in Freshkiss3D can either take into account the effect of temperature, salinity and particle suspensions. See the parameters of the ``fk.Tracer`` class in the :ref:`api`.

I.3 - Boundary conditions and initial conditions
--------------------------------------------------------------------------------

To complete the **Navier-Stokes system with free surface** boundary conditions
and initial conditions have to be provided.
Bottom and free surface kinematic and dynamic boundary conditions are defined in Freshkiss3D
by:

.. math::
    \begin{cases}
    \partial_t \eta + {\bf U}_s . {\bf n}_s = 0 \\
    {\bf U}_b . {\bf n}_b = 0
    \end{cases}

and:

.. math::
    \begin{cases}
    \boldsymbol{\Sigma}. \mathbf{n}_s = -p^{atm}.\mathbf{n}_s + \tau_w . \mathbf{t}_s \\
    (\boldsymbol{\Sigma}. \mathbf{n}_b).  \mathbf{t}_b = - \kappa_b(h,\mathbf{u}) \mathbf{u} .\mathbf{t}_b
    \end{cases}

with 
:math:`{\bf U}_s=U(t,x,y,z=\eta)`,
:math:`{\bf U}_b=U(t,x,y,z=Z)`,
:math:`{\bf n}_s` and :math:`{\bf n}_b` the unit outwards normals at the free surface and bottom,
:math:`p^{atm}` the atmospheric pressure,
:math:`\tau_w` the wind shear stress at the free surface and
:math:`\kappa_b(h,\mathbf{u})` the friction coefficient.

.. note::
    In Freshkiss3D (cf. :ref:`api`)

    * :math:`\tau_w` is optional and can be defined with the ``fk.WindEffect`` external effect class.
    * :math:`\kappa_b(h,\mathbf{u})` is optional and can be defined with the ``fk.Friction`` class.

On horizontal boundaries several types of boundary condition can be prescribed depending
on the flow regime. The way of imposing these conditions depends on the layer-averaged model's
formulation and will not be detailed here. Additionally one has also to prescribe initial conditions
to complete the model.

II - Layer-averaged models
================================================================================

Layer-averaged models are the core of Freshkiss3D and allow a 3D modelling of free
surface flows. Depending on the initial hypotheses 
(i.e. user selected options) the models solved by Freshkiss3D behind
the curtain can vary. In this section we briefly present some mathematical aspects of 
the vertical discretization as well as the different models obtained depending on various 
hypotheses. Some hints on how to select specific options are given with links to the :ref:`api`.

II.1 - Layerwise discretization
--------------------------------------------------------------------------------

We consider a discretization of the fluid domain by layers where the
layer :math:`\alpha` contains the points of coordinates :math:`(x,y,z)` with :math:`z \in
L_\alpha(t,x,y)= (z_{\alpha-1/2},z_{\alpha+1/2})` and
:math:`\{z_{\alpha+1/2}\}_{{\alpha} \in \{1,\dots,N \}}` is defined for :math:`\alpha \in \{1,\dots,N \}` by:

.. math::
    \left\{\begin{array}{l}
    z_{\alpha+1/2}(t,x,y) = z_b(x,y) + \sum_{j=1}^\alpha h_j(t,x,y) \\
    h_\alpha(t,x,y) = z_{\alpha+1/2}(t,x,y) - z_{\alpha-1/2}(t,x,y)=l_\alpha h(t,x,y) \quad \text{with} \quad \sum_{\alpha=1}^N l_\alpha =1.
    \end{array}\right.

.. figure:: _static/images/free_discrete.png
    :width: 550px
    :align: center

    Fluid domain and layerwise discretization

Using the previous notations, let us consider the space
:math:`\Bbb{P}_{0}` of piecewise constant functions defined by

.. math::
    \Bbb{P}_{0} = \left\{ \mathbb{1}_{z \in L_\alpha(t,x,y)}(z),\quad
      \alpha\in\{1,\ldots,N\}\right\},

where :math:`\mathbb{1}_{z \in L_\alpha(t,x,y)}(z)` is the characteristic
function of the layer :math:`L_\alpha(t,x,y)`. Using this formalism, the
projection of :math:`u`, :math:`v`, :math:`w` and :math:`T` on :math:`\Bbb{P}_{0}`
is a piecewise constant function defined by:

.. math::
    X(t,x,y,\{z_\alpha\})  = \sum_{\alpha=1}^N
    \mathbb{1}_{[z_{\alpha-1/2},z_{\alpha+1/2}]}(z)X_\alpha(t,x,y),

Note that the free surface changing over time, the :math:`\Bbb{P}_{0}` space as well
as the piecewise constant function projections are time dependent.

II.2 - The layer-averaged Euler system
--------------------------------------------------------------------------------

Let us first neglect the viscosity and assume that the density is constant (i.e. :math:`\rho =\rho_0`).
Using the space :math:`\Bbb{P}_0` and the above 
decomposition, the Galerkin approximation of the incompressible and hydrostatic
Euler equations leads to the system:

.. important:: **(LE.1) The layer-averaged Euler system**:

    for :math:`\alpha \in \{1,\dots,N \}`:

    .. math::
        \begin{eqnarray}
        &&\sum_{\alpha=1}^N \frac{\partial  h_\alpha  }{\partial t } + \sum_{\alpha=1}^N  \nabla_{x,y} .   (h_\alpha  {\bf u} _\alpha)  =0\\
        &&\frac{\partial  h_{\alpha}{\bf u}_{\alpha}}{\partial t} +  \nabla_{x,y} .\left(h_{\alpha} {\bf u}_{\alpha} \otimes {\bf u}_{\alpha} \right) + \nabla_{x,y} \bigl(\frac{g}{2} h h_{\alpha}\bigr) = -gh_\alpha \nabla_{x,y} z_b \nonumber\\
        && \qquad + {\bf u}_{\alpha+1/2} G_{\alpha+1/2} - {\bf u}_{\alpha-1/2} G_{\alpha-1/2}
        \end{eqnarray}

    The quantity :math:`G_{\alpha+1/2}` (resp. :math:`G_{\alpha-1/2}`) corresponds to the mass
    exchange accross the interface :math:`z_{\alpha+1/2}` (resp. :math:`z_{\alpha-1/2}`)
    and :math:`G_{\alpha+1/2}` is defined by:

    .. math::
        \begin{eqnarray}
        %G_{1/2} & = & 0\nonumber\\
        G_{\alpha+1/2}
        %& = &\frac{\partial  z_{\alpha+1/2}}{\partial t} + u_{\alpha+1/2} \frac{\partial  z_{\alpha+1/2}}{\partial x} + v_{\alpha+1/2} \frac{\partial  z_{\alpha+1/2}}{\partial y} - w(x,y,z_{\alpha+1/2},t)\nonumber\\
        & = & \sum_{j=1}^\alpha \left(\frac{\partial h_j}{\partial t }
        +  \nabla_{x,y} . (h_j {\bf u}_j) \right)
         =  -\sum_{j=1}^N \Bigl(\sum_{p=1}^\alpha l_p - \mathbb{1}_{j\leq
              \alpha}\Bigr)\nabla_{x,y} . (h_j {\bf u}_j),
        \end{eqnarray}

    for :math:`\alpha=1,\ldots,N`. The velocities at the interfaces :math:`{\bf u}_{\alpha+1/2}` are defined by:

    .. math::
        {\bf u}_{\alpha+1/2} =
        \left\{\begin{array}{ll}
        {\bf u}_\alpha & if \;G_{\alpha+1/2} \leq 0\\
        {\bf u}_{\alpha+1} & if \;G_{\alpha+1/2} > 0
        \end{array}\right.

For smooth solutions, the system **(LE.1)** satisfies an energy balance (see :ref:`references` for complete proof).
Note that the vertical velocity is not an unknown of the system but can be computed considering the 
divergence free condition and we obtain the following piecewise approximation of the vertical velocity:

.. math::
    w_{\alpha} = k_{\alpha} - z_{\alpha} \nabla_{x,y}.{\bf u}_{\alpha} \quad \text{with:} \quad \begin{cases}
    k_1 = \nabla_{x,y} . (z_b{\bf u}_{1})\\
    k_{\alpha+1} = k_{\alpha} + \nabla_{x,y} . (z_{\alpha+1/2} ({\bf u}_{\alpha+1} - {\bf u}_{\alpha} ))
    \end{cases}

Also note that when the number of layers is equal to one, the classical Saint-Venant system with topography is obtained:

.. important:: **(SV) The Saint-Venant system**:

    .. math::
        \begin{eqnarray}
        &&\frac{\partial h}{\partial t} + \nabla_{x,y} . (h {\bf u})= 0 \\
        &&\frac{\partial h {\bf u}}{\partial t} + \nabla_{x,y} . \left(h {\bf u} \otimes {\bf u} + g \frac{h^2}{2} \right) = -gh \nabla_{x,y} z_b
        \end{eqnarray}

.. note::
    In Freshkiss3D (cf. :ref:`api`)

    * **The layer-averaged Euler system is solved by default when no other option is selected**
    * Layers are defined with the ``fk.Layer`` class which is initialized with a user defined number of layers:``NL``.


Applying the same decomposition to the tracer equation gives:

.. important:: **(LT) The layer-averaged tracer equation**:

    for :math:`\alpha \in \{1,\dots,N \}`:

    .. math::
        \begin{eqnarray}
        \frac{\partial  h_\alpha T_\alpha }{\partial t } + \nabla_{x,y} . (h_\alpha  {\bf u}_\alpha T_\alpha)  = T_{\alpha+1/2} G_{\alpha+1/2} - T_{\alpha-1/2} G_{\alpha-1/2}
        \end{eqnarray}

    where :math:`G_{\alpha+1/2}` and :math:`G_{\alpha-1/2}` are defined as in **(II.2)**.


II.3 - The layer-averaged Euler system with variable density
--------------------------------------------------------------------------------

As mentioned before, the variable density model relies on the Boussinesq assumption, 
which means that the variable density only appears in the gravity term.
When the density is not constant, integrating the vertical pressure gradient leads to a slighly
different model as one cannot assume :math:`p=\rho_0 g (\eta-z)` anymore. Now
:math:`p = g \int_{z}^{\eta} \rho dz`. Again, using piecewise approximation one can 
write the pressure as:

.. math::
    p(x,y,z,t) = g \left( \rho_{\alpha} (z_{\alpha+1/2}-z) + \sum_{j=\alpha+1}^{N} \rho_j h_j \right)

The Galerkin approximation of the incompressible and hydrostatic
Euler equations with variable density leads to the system:

.. important:: **(LE.2) The layer-averaged Euler system with variable density**:

    for :math:`\alpha \in \{1,\dots,N \}`:

    .. math::
        \begin{eqnarray}
        &&\sum_{\alpha=1}^N \frac{\partial  h_\alpha  }{\partial t } + \sum_{\alpha=1}^N  \nabla_{x,y} .   (h_\alpha  {\bf u} _\alpha)=0\\
        &&\frac{\partial  h_{\alpha}{\bf u}_{\alpha}}{\partial t} +  \nabla_{x,y} .\left(h_{\alpha} {\bf u}_{\alpha} \otimes {\bf u}_{\alpha} \right) + \nabla_{x,y} (h_{\alpha} p_{\alpha}) = {\bf u}_{\alpha+1/2} G_{\alpha+1/2} - {\bf u}_{\alpha-1/2} G_{\alpha-1/2}\\
        && \qquad + p_{\alpha+1/2} \nabla_{x,y} z_{\alpha+1/2} - p_{\alpha-1/2} \nabla_{x,y} z_{\alpha-1/2}\\
        &&\frac{\partial  h_\alpha T_\alpha }{\partial t } + \nabla_{x,y} . (h_\alpha  {\bf u}_\alpha T_\alpha)  = T_{\alpha+1/2} G_{\alpha+1/2} - T_{\alpha-1/2} G_{\alpha-1/2}
        \end{eqnarray}

    where :math:`G_{\alpha+1/2}` (resp. :math:`G_{\alpha-1/2}`),
    :math:`z_{\alpha+1/2}` (resp. :math:`z_{\alpha-1/2}`),
    :math:`{\bf u}_{\alpha+1/2}` (resp. :math:`{\bf u}_{\alpha-1/2}`) are defined as in **(LE.1)**,
    :math:`\rho_{\alpha}=f(T_{\alpha})` with :math:`f` a state law and
    :math:`p_{\alpha}`, :math:`p_{\alpha+1/2}` are defined by:

    .. math::
        \begin{cases}
        p_{\alpha} = \frac{1}{h_{\alpha}} \int_{z_{\alpha-1/2}}^{z_{\alpha+1/2}} p(x,y,z,t)dz = g \left( \frac{\rho_{\alpha} h_{\alpha}}{2} + \sum_{j=\alpha +1}^{N} \rho_j h_j \right)\\
        p_{\alpha+1/2} = p(x,y,z_{\alpha+1/2},t) = g \sum_{j=\alpha +1}^{N} \rho_j h_j
        \end{cases}


.. note:: 
    In Freshkiss3D (cf. :ref:`api`)

    * As mentionned before, variable density model can be activated by setting the ``variable_density`` key of ``numerical_parameters`` to ``True``.
    * The ``ipres`` key of ``numerical_parameters`` allows the user to select the pressure term formulation of **(LE.2)** instead of **(LE.1)** even if :math:`\rho` is set to constant, in which case :math:`\rho_{\alpha}` is replaced in all equations by :math:`\rho_0` and the tracer is optional. Note that when ``boussinesq`` is prescribed, ``ipres`` is automatically set to ``True`` and tracer is mandatory.


II.4 - The layer-averaged Navier-Stokes system
--------------------------------------------------------------------------------

Let us now apply the same process to the complete Navier-Stokes with free surface system **(NS)**, 
so that viscosity and friction terms are added to the system.
Note that it can be done with or without variable density but only
the constant :math:`\rho` case is presented for sake of simplicity.
From now on we denote :math:`\bf \Sigma` as the x,y projection of the 3D stress tensor.

.. important:: **(LNS.1) The layer-averaged Navier-Stokes system**:

    for :math:`\alpha \in \{1,\dots,N \}`:

    .. math::
        \begin{eqnarray}
        &&\sum_{\alpha=1}^N \frac{\partial  h_\alpha  }{\partial t } + \sum_{\alpha=1}^N  \nabla_{x,y} .   (h_\alpha  {\bf u} _\alpha)  =0\\
        &&\frac{\partial  h_{\alpha}{\bf u}_{\alpha}}{\partial t} +  \nabla_{x,y} .\left(h_{\alpha} {\bf u}_{\alpha} \otimes {\bf u}_{\alpha} \right) + \nabla_{x,y} \bigl(\frac{g}{2} h h_{\alpha}\bigr) = -gh_\alpha \nabla_{x,y} z_b \nonumber\\
        && \qquad + {\bf u}_{\alpha+1/2} G_{\alpha+1/2} - {\bf u}_{\alpha-1/2} G_{\alpha-1/2} + \nabla_{x,y} .   \bigl( h_\alpha {\bf \Sigma}_\alpha \bigr) \nonumber\\
        && \qquad - {\bf \Sigma}_{\alpha+1/2} \nabla_{x,y} z_{\alpha+1/2} + {\bf \Sigma}_{\alpha-1/2} \nabla_{x,y} z_{\alpha-1/2} \nonumber\\
        && \qquad +2\mu_{\alpha+1/2} \frac{{\bf u}_{{\alpha}+1}-{\bf u}_{\alpha}}{h_{{\alpha}+1}+h_{\alpha}} - 2\mu_{\alpha-1/2} \frac{{\bf u}_{\alpha}-{\bf
        u}_{{\alpha}-1}}{h_{\alpha}+h_{{\alpha}-1}} -\kappa_{\alpha}{\bf u}_{\alpha}
        \end{eqnarray}

    where :math:`G_{\alpha+1/2}` (resp. :math:`G_{\alpha-1/2}`),
    :math:`z_{\alpha+1/2}` (resp. :math:`z_{\alpha-1/2}`),
    :math:`{\bf u}_{\alpha+1/2}` (resp. :math:`{\bf u}_{\alpha-1/2}`) are defined as in **(LE.1)** and the stress tensor terms are defined by:

    .. math::
        \begin{eqnarray}
        {\bf \Sigma}_{\alpha+1/2} & = & \begin{pmatrix}
        \Sigma_{xx,\alpha+1/2} & \Sigma_{xy,\alpha+1/2}\\
        \Sigma_{yx,\alpha+1/2} & \Sigma_{yy,\alpha+1/2}
        \end{pmatrix}\\
        \Sigma_{xx,\alpha+1/2} & = & \frac{\mu_{\alpha+1/2}}{h_{\alpha+1}+h_\alpha}\bigl( h_{\alpha}\frac{\partial
          u_{\alpha}}{\partial x} + h_{\alpha+1}\frac{\partial u_{\alpha+1}}{\partial x}  \bigr) - 2 \mu_{\alpha+1/2}\frac{\partial z_{\alpha+1/2}}{\partial x} \frac{u_{\alpha+1} - u_\alpha}{h_{\alpha+1} + h_\alpha}\\
        \Sigma_{xy,\alpha+1/2} & = & \frac{\mu_{\alpha+1/2}}{h_{\alpha+1}+h_\alpha}\bigl( h_{\alpha}\frac{\partial
          u_{\alpha}}{\partial y} + h_{\alpha+1}\frac{\partial u_{\alpha+1}}{\partial y}  \bigr) - 2 \mu_{\alpha+1/2}\frac{\partial z_{\alpha+1/2}}{\partial y} \frac{u_{\alpha+1} - u_\alpha}{h_{\alpha+1} + h_\alpha}\\
        {\bf \Sigma}_{\alpha} & = & \begin{pmatrix}
        \Sigma_{xx,\alpha} & \Sigma_{xy,\alpha}\\
        \Sigma_{yx,\alpha} & \Sigma_{yy,\alpha}
        \end{pmatrix} = \frac{{\bf \Sigma}_{\alpha+1/2}+ {\bf \Sigma}_{\alpha-1/2}}{2}
        \end{eqnarray}

    with the friction and viscosity coefficients given by:

    .. math::
        \begin{equation}
        \kappa_{\alpha}=
        \left\{
        \begin{array}{l}
        \kappa \quad {\rm if} \quad \alpha =1\\
        0 \quad {\rm if} \quad\alpha \ne 1
        \end{array}
        \right.
        \quad
        \mu_{\alpha+1/2}=
        \left\{
        \begin{array}{l}
        0 \quad {\rm if} \quad \alpha =0\\
        \mu \quad {\rm if} \quad\alpha \in \{1,...,N-1 \} \\
        0\quad {\rm if} \quad\alpha=N
        \end{array}\right.
        \end{equation}

    with :math:`\kappa` the Navier friction coefficient and :math:`\mu` the dynamic viscosity

For smooth solutions, the system **(LNS.1)** satisfies an energy balance (see :ref:`references` for complete proof). The viscous terms in the layer-averaged Navier-Stokes sytem are difficult to discretise so
only a simplified version of this model is implemented into Freshkiss3D (see :ref:`references` section for details).

Also note that when the number of layers is equal to one, the viscous Saint-Venant system with topography is obtained. 
In this case, the friction coefficient can be chosen amongst a wide variety of empirical laws. Note that only the constant :math:`\kappa` (Navier friction law) is relevant when :math:`NL>1`.

.. important:: **(VSV) The viscous Saint-Venant system**:

    .. math::
        \begin{eqnarray}
        &&\frac{\partial h}{\partial t} + \nabla_{x,y} . (h {\bf u})= 0 \\
        &&\frac{\partial h {\bf u}}{\partial t} + \nabla_{x,y} . \left(h {\bf u} \otimes {\bf u} + g \frac{h^2}{2} \right) = -gh \nabla_{x,y} z_b - \kappa(h,{\bf u}).{\bf u} + \nabla_{x,y} . (h {\bf \Sigma})
        \end{eqnarray}

    with :math:`\kappa(h,{\bf u})`:

    * **Navier**: :math:`\kappa(h,{\bf u}) = \kappa`
    * **Chézy**: :math:`\kappa(h,{\bf u}) = \frac{g |{\bf u}|}{C_h^2}`
    * **Strikler**: :math:`\kappa(h,{\bf u}) = \frac{g}{S_t^2} \left( \frac{|{\bf u}|}{h^{1/3}} \right)`
    * **Manning**: :math:`\kappa(h,{\bf u}) = g M_a^2  \left( \frac{|{\bf u}|}{h^{1/3}} \right)`


.. note::
    In Freshkiss3D (cf. :ref:`api`)

    * Viscosity can be toggled by setting the ``horizontal_viscosity`` and ``vertical_viscosity`` keys of the ``physical_parameters`` dictionary.
    * Friction can be toggled by setting the ``friction_coeff`` and ``friction_law`` keys of the ``physical_parameters`` dictionary. Note that friction can also be added to the layer-averaged Euler system when viscosity is set to zero.
    * ``friction_law`` is set to ``'Navier'`` by default
    * When ``NL=1`` user can choose between several ``friction_law`` including 

.. important:: **The Saint-Venant system with coriolis force**:

    .. math::
        \frac{\partial  h  }{\partial t } + \nabla_{x,y}\cdot   (h  {\bf u} )  =0, \\
        \frac{\partial (h{\bf u})}{\partial t} +  \nabla_{x,y} \cdot\left(
        h {\bf u} \otimes {\bf u} \right) +
        \nabla_{x,y} \bigl(\frac{g}{2} h^2\bigr) = -gh \nabla_{x,y} z_b
        -\Omega h{\bf u}^\perp.


.. note::
    In Freshkiss3D (cf. :ref:`api`)
 
    * Coriolis force might be applied by setting the ``coriolis_effect`` keys of the ``physical_parameters`` dictionary. The user can choose to pass :math:`\Omega(x, y)` to ``Coriolis_coeff`` attribut of ``fk.primitive`` by an array or a constant or a function.




III - Numerical scheme
================================================================================

The layer-averaged Navier-Stokes system has the form:

.. math::
    \partial_t U +  \nabla_{x,y} . F(U) = S_b(U) + S_e(U) + S_{v,f}(U)

where the state vector is defined by:

.. math::
    U=\left(h, q_{x,1}, \ldots, q_{x,N}, q_{y,1}, \ldots, q_{y,N}\right)^T

with :math:`q_{x,\alpha}=h_\alpha u_{\alpha}`, :math:`q_{y,\alpha}=h_\alpha v_{\alpha}`.
:math:`F(U)` stands for the fluxes of the conservative part,
and :math:`S_b(U)`, :math:`S_e(U)` and :math:`S_{v,f}(U)` are the source terms, 
representing respectively the topography effects,
the momentum exchanges and the viscous and friction effects.

III.1 - Time discretization
--------------------------------------------------------------------------------

Let us denote discrete times by :math:`t^n` and :math:`t^{n+1}=t^n+\Delta t^n`. For
the time discretisation of the layer-averaged Navier-Stokes system 
a time splitting technique is adopted to treat the conservative part, the vertical exchanges
and the viscous part successively:

.. math::
    \begin{cases}
    U^* = U^n - \Delta t^n \left( \nabla_{x,y} . F(U^n) - S_b(U^n) \right) \\
    \tilde{U} = U^* + \Delta t^n S_e(\tilde{U}, G(U^{n}))\\
    U^{n+1} = \tilde{U} + \Delta t^n S_{v,f}(U^{n+1}, \tilde{U})
    \end{cases}

The conservative part is fully explicit and includes the topography source term which 
is treated with a well-balanced hydrostatic reconstruction scheme. 
The vertical exchanges terms are semi-implicit, :math:`G` being computed with :math:`U^{n}`.
The viscous and friction effects are semi-implicit as well, the horizontal viscosity being fully
explicit and the vertical viscosity and friction (treated at the same time) fully implicit. 
The coriolis force is an semi-implicit.

Noting :math:`U^{n+1} = U^n + \Delta t^n f(U^n)` the fully explicit first order scheme, 
the second order scheme extension implemented in Freshkiss3D, based on a modified version
of the Heun scheme, is:

.. math::
    \begin{cases}
    U^{(1)} =  U^n + \Delta t_1^n f(U^n) \\
    U^{(2)} =  U^{(1)} + \Delta t_2^{n} f(U^{(1)}) \\
    U^{n+1} =  (1 - \gamma) U^n + \gamma U^{(2)}
    \end{cases}

with:

.. math::
    \Delta t^n = \frac{2\Delta t_1^n\Delta t _2^n}{\Delta t_1^n + \Delta t_2^n}, \quad \gamma = \frac{(\Delta t^n)^2}{2 \Delta t_1^n \Delta t_2^n}

Since :math:`\gamma \geq 0`, :math:`U^{n+1}` is a convex combination of :math:`U^n`
and :math:`U^{(2)}` so the scheme preserves the positivity.
For the previous relations :math:`\Delta t_1^{n}` and :math:`\Delta t_2^{n}` respectively satisfy the
CFL conditions associated with :math:`U^n` and :math:`U^{(1)}`.

.. note::
    In Freshkiss3D (cf. :ref:`api`)

    * Time discretization parameters are passed to ``fk.Problem`` class via the ``fk.Simutime`` class.
    * By default time step :math:`\Delta t^n` is computed with a CFL condition to ensure stability. The default CFL number is set to 0.9 but can be changed with the parameter ``cfl_adjuster`` key of the ``numerical_parameters`` dictionary. Although it is not recommended, the user can force a constant time step by passing ``constant_dt`` to ``fk.Simutime``.
    * In ``fk.Simutime`` you can set the ``second_order`` to ``True`` to select second order time scheme.
    * By default semi-implicit methods are used for :math:`S_e` and :math:`S_{v,f}` but the user can switch back to fully explicit methods by setting ``implicit_exchanges`` and ``implicit_vertical_viscosity`` parameters to ``False`` in the ``fk.Problem`` class.

III.2 - Space discretization
--------------------------------------------------------------------------------

Let :math:`\Omega` denote the computational domain with boundary
:math:`\Gamma`, which we assume is polygonal. Let :math:`T_h` be a triangulation
of :math:`\Omega` for which the vertices are denoted by :math:`P_i` with :math:`S_i` the
set of interior nodes and :math:`G_i` the set of boundary nodes.
For the space discretization of the system, we use an hybrid method. A finite volume
technique is used for the conservative part and exchange terms and a :math:`\Bbb{P}_1`  finite
element approach on :math:`T_h` for the viscous part.

We recall now the general formalism of finite volumes on unstructured
meshes. The dual cells :math:`C_i` are obtained by joining the centers of mass of the
triangles surrounding each vertex :math:`P_i`. We use the following
notations:

* :math:`K_i`, set of subscripts of nodes :math:`P_j` surrounding :math:`P_i`,
* :math:`|C_i|`, area of :math:`C_i`,
* :math:`\Gamma_{ij}`, boundary edge between the cells :math:`C_i` and :math:`C_j`,
* :math:`L_{ij}`, length of :math:`\Gamma_{ij}`,
* :math:`{\bf n}_{ij}`, unit normal to :math:`\Gamma_{ij}`, outward to :math:`C_i` (:math:`{\bf
  n}_{ji}= -{\bf n}_{ij}`).

If :math:`P_i` is a node belonging to the boundary :math:`\Gamma`, we join the centers of mass of the triangles adjacent to the
boundary to the middle of the edge belonging to :math:`\Gamma` and we denote

* :math:`\Gamma_i`, the two edges of :math:`C_i` belonging to :math:`\Gamma`,
* :math:`L_i`, length of :math:`\Gamma_i` (for sake of simplicity we assume in
  the following that :math:`L_i = 0` if :math:`P_i` does not belong to :math:`\Gamma`),
* :math:`{\bf n}_i`, the unit outward normal defined by averaging the two adjacent normals.

.. figure:: _static/images/mesh_2d.png
    :width: 500px
    :align: center

    Dual cell :math:`C_i` and bondary cell.

.. note::
    In Freshkiss3D (cf. :ref:`api`)

    * The triangular mesh is passed to the ``fk.Problem`` class via the ``fk.TriangularMesh`` class.
    * The ``fk.TriangularMesh`` class can be defined from a file with ``fk.TriangularMesh.from_msh_file('file.msh')``
    * Freshkiss3D mesh IO supports 3 different mesh formats for your mesh file: FreeFem++ ``.msh``, Medit ``.mesh`` and Gmsh ``.msh`` (format msh2).
    * See the **Meshes** section of the :ref:`gallery_of_examples` to learn how to manipulate meshes with Freshkiss3D


We define the piecewise constant functions :math:`U^n(x,y)` on cells :math:`C_i` 
corresponding to time :math:`t^n` and :math:`z_b(x,y)` as:

.. math::
	U^n(x,y)= U^n_i,\quad z_b(x,y)= z_{b,i},\quad\mbox{ for } (x,y)\in C_i

with :math:`U^n_i = (h^n_i,q^n_{x,1,i},\ldots,q^n_{x,N,i},q^n_{y,1,i},\ldots,q^n_{y,N,i})^T` i.e. :math:`U^n_i \approx \frac{1}{|C_i|}\int_{C_i} U(t^n,x,y)dxdy`, :math:`z_{b,i} \approx \frac{1}{|C_i|}\int_{C_i} z_b(x,y)dxdy`. We will also use the notation :math:`U^n_{\alpha,i} \approx \frac{1}{|C_i|}\int_{C_i}U_\alpha(t^n,x,y)dxdy`.
A finite volume scheme for solving the conservative part is: 

.. math::
	U^{*}_i=U^n_i- \sum_{j\in K_i} \sigma_{i,j}{\cal F}_{i,j} - \sigma_i {\cal F}_{e,i}


.. math::
    \sum_{j\in K_i} L_{i,j} {\cal F}_{i,j} \approx \int_{C_i}
    \nabla_{x,y} . F(U^n) dx dy,
    
with :math:`\sigma_{i,j} = \frac{\Delta t^nL_{i,j}}{|C_i|},\quad \sigma_{i} = \frac{\Delta t^nL_{i}}{|C_i|}`. Here we consider first-order explicit schemes where:

.. math::
    {\cal F}_{i,j} = F(U^n_i,U^n_j,z^n_{b,i}-z^n_{b,j},{\bf n}_{i,j}) \quad \text{and} \quad {\cal F}_{e,i} = F(U^n_i,U^n_{e,i},{\bf n}_{i})

The term :math:`\cal{F}_{i,j}` denotes an interpolation of the normal component of the flux 
:math:`F(U^n). {\bf n}_{i,j}` along the edge :math:`C_{i,j}`. The functions
:math:`F(U^n_i,U^n_j,z^n_{b,i} - z^n_{b,j},{\bf n}_{i,j})\in {\cal R}^{2N+1}` are the numerical 
fluxes which are computed using a kinetic interpretation of the system (see :ref:`references`).
The source terms discretization (vertical exchanges, viscosity,...) as well as the second order scheme
have been omitted but are detailed in reference papers (see :ref:`references`).

.. note::
    In Freshkiss3D (cf. :ref:`api`)

    * In ``fk.Problem`` you can either set ``second_order`` parameter or ``second_order`` key of ``numerical_parameters`` dictionary to ``True`` in order to select the second order space scheme.

III.3 - Initial conditions
--------------------------------------------------------------------------------

Initial conditions are needed to define initial state vector:

.. math::
    U^0=\left(h^0, q^0_{x,1}, \ldots, q^0_{x,N}, q^0_{y,1}, \ldots, q^0_{y,N}\right)^T

as well as initial values for topography: :math:`z^0_b`. By default all variables are set to zero in Freshkiss3D if no initial values are provided.

.. note::
    In Freshkiss3D (cf. :ref:`api`):

    * The topography is initialized when calling the ``fk.Layer`` class presented above.
    * The state vector containing ``U``, ``V``, ``H``, ``HL``, ``QX``, and ``QY`` variables is defined with the ``fk.Primitive`` class. Initialization takes place when declaring this class for the first time using ``init`` parameters, the types of which can be constant float, vector or function.
    * As mentioned before the tracer ``T`` is defined and initialized in ``fk.Tracer`` class.


III.4 - Boundary conditions
--------------------------------------------------------------------------------

In order to complete the numerical scheme, one has to provide boundary conditions, i.e. to define the fluxes :math:`{\cal F}_{e,i} = F(U^n_i,U^n_{e,i},{\bf n}_{i})`. 

.. note::
    In Freshkiss3D (cf. :ref:`api`), boundary conditions are selected with the following class:

    * ``fk.FluvialFlowRate``: fluvial regime fluid boundary condition with given flow rate.
    * ``fk.FluvialHeight``: fluvial regime fluid boundary condition with given water depth.
    * ``fk.TorrentialOut``: torrential regime fluid boundary condition with free outlet.
    * ``fk.Slide``: solid wall boundary with slip condition :math:`{\bf u} \cdot {\bf n} = 0`.
    * ``fk.DampingOut``: damp waves with coefficient :math:`\theta` ; :math:`\partial_t \bf X = -\theta ({\bf X} - {\bf X}^0)`.
    * Time varying boundary conditions can be passed to ``fk.FluvialFlowRate`` and ``fk.FluvialHeight`` using ``fk.TimeDependentFlowRate`` and ``fk.TimeDependentHeight``.