Uhaina is based on the high order finite elements library Aerosol, which provides arbitrary high order continuous and discontinuous elements on hybrid meshes (triangles and quadrangle in two dimensions).
The discretization of the hyperbolic part of the equations presented in section Governing equations is done using the discontinuous Galerkin method.
Discontinuous Galerkin discretization
The discontinuous Galerkin methods was introduced in [1] for neutrons transport equation. Since its popularity grew, and it is now extended to many other applications as in [2] and [3]. To compute the flux function, we use a classical numerical flux such as Lax Friedrich flux or HLL flux.
Lax Friedrich flux
The Lax Friedrich flux is a well known numerical flux, extremely simple to implement.
It consist in a centered part and a diffusive one :
H_{LF} = \dfrac{F^{-} + F^{+}}{2} - \dfrac{ | \lambda | }{2} ( \check{W}_{h}^{+} - \check{W}_{h}^{-} )
Figure 1 : Riemann problem at boundary.
Boundary conditions
For now, the boundary condition are imposed using the same numerical flux describe earlier, where the boundary state is given by the user. The main difficulty here is to prescribe the boundary state outside of the domain (U_{b}
). In fact, it often happens that the available data is not complete. Also we want our boundary condition to prevent reflection of outgoing wave into the domain.
Classical subcritical boundary
Considering the case where not all the variables of the system are known, we can reconstruct the boundary state using the Riemann invariant R^{+} = u + 2c = cte
.
For example, if we know the water eight h
, we can compute the velocity:
u_{bound} = u_{inner} - 2 (\sqrt{gh_{bound}} - \sqrt{gh_{inner}})
If this kind of boundary condition allow to impose correct water depth, it is highly reflective.
LODI boundary condition
A characteristic based boundary condition was implemented referred as LODI for \textit{local associated one-dimensional inviscid}, looking at the work in [4]. Here are detail of the method.
Considering the one dimensional shallow water equations, we derived the following characteristic system :
\begin{aligned}
& \partial_{t} (u - 2c) + (u - c) \partial_{x} (u - 2c) = 0 \\
& \partial_{t} (u + 2c) + (u + c) \partial_{x} (u + 2c) = 0
\end{aligned}
Defining the amplitudes of characteristic wave associated with each characteristic velocity :
\begin{aligned}
& L_{1} = \lambda_{1} \partial_{x} (u - 2c) \\
& L_{2} = \lambda_{2} \partial_{x} (u + 2c)
\end{aligned}
By simple combination of the characteristic system, we derived the following system :
\begin{aligned}
& \partial_{t}(u) + \dfrac{1}{2} (L_{1} + L_{2}) = 0 \\
& \partial_{t}(c) + \dfrac{1}{4} (L_{2} - L_{1}) = 0
\end{aligned}
To prevent wave to enter the domain we simply kill it's amplitude by setting L_{1} = 0
. The boundary state is then obtain by solving the previous system. Unfortunately, if this method gives non reflective boundary condition it does not allow correct imposition of water depth.
Riemann invariant extrapolation
This method consist in computing the boundary state using Riemann invariant :
\begin{aligned}
& h_{bound} = \frac{1}{g} \left(\frac{R^{+} - R^{-}}{4}\right)^{2} \\
& u_{bound} = \frac{1}{2}(R^{+} + R^{-})
\end{aligned}
The outgoing Riemann invariant R^{+}
is simply computed using interior state at boundary. For the incoming one, R^{-}
, we use the imposed water depth with velocity reconstructed using the reference state at boundary :
u_{R^{-}} = u_{0} - 2 (\sqrt{gh_{imposed}} - \sqrt{gh_{0}})
This method allow correct imposition of the water depth and non-reflexion of outgoing waves.
Non-reflective characteristic boundary
To be continue
Numerical issues
The numerical resolution of the nonlinear shallow water or Green-Naghdi equations raises a number of problems.
Dry cells
A well known problem concern the ability to handle dry areas during the computation. Indeed, the well balanced Shallow Water equations is valid under the condition h > 0
. Numerically, this issue is tackle by forcing this quantity to be greater than an tolerance value \epsilon
, as done in [5] or [6] for example. When h
reach \epsilon
the others quantities are set to zero. Unfortunately, this is generally not sufficient to ensure stability of the method.
In Uhaina we prevent this behavior by applying a positivity limiter describe in [6]. It consist in relaxing the solution around cell average value if the cell exhibit negative h
value at Zhang and Shu quadrature points. This formula is construct by transformation of the tensor product of a Gauss-Legendre and Gauss-Lobatto quadrature.
For each element, we compute the minimum value of h
in each points in S_{Zhang}
the set of quadrature points:
h_{min} = \min_{x \in S_{Zhang}} h(x)
If h_{min} \leq \epsilon
we perform a relaxation of the solution state around \bar{h} = \int_{\Omega_{e}} h d \Omega
and \bar{q} = \int_{\Omega_{e}} q d \Omega
:
\begin{aligned}
& \tilde{h}(x) = \Theta (h(x) - \bar{h}) + \bar{h}, \quad \forall x \in S_{Gauss} \\
& \tilde{q}(x) = \Theta (q(x) - \bar{q}) + \bar{q}, \quad \forall x \in S_{Gauss}
\end{aligned}
With
\Theta = \min \left(1, \dfrac{\epsilon - \bar{h}}{h_{min} - \bar{h}} \right)
Doing so the scheme remain conservative and the positivity of the solution is ensure at Gauss quadrature points for faces integrals.
Remark : the positivity of h
is not ensure at Gauss quadrature points used for volume integration, which can be problematic for cell flux computation depending on how we compute the debit on almost dry areas. To handle this issue we first add the Gauss quadrature points when computing h_{min}
. Recently we modified the way the debit is compute, by first calculate the velocity vector depending on the h
value, so that :
u = 0 \quad if \quad h < \epsilon
Remark : additionally, when we perform the relaxation on the debit, we define an other parameter \epsilon_{2} > \epsilon
so that :
\tilde{q} = \dfrac{\tilde{h} - \epsilon}{\epsilon_{2} - \epsilon} \tilde{q}
We do so, in order to prevent excessive velocity value near the shoreline (where h
is close to \epsilon
) that could be harmful when evaluating wave speed : \lambda = max(u \pm c)
.
Remark : we first tried to solve shoreline instability with a method that consist in locally decreasing the order of the approximation by killing high order mode on modal basis.
This method, which require to perform two change of basis for each concerned element, create oscillations around limited areas. I explained this by the fact that we strongly modify the solution, inducing a "jump" at the interface between two elements.
This solution is not used anymore in the code.
Shock limitation
An other stability issue arise when dealing with shock, with high order approximation of the solution. In fact the scheme becomes unstable in discontinuities. This behavior is well know as the Gibbs phenomenon. In order to make the code more robust and to be able to compute breaking wave, we add a limiter to prevent oscillations. The method used is the one proposed by [7], which consist in adding an artificial viscosity to the well balanced Shallow water system:
\dfrac{\partial W}{\partial t} + \nabla \cdot \widetilde{F}(W, z_{b}) = \widetilde{S}(W, z_{b}) + \nu \Delta W
based on the amount of entropy of the solution.
Let define the energy density of a column of fluid :
E = \dfrac{1}{2} \left( h u^{2} + g(\eta^{2} + z^{2}) \right)
and the associate conservation equation :
D = \dfrac{\partial E}{\partial t} + \nabla \cdot f(E) \leq 0
(E, f)
is a satisfying mathematical entropy pair for the shallow water system and the could be used for limitation.
From a physical perspective, the residual D
is always zero except in discontinuities where it is negative, which make it an effective criterion to detect shocks.
Based on this quantity and on [7], we define the artificial viscosity :
\nu (D) = \alpha \nu_{max} min \left(1.0, \delta x \dfrac{|D|}{\beta N} \right)
Where \nu_{max}
is of the same order as the dissipation of a first order numerical scheme :
\nu_{max} = \delta x max |\lambda|
\beta
and \alpha
are tunable parameters. They can be modified in Uhaina through the configuration file (see Configuration files). N
a scaling parameter who aims to normalize |D|
. We proposed to take N = g^{\frac{3}{2}} h^{\frac{5}{2}}
, the theoretical dissipation in breaking wave.
We then add a viscous term \nu \Delta W
to the right hand side of the well-balanced ShallowWater equation, which is discretized using the discontinuous Galerkin method presented above :
\sum_{\Omega_{e} \in M_{h}} \int_{\Omega_{e}}{\nu \nabla v_{h} \cdot \nabla W_{h}}
Remark : one could remark that we remove the face integral that came with the integration by part and do not consider any penalization operator. In fact, as we are interested only in adding dissipation to the system, we are not concern by the numerical consistency of this term.
Remark : the evaluation of D
is done once at each time step independently of the order of the Runge-Kutta scheme. It is performed using the same discontinuous Galerkin discretization as for the shallow water system and a second order BDF scheme.
References
[1] T.R. Reed, Wm H and Hill. “Triangular Mesh Method for the Neutron Transport Equation”. Los Alamos Report LA-UR-73- 479 (1973) (cit. on p. 6).
[2] F. Bassi and S. Rebay. “A High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier–Stokes Equations”. Journal of Computational Physics 131.2 (1997), pp. 267–279 (cit. on p. 6).
[3] B. Cockburn and C. W. Shu. “Runge – Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems”. Journal of Scientific Computing 16.3 (2001), p. 89 (cit. on pp. 6, 19).
[4] T. J. Poinsot and Sanjiva K. Lele. “Boundary conditions for direct simulations of compressible viscous flows”. Journal of Computational Physics 101 (1992), pp. 104–129 (cit. on p. 11).
[5] A. Duran “Numerical simulation of depth-averaged flow models : a class of Finite Volume and discontinuous Galerkin approaches.” PhD thesis. 2014 (cit. on pp. 13, 32).
[6] X. Zhang, Y. Xia, and C-W Shu. “Maximum-Principle-Satisfying and Positivity-Preserving High Order Discontinuous Galerkin Schemes for Conservation Laws on Triangular Meshes”. Journal of Scientific Computing (2012), pp. 29–62 (cit. on p. 13).
[7] J. L. Guermond, R. Pasquetti, and B. Popov. “Entropy viscosity method for nonlinear conservation laws”. Journal of Computational Physics 230.11 (2011), pp. 4248–4267 (cit. on pp. 14, 15).