/*!
\page numericalPage The numerical methods 


CoreFlows proposes a variety of finite volume methods (see \ref leveque for an introduction). The method  can be explicit or implicit (enum \ref TimeScheme), and the convection fluxes can be approximated using the  \ref roe,  \ref vfroe or \ref vffc formulations (enum \ref NonLinearFormulation).

The basic method for fluid models is the \ref roe scheme with entropic correction (see \ref kieu) and source upwinding (see\ref wbsourceterms). In order to increase precision it is possible to use centered, staggered, pressureCorrection or lowMach schemes through the enum \ref SpaceScheme.

The finite volume discretisation allows an easy handling of general geometries and meshes generated by \ref salome .

Explicit schemes are used in general for fast dynamics solved with small time steps while implicit schemes allow the use of large time step to quickly reach the stationary regime. The implicit schemes result in nonlinear systems that are solved using a Newton type method.

The upwind scheme is the basic scheme but options are available to use a centered scheme (second order in space) or entropic corrections.

Our models can be written in generic form as a nonlinear system of balance laws:

\f{eqnarray*} 
 
\frac{\partial U}{\partial t} + \nabla \cdot \left(\mathcal{F}^{conv}(U)\right)+\nabla \cdot \left(\mathcal{F}^{diff}(U)\right) = 0, \label{equationnNS}
\f}


where 
-\f$U\f$ is the vector of conservative unknowns, 
-\f$\mathcal{F}^{conv}\f$ is the convective flux 
- and \f$\mathcal{F}^{diff}\f$ the diffusive flux.

We decompose the computational domain into \f$N\f$ disjoint cells \f$C_i\f$ with volume \f$v_i\f$.
+ Two neighboring cells \f$C_i\f$ and \f$C_j\f$ have a common boundary \f$\partial C_{ij}\f$ with area \f$s_{ij}\f$.
+ We denote \f$N(i)\f$ the set of neighbors of a given cell \f$C_i\f$ and \f$\vec n_{ij}\f$ the exterior unit normal vector of \f$\partial C_{ij}\f$ . 

Integrating the system (\ref NStokesEq) over \f$C_{i}\f$ and setting \f$U_i(t)=  \frac{1}{v_i} \int_{C_i} U(x,t) dx\f$, the semi-discrete equations can be written:

\f{eqnarray*}
 \frac{\mathrm{d} U_i}{\mathrm{d} t} + \sum_{j \in N(i)} \frac{s_{ij}}{v_i}\left(\overrightarrow \Phi^{conv}_{ij} +  \overrightarrow{\Phi}^{diff}_{ij}\right)= S_i(U,x).\label{eq:numer scheme}
\f}
with: 
- \f$\overrightarrow{\Phi}^{con}_{ij}= \frac{1}{s_{ij}}\int_{\partial C_{ij}}\mathcal F^{conv}(U).\vec n_{ij}ds \f$,
- \f$\overrightarrow \Phi^{diff}_{ij}= \frac{1}{s_{ij}}\int_{\partial C_{ij}}\mathcal {F}^{diff}(U).\vec n_{ij}ds \f$

To approximate the convection numerical flux \f$\overrightarrow \Phi^{conv}_{ij}\f$ we solve an  approximate Riemann problem at the interface \f$\partial C_{ij}\f$. There are three possible formulations for the convection fluxes. 
- Using the \ref roe local linearisation of the fluxes, we obtain the following formula:
\f{eqnarray*}
\overrightarrow{\Phi}^{conv, Roe}_{ij}&= &\frac{\mathcal{F}^{conv}(U_i) + \mathcal{F}^{conv}(U_j)}{2} \vec{n}_{ij}- \mathcal{D}(U_i,U_j) \frac{U_j-U_i}{2}\\
&=&\mathcal{F}^{conv}(U_i) \vec{n}_{ij} + A^-(U_i,U_j) (U_j - U_i),\label{eq:flux roe}
\f}
- Using the \ref vfroe local linearisation of the fluxes, we obtain the following formula:
\f{eqnarray*}
\overrightarrow{\Phi}^{conv, VFRoe}_{ij}&= &\mathcal{F}^{conv}\left(\frac{U_i + U_j}{2} - \mathcal{D}(U_i,U_j) \frac{U_j-U_i}{2}\right)\vec{n}_{ij},\label{eq:flux vfroe}
\f}
- Using the \ref vffc local linearisation of the fluxes, we obtain the following formula:
\f{eqnarray*}
\overrightarrow{\Phi}^{conv, VFFC}_{ij}&= &\frac{\mathcal{F}^{conv}(U_i) + \mathcal{F}^{conv}(U_j)}{2} \vec{n}_{ij}- \mathcal{D}(U_i,U_j) \frac{\mathcal{F}^{conv}(U_j)-\mathcal{F}^{conv}(U_i)}{2} \vec{n}_{ij},\label{eq:flux vffc}
\f}
where 
- \f$\mathcal{D}\f$ is an upwinding matrix,
- \f$A(U_i,U_j)\f$ the Roe matrix
- and \f$A^-= \frac{A - \mathcal{D}}{2}\f$.

 The choice \f$\mathcal{D}= 0\f$ gives the \ref centered upwinding, \f$\mathcal{D}= |A|\f$ for the \ref roe formulation and  \f$\mathcal{D}= sign(A)\f$ for the \ref vfroe and \ref vffc formulations give the full \ref upwind upwinding. The \ref lowMach, \ref pressureCorrection and \ref staggered upwindings allow more precision for  Mach number flows. 

The diffusion numerical flux \f$\overrightarrow\Phi_{ij}^{diff}\f$ is approximated on structured meshes using the formula:
\f{eqnarray*}
\overrightarrow \Phi_{ij}^{diff}= D (\frac{U_i+U_j}{2},\vec{n}_{ij})(U_j-U_i),\label{eq:flux diff}
\f}

where 
-\f$D(U,\vec{n}_{ij})=\nabla\mathcal{F}^{diff}(U).\vec{n}_{ij}\f$ is the matrix of the diffusion tensor.
\f$(\ref{eq:flux diff})\f$ is not accurate for highly non structured or non conforming meshes. However, since we are mainly interested in convection driven flows, we do not ask for a very precise scheme.

Finally, since \f$\sum_{j \in N(i)}\mathcal {F}^{conv}(U_i). \vec{n}_{ij}=0\f$, using \f$(\ref{eq:flux roe})\f$ and \f$(\ref{eq:flux diff})\f$ the equation \f$(\ref{eq:numer scheme})\f$ of the semi-discrete scheme becomes:
\f{eqnarray*}
\frac{\mathrm{d} U_{i}}{\mathrm{d} t} + \sum_{j\in N(i)} {\frac{s_{ij}}{v_i}\{(A^-+ D)(U_i,U_j)\}(U_j-U_i)} = S_i(U,x),\label{eq:reduced scheme} 
\f}

The source term in \f$(\ref{eq:reduced scheme})\f$ can be approximated using either a 
\f{eqnarray*}
 \textrm{ Centered source } S_i=S(U_i)\nonumber
\f}
or an
\f{eqnarray*}
 \textrm{ Upwind source } S_i&=&\frac{1}{2}(Id-signe(A^{Roe}_{i,i+1}))\frac{S(U_i)+S(U_{i+1})}{2}\\
			      &&+\frac{1}{2}(Id+signe(A^{Roe}_{i-1,i}))\frac{S(U_{i-1})+S(U_i)}{2}.\nonumber
\f}



\section ExplSchemes Explicit schemes


In explicit schemes, in order to compute the values \f$U_i^{n+1}\f$, the fluxes \f$\Phi^{conv}_{ij}\f$, \f$\Phi^{diff}_{ij}\f$ and the source term \f$S(U,x)\f$ in \f$(\ref{eq:numer scheme})\f$ are evaluated at time \f$n\f$ :
\f{eqnarray*}
\frac{U_{i}^{n+1} - U_{i}^{n}}{\Delta t} &+& \sum_{j\in N(i)} \frac{s_{ij}}{v_i}\left(\frac{1}{2}(\mathcal{F}^{conv}(U_i^n) + \mathcal{F}^{conv}(U_j^n)). \vec{n}_{ij}- \mathcal{D}(U_i^n,U_j^n,\vec{n}_{ij}) \frac{U_j^n-U_i^n}{2}\right)\\
&+&\frac{s_{ij}}{v_i} D (\frac{U_i+U_j}{2},\vec{n}_{ij})(U_j-U_i)= S(U^n,x_i), \nonumber
\f}
or equivalently using \f$(\ref{eq:flux roe})\f$ and \f$(\ref{eq:flux diff})\f$
\f{eqnarray*}
\label{explicitscheme}
\frac{U_{i}^{n+1} - U_{i}^{n}}{\Delta t} + \sum_{j\in N(i)} {\frac{s_{ij}}{v_i}\{(A^-+ D)(U_i^{n},U_j^{n},\vec{n}_{ij})\}(U^{n}_j-U^{n}_i)} =  S(U^n,x_i). 
\f}

From the system \f$(\ref{explicitscheme})\f$ we can obtain \f$U_i^{n+1}\f$ easily using matrix-vector products and vector sum. 
However the time step of explicit methods is constrained by the CFL condition for stability reasons.

\section ImplSchemes Implicit schemes


In implicit schemes, in order to compute the values \f$U_i^{n+1}\f$, the fluxes \f$\Phi^{conv}_{ij}\f$, \f$\Phi^{diff}_{ij}\f$ and the source term \f$S(U,x)\f$ in \f$(\ref{eq:numer scheme})\f$ are evaluated at time \f$n+1\f$ :
\f{eqnarray*}
\frac{U_{i}^{n+1} - U_{i}^{n}}{\Delta t} &+& \sum_{j\in N(i)} \frac{s_{ij}}{v_i}\left(\frac{1}{2}(\mathcal{F}^{conv}(U_i^{n+1}) + \mathcal{F}^{conv}(U_j^{n+1})). \vec{n}_{ij}- \mathcal{D}(U_i^{n+1},U_j^{n+1},\vec{n}_{ij}) \frac{U_j^{n+1}-U_i^{n+1}}{2}\right)\\
&+&\frac{s_{ij}}{v_i} D (\frac{U_i+U_j}{2},\vec{n}_{ij})(U_j-U_i) = S(U^{n+1},x_i),
\f}
or equivalently using \f$(\ref{eq:flux roe})\f$ and \f$(\ref{eq:flux diff})\f$
\f{eqnarray*}
\label{implicitscheme}
\frac{U_{i}^{n+1} - U_{i}^{n}}{\Delta t} + \sum_{j\in N(i)} {\frac{s_{ij}}{v_i}\{(A^-+ D)(U_i^{n+1},U_j^{n+1},\vec{n}_{ij})\}(U^{n+1}_j-U^{n+1}_i)} =  S(U^{n+1},x_i). 
\f}

The system \f$(\ref{implicitscheme})\f$ is nonlinear. The computation of \f$U_i^n\f$ is more expensive but we can expect to use larger time steps than would be allowed with the explicit scheme.

We use the following Newton iterative method to obtain the required solutions:
\f{eqnarray*}
\frac{\delta U_i^{k+1}}{\Delta t}& + &\sum_{j \in N(i)} \frac{s_{ij}}{v_i} \left[( A^-+ D)(U_i^k,U_j^k) \right] \left(\delta U_j^{k+1}- \delta U_i^{k+1} \right)\\
& = &- \frac{U^k_i-U^n_i}{\Delta t} - \sum_{j \in N(i)} \frac{s_{ij}}{v_i} \left[( A^-+ D)(U_i^k,U_j^k)\right] (U^k_j-U^k_i),
\f}
where :
- \f$\delta U_i^{k+1} = U_i^{k+1} - U_i^{k}\f$ is the variation of the \f$k\f$-th iterate that approximate the solution at time \f$n+1\f$. 

Defining the unknown vector \f$\mathcal U = (U_1,\dots,U_N)^t\f$, each Newton iteration for the computation of \f$\mathcal U\f$ at time step \f$n+1\f$ requires the numerical solution of the following linear system:
\f{eqnarray*}
 \label{eq: linearsystem}
 \mathcal A(\mathcal U^k)\delta \mathcal U^{k+1} =  b(\mathcal U^n, \mathcal U^{k}).
\f}



\section NumNSEqs Numerical scheme for the Navier-Stokes equations

For the Navier-Stokes equations \f$U=(\rho, \vec q, \rho E)^t\f$ and the fluxes in \f$(\ref{Navierstokes})\f$ write
\f{eqnarray*}
\mathcal{F}^{conv}(U)= \left(\begin{array}{c}
\vec{q} \\
\vec{q} \otimes \frac{\vec{q}}{\rho} + p {I}_d \\
\left( \rho E + p \right) \frac{\vec{q}}{\rho} 
\end{array}
 \right) ,\;
\mathcal{F}^{diff}(U)=\left(\begin{array}{c}
0\\
-\nu \vec \nabla (\frac{\vec{q}}{\rho}) \\
-\lambda \vec \nabla T 
\end{array}
 \right).
\f}

For the Euler equations, we can build the \ref roe matrix \f$A(U_i,U_j)\f$ explicitly using the Roe averaged state \f$U_{Roe}(U_i,U_j)=(\tilde\rho, \tilde \rho\tilde u, \tilde\rho \tilde E=\tilde \rho\tilde H-\tilde p)^t\f$ defined by
 \f{eqnarray*}
\tilde{\rho}&=&\sqrt{\rho_{i}\rho_{j}}\\
\tilde{u}&=&\frac{\sqrt{\rho_{i}}u_{i}+\sqrt{\rho_{j}}u_{j}}{\sqrt{\rho_{i}}+\sqrt{\rho_{j}}}\\
\tilde{H}&=&\frac{\sqrt{\rho_{i}}H_{i}+\sqrt{\rho_{j}}H_{j}}{\sqrt{\rho_{i}}+\sqrt{\rho_{j}}}.
\f}
The Roe matrix writes (see \ref leveque)
\f{eqnarray*}
A_{Roe}(U_i,U_j)=\nabla\mathcal{F}^{conv}(U_{Roe}(U_i,U_j))\vec{n}_{ij}=
\left(\begin{array}{ccc}
       0  & 1 & 0\\
       \tilde{\chi}+(\frac{1}{2}\tilde{\kappa} -1)\tilde{u}^2 & (2-\tilde{\kappa})\tilde{u} & \tilde{\kappa}\\
       (\tilde{\chi}+\frac{1}{2}\tilde{\kappa} \tilde{u}^2- \tilde{H})\tilde{u} & \tilde{H}-\tilde{\kappa} \tilde{u}^2 & (\tilde{\kappa}+1)\tilde{u}
      \end{array}\right)\label{eq:matrice de roe F4}
\f}

The diffusion numerical flux \f$\overrightarrow\Phi_{ij}^{diff}\f$ is approximated with the formula:
\f{eqnarray*}
\overrightarrow \Phi_{ij}^{diff}= D (\frac{U_i+U_j}{2})(U_j-U_i)
\f}
 with the matrix 
\f{eqnarray*}
D(U)= 
\left(
\begin{array}{ccc}
 0&\vec 0& 0\\
\frac{\nu \vec q}{\rho^2}& \frac{-\nu}{\rho} I_d&0\\
\frac{\lambda}{c_v}\left(\frac{c_vT}{\rho}-\frac{||\vec q||^2}{2\rho^3}\right)&\quad \frac{{\vec q}^{\;t} \lambda}{\rho^2 c_v}&\quad-\frac{\lambda}{c_v  \rho}
\end{array}
\right) , \f} 
where \f$c_v\f$ is the heat capacity at constant volume.



*/