[
  {
    "objectID": "index.html",
    "href": "index.html",
    "title": "CERF : Computation, ERosion and Flow V1.0",
    "section": "",
    "text": "Introduction\nThis document is a user guide for CERF. CERF is a finite-volume computational code dedicated to the simulation of free-surface flows. CERF is a research code developed at the Institut de Mathématiques de Toulon under the supervision of Frédéric Golay. CERF solves the Saint-Venant equations and the isothermal two-fluid Euler equations. These hyperbolic systems of equations allow the use of an explicit and parallel time integration scheme. The data structure of CERF is based on a BB-AMR (Block Based Adaptive Mesh Refinement) approach, which makes it possible to handle meshes that adapt in both time and space. CERF is a parallel computational code that uses the MPI library for communication between processes. CERF is written in Fortran 90. In its V1.0 version, the CERF mesh consists exclusively of hexahedra.\nThis document is organized as follows. In the first chapter, we present the physical models used in CERF to simulate free-surface flows. Erosion models are currently under development. We also describe the numerical schemes used. The second chapter is devoted to the installation and use of CERF, as well as the description of its input and output files. We explain how to configure a simulation and how to visualize the results. Finally, the third chapter consists of a set of documented examples.\n\n\n\n\n\n\nImportantWarning\n\n\n\nMy command of the English language is, shall we say… rough. At best, it’s understandable “Globish.” As a result, this documentation has been automatically translated. I apologize in advance if any inaccuracies—which I may not be able to detect—appear from time to time.\n\n\nThis document was written with Quarto.\nThis software is governed by the CeCILL-B license, subject to French law and complying with the principles of free software distribution.",
    "crumbs": [
      "Introduction"
    ]
  },
  {
    "objectID": "intro.html",
    "href": "intro.html",
    "title": "License",
    "section": "",
    "text": "Copyright Institut de Mathématiques de Toulon, (2025).\nfrederic.golay@univ-tln.fr\nThis software is a finite-volume computational code dedicated to the simulation of free-surface flows.\nThis software is governed by the CeCILL-B license, subject to French law and complying with the principles of free software distribution. You may use, modify and/or redistribute this program under the terms of the CeCILL-B license as distributed by CEA, CNRS and INRIA at the following URL “https://cecill.info/licences/Licence_CeCILL-B_V1-en.html”.\nAs a counterpart to the access to the source code and rights to copy, modify and redistribute granted by the license, users are provided only with a limited warranty. For the same reason, only limited liability is incurred by the software’s author, the holder of the economic rights, and the successive licensors.\nIn this respect, the user’s attention is drawn to the risks associated with loading, using, modifying and/or developing or reproducing the software by the user, given its free-software nature, which may make it complex to manipulate, and which therefore means that it is reserved for developers and experienced professionals having in-depth computer knowledge. Users are therefore encouraged to load and test the software’s suitability for their needs under conditions that ensure the security of their systems and/or data, and, more generally, to use and operate it under the same conditions of security.\nThe fact that you are presently reading this means that you have had knowledge of the [CeCILL-B] license and that you accept its terms.",
    "crumbs": [
      "License"
    ]
  },
  {
    "objectID": "chapitre1.html",
    "href": "chapitre1.html",
    "title": "1  Modeling",
    "section": "",
    "text": "1.1 The finite-volume method\nThis chapter describes the physical models used in CERF to simulate flows (without erosion in the version V1.0), as well as how these models are implemented in the code. It is divided into three sections. The first section presents the finite-volume method in the general case. The second section describes the Saint-Venant shallow-water flow model. The third section presents the Eulerian two-fluid flow model used in CERF.\nThe finite volume method is a numerical method used to solve partial differential equations. This method is used in CERF to solve the Saint-Venant equations and the two-fluid flow equations.",
    "crumbs": [
      "<span class='chapter-number'>1</span>  <span class='chapter-title'>Modeling</span>"
    ]
  },
  {
    "objectID": "chapitre1.html#sec-vf",
    "href": "chapitre1.html#sec-vf",
    "title": "1  Modeling",
    "section": "",
    "text": "1.1.1 Principle of the finite volume method\nIn this section, we briefly recall the finite volume approximation, in dimension \\(n_{dim}\\), of a system of conservation equations of the form: \\[ \\frac{\\partial \\overrightarrow{w}}{\\partial t} + \\nabla \\cdot \\overline{\\overline{F}}(\\overrightarrow{w}) = \\overrightarrow{S}(\\overrightarrow{w})  \\tag{1.1}\\] where \\(\\overrightarrow{w}(\\overrightarrow{x},t)\\) is the vector of conservative variables of size \\(n_{dof}\\) depending on space and time, \\(\\overline{\\overline{F}}\\) the flux matrix of size \\(n_{dof} \\times n_{dim}\\), \\(\\overrightarrow{S}\\) the vector of source terms of size \\(n_{dof}\\), and \\(\\nabla \\cdot\\) the divergence operator.\nThe system is solved in \\(\\mathbb{R}^{n_{dim}}\\). The computational domain \\(\\Omega \\subset \\mathbb{R}^{n_{dim}}\\) is divided into a set of control volumes, also called cells, \\(\\Omega = \\cup_{k} C_k\\). We begin by integrating Equation 1.1 over a cell \\(C_k\\) using Green’s formula:\n\\[\\int _{C_k} \\frac{\\partial \\overrightarrow{w}}{\\partial t} d \\Omega + \\int _{C_k} \\nabla \\cdot \\overline{\\overline{F}} d \\Omega = \\int _{C_k} \\overrightarrow{S}(\\overrightarrow{w}) d \\Omega\\]\n\\[\n\\int _{C_k} \\frac{\\partial \\overrightarrow{w}}{\\partial t} d \\Omega + \\sum _{a} \\int _{\\partial C_{k/a}} \\overline{\\overline{F}} \\cdot  \\overrightarrow{n}_{k/a} d \\partial\\Omega = \\int _{C_k} \\overrightarrow{S}(\\overrightarrow{w}) d \\Omega \\tag{1.2}\\]\nwhere \\(\\overrightarrow{n}_{k/a}\\) denotes the unit normal vector on the boundary \\(\\partial C_{k/a}\\) between cell \\(C_k\\) and the neighboring cell \\(C_a\\). The summation is carried out over all the neighboring cells of \\(C_k\\).\nThe finite volume method then consists in approximating the solution by piecewise constant functions per cell, such that by setting \\[\\overrightarrow{w}_k(t) \\approx \\frac {1}{|C_k|}\\int _{C_k} \\overrightarrow{w}(\\overrightarrow{x},t) d \\Omega\\] where \\(|C_k|\\) represents the volume of cell \\(C_k\\), we approximate Equation 1.2 by \\[\n|C_k| \\frac{\\partial \\overrightarrow{w}_k(t)}{\\partial t} +  \\sum _{a}  |\\partial C_{k/a}| \\: \\widehat{\\overrightarrow{F}} \\left(\\overrightarrow{w}_k(t),\\overrightarrow{w}_a(t),\\overrightarrow{n}_{k/a} \\right) = |C_k|S_k(t) \\tag{1.3}\\]\nwhere \\(S_k(t)\\approx \\frac {1}{|C_k|}\\int _{C_k}S(\\overrightarrow{x},t) d \\Omega\\) represents the average value of the source term vector over the cell, \\(|\\partial C_{k/a}|\\) the surface area of the face between cells \\(C_k\\) and \\(C_a\\), and \\(\\widehat{\\overrightarrow{F}}\\) an approximation of the flux.\nSince the unknowns \\(\\overrightarrow{w}\\) are not defined at the interface between cells, the flux cannot be computed directly. It is therefore approximated by a numerical flux depending on the interface and on the values of the variables in the adjacent cells. We denote:\n\\[\\overrightarrow{F}_{k/a}(t) = \\widehat{\\overrightarrow{F}} \\left(\\overrightarrow{w}_k(t),\\overrightarrow{w}_a(t),\\overrightarrow{n}_{k/a} \\right) \\approx \\frac{1}{|\\partial C_{k/a}|}\\int _{\\partial C_{k/a}}  \\overline{\\overline{F}} \\cdot \\overrightarrow{n}_{k/a} ds\\]\nThere are many possible approaches: Rusanov, HLLC, etc. In CERF, the numerical flux is computed using Godunov’s method, based on the exact solution of the Riemann problem at the interface \\(k/a\\) (for more details, see Section 1.2.5 and Section 1.3.5).\nSecond-order accuracy in space can be achieved using the MUSCL (Monotone Upstream Scheme for Conservation Laws) method to reconstruct the primitive variables at the interfaces. From the values of the variables on the neighboring cells, the gradient of the primitive variables is estimated. Let \\({w_p}_i\\) be the \\(i\\)-th component of the primitive variables and \\({w_p}_i^k\\) its value at the center of the cell; the gradient over cell \\(C_k\\) can be estimated as: \\[ |C_k| \\widetilde{\\nabla} {w_p}_i = \\int _{C_k} \\nabla {w_p}_i d\\Omega =\\sum _{a} \\int _{\\partial C_{k/a}} w_{p_i} \\overrightarrow{n}_{k/a} d\\partial\\Omega =\\sum_{a} |\\partial C_{k/a}| \\frac{1}{2}\\left( {w_p}_i^k + {w_p}_i^a  \\right) \\overrightarrow{n}_{k/a}\\] and the primitive variables of cell \\(C_k\\) with center \\(K\\) can then be reconstructed at the centers \\(F\\) of the cell faces such that: \\[w_{p_i}^F = w_{p_i}^k + \\widetilde\\nabla w_{p_i} \\cdot \\overrightarrow{KF}\\]\nTo achieve second-order accuracy, \\(\\widehat{\\overrightarrow{F}} \\left(\\overrightarrow{w}_k(t),\\overrightarrow{w}_a(t),\\overrightarrow{n}_{k/a} \\right)\\) is then replaced by \\(\\widehat{\\overrightarrow{F}} \\left(\\overrightarrow{w}_k^F(t),\\overrightarrow{w}_a^F(t),\\overrightarrow{n}_{k/a} \\right)\\)\n\n\n\n\n\n\nNote\n\n\n\nIn general, a limiter must be combined with the MUSCL method. In CERF, the Barth-Jespersen limiter is used: \\(Min(w_{p_F},w_{p_F}) \\leq  w_{p_i}^F \\leq Max(w_{p_F},w_{p_F})\\).\n\n\n\n\n\n1.1.2 Time discretization\nTo facilitate the parallelization of the code, CERF uses only explicit schemes: the first- or second-order Runge-Kutta scheme. The time interval \\(\\left[0,T\\right]\\) is discretized such that \\(t_{n+1} = t_n + \\delta  t_n\\), where \\(\\delta t_n\\) represents the time step at time \\(t_n\\). The first-order Runge-Kutta method, also known as the first-order explicit Euler method, consists in rewriting Equation 1.3 at time \\(t_n\\) by approximating the time derivative as a finite difference, such that: \\[\n|C_k|\\frac{ \\overrightarrow{w}_k(t_{n+1}) - \\overrightarrow{w}_k(t_{n})}{\\delta t_n} +  \\sum _{a}  |\\partial C_{k/a}| \\: \\overrightarrow{F}_{k/a}(t_n) = |C_k|\\overrightarrow{S}_k(t_n)\\]\nBy denoting \\(\\overrightarrow{w}_k(t_{n})=\\overrightarrow{w}_k^{n}\\), the new value of the conservative variables at time \\(t_{n+1}\\) is thus computed explicitly: \\[\n\\overrightarrow{w}_k^{n+1} = \\overrightarrow{w}_k^{n}+\\frac{\\delta t_n}{|C_k|}  \\sum _{a}  |\\partial C_{k/a}| \\: \\overrightarrow{F}_{k/a}^n - \\delta t_n \\overrightarrow{S}_k^n\\]\nThe second-order Runge-Kutta method consists in computing the solution at time \\(t_{n+1}\\) in two steps. The first step consists in computing the solution at time \\(t_{n+\\frac{1}{2}}\\) using the first-order explicit Euler scheme. The second step consists in computing the solution at time \\(t_{n+1}\\) using the first-order explicit Euler scheme together with the solution obtained at time \\(t_{n+\\frac{1}{2}}\\). \\[  \\overrightarrow{w}_k^{n+\\frac{1}{2}} = \\overrightarrow{w}_k^{n}+\\frac{\\delta t_n}{2 |C_k|}   \\sum _{a}  |\\partial C_{k/a}| \\: \\overrightarrow{F}_{k/a}^n - \\frac{\\delta t_n}{2} \\overrightarrow{S}_k^n \\] \\[  \\overrightarrow{w}_k^{n+1} = \\overrightarrow{w}_k^{n}+\\frac{\\delta t_n}{|C_k|}   \\sum _{a}  |\\partial C_{k/a}| \\:  \\overrightarrow{F}_{k/a}^{n+\\frac{1}{2}} - \\delta t_n\\overrightarrow{S}_k^{n+\\frac{1}{2}} \\]\n\n\n\n\n\n\nImportant 1.1: Courant-Friedrichs-Lewy stability condition\n\n\n\nExplicit methods are simple to implement, but their stability is constrained by the Courant-Friedrichs-Lewy (CFL) condition. The CFL condition is a necessary condition to guarantee the stability of the numerical scheme. This condition depends on the wave propagation speed in the medium and on the size of the cells. In practice, it is determined for each hyperbolic model considered. In CERF, the cell size \\(d_k\\) is determined by the ratio between the volume of the cell and the surface area of its boundary: \\[ d_k = \\frac{|C_k|}{\\sum_{a}|\\partial C_{k/a}|}\\]\nIn the 1D or 2D case, symmetry faces are not taken into account when computing the boundary surface area. The CFL condition is then given by: \\[ \\delta t_n \\leq \\alpha_{CFL}  \\frac{\\min_{k}d_k}{\\max_{k} \\lambda_k}\\]\nwhere \\(\\alpha_{CFL}\\) is a coefficient to be defined by the user and \\(\\lambda_k\\) is the wave propagation speed in cell \\(C_k\\), depending on the model used (Note 1.1 for Saint-Venant and Note 1.2 for two-fluid flow).\n\n\n\n\n\n1.1.3 Spatial discretization\n\n\n\n\n\n\nImportant\n\n\n\nIn its V1.0 version, CERF uses exclusively hexahedral cells and therefore quadrangular faces! Generalization is currently under development.\n\n\n\n\n\n\n\n\nNote\n\n\n\nCERF allows the use of dynamic meshes, meaning that the mesh can be regularly refined or coarsened (AMR: Adaptive Mesh Refinement).\n\n\nCERF is a natively three-dimensional code. In its initial version, CERF uses only hexahedral cells and therefore quadrangular faces. Cells and faces are not necessarily regular. Nevertheless, as illustrated in Figure 1.1, two-dimensional simulations can be performed by creating a mesh with a single cell in the \\(z\\) direction and imposing “mirror”-type conditions on the faces in the \\(x,y\\) plane. Similarly, one-dimensional simulations can be performed by creating a mesh with a single cell in the \\(y\\) and \\(z\\) directions and imposing “mirror”-type conditions on the faces in the \\(x,z\\) and \\(x,y\\) planes.\n\n\n\n\n\n\n\nFigure 1.1: Three-, two-, and one-dimensional meshes with hexahedral cells\n\n\n\n\n\n\n1.1.3.1 Meshing strategy\nCERF provides mesh refinement and coarsening tools (AMR: Adaptive Mesh Refinement). To avoid penalizing computation time, the mesh is not modified at every time step, but at instants determined by the user. This approach makes it possible to maintain good mesh quality while limiting computational cost. This method is referred to as BB-AMR (Block Based Adaptive Mesh Refinement). The various steps of the meshing strategy used in CERF are described below.\n\n\n\n\n\n\n\nFigure 1.2: Initial conforming mesh: the “Master Mesh”.\n\n\n\n\nStep 1: For the sake of simplicity, the mesh shown in Figure 1.2 consists of regular quadrangles, which can represent a hexahedral mesh viewed along \\(z\\) with a single cell in the \\(z\\) direction. As illustrated in Figure 1.2, we start from an initial conforming mesh made up solely of hexahedra of arbitrary shape. In a simple case, the mesh can be created by CERF; otherwise it can be imported in GMSH V2.2 format. This mesh is referred to as the master mesh. It is essential that the mesh be conforming in order to easily define the faces connecting the cells. The cells of this initial mesh are referred to as blocks.\n\n\n\n\n\n\n\nFigure 1.3: Numbering of the initial mesh according to Morton codes.\n\n\n\n\nStep 2: The blocks of the initial mesh are ordered by numbering them according to the Z-order or Morton code (Figure 1.3). With this approach, known as a Space Filling Curve (Section 1.1.3.2), the numbering makes it possible to traverse the blocks, even in three dimensions, following a numbering such that successive blocks in the numbering are, most of the time, neighbors in space.\n\n\n\n\n\n\n\nFigure 1.4: Definition of the block refinement level.\n\n\n\n\nStep 3: The refinement level of the blocks is then defined, as illustrated in Figure 1.4. Initially, the refinement level is defined by the user. During the computation, the refinement level is defined by a criterion Section 1.1.3.3. As verified in Ersoy, Golay, and Yushchenko (2013) and Thomas Altazin and Yushchenko (2016), the difference in refinement level between two blocks must not exceed \\(1\\). During refinement, the block is divided by \\(2\\) in each direction, for each level required. Thus, denoting by \\(Niv\\) the required refinement level (where \\(Niv=0\\) represents a block that is not refined) and \\(dim=1,2,3\\) the dimension of the problem, the number of cells in a block will be: \\(2^{dim \\times Niv}\\). In the example shown in Figure 1.4, block \\(1\\) is not refined, block \\(2\\) is refined by \\(1\\) level, and block \\(4\\) by \\(2\\) levels. This leads, as illustrated in Figure 1.5 for the case \\(n_{dim}=2\\), to \\(1\\) cell in block \\(1\\), \\(4\\) cells in block \\(2\\), and \\(16\\) cells in block \\(4\\).\n\n\n\n\n\n\n\nFigure 1.5: Decomposition of the mesh into domains.\n\n\n\n\nStep 4: Since the time integration schemes in CERF are explicit, it is straightforward to parallelize the solution algorithm. To do this, the mesh must be split into several domains, ideally with a balanced number of cells and minimal interfaces, in order to minimize data exchange between domains. Many algorithms exist, some of them quite complex, that perform this operation efficiently. In CERF, we chose a compromise between simplicity and efficiency. We opted for a decomposition driven by Morton codes. Knowing the number of cells per block, the blocks are traversed following the Morton numbering and assigned to a domain. The number of cells per domain is kept balanced. As illustrated in Figure 1.5, the initial mesh is split into 3 domains made up of 25, 21, and 18 cells. Each computational domain is then assigned to a computing process.\n\n\n\n\n\n\n\nFigure 1.6: Mesh generation.\n\n\n\n\nStep 5: The final step consists in generating the mesh. To exchange information between domains, CERF creates ghost cells. These cells are copies of the neighboring cells of the domain. They are used to compute the numerical fluxes at the interfaces between domains. As illustrated in Figure 1.6, ghost cells are shown in gray. Note that the mesh generated in this way is not conforming! For example, the first cell at the bottom left of domain \\(1\\) in Figure 1.6 is surrounded by 6 faces (in the \\(x,y\\) plane). The refinement of the blocks defining the mesh naturally depends on the dimension of the problem being treated (Figure 1.7). Finally, if the mesh is modified during the computation, it is necessary to project the conservative variables onto the new mesh. In a remeshing operation, the level can only increase or decrease by at most \\(1\\), so that the projection step is relatively simple.\n\n\n\n\n\n\n\nFigure 1.7: Refinement depending on the dimension of the problem.\n\n\n\n\n\n\n\n1.1.3.2 Morton code\nSpace-filling curves are curves that make it possible to traverse a higher-dimensional space. First introduced by Peano, then by Hilbert or Lebesgue, they are used for example to index data, and in particular in the construction of quadtrees or octrees. In CERF, numbering by Morton code, also known as the Z-order curve, is used. This numbering makes it possible to index a set of data (points, cells, faces, etc.) such that data that are neighbors in space are also neighbors in the indexing. This property is used in CERF to split the mesh into computational domains by indexing the blocks of the master mesh. Consider a point with coordinates \\((x,y,z)\\) in a three-dimensional space. To index this point, the coordinates are first converted into positive integers. To do this, the smallest coordinates of the study domain \\((x_{min},y_{min},z_{min})\\) are defined, along with a distance \\(d_{morton}\\) defining the grid of the domain, for example one tenth of the size of the smallest cell in the domain. The integer coordinates of the point \\((ix,iy,iz)\\) are then defined by: \\[ix =  \\frac{x-x_{min}}{d_{morton}} ,\\: iy =  \\frac{y-y_{min}}{d_{morton}} ,\\: iz =  \\frac{z-z_{min}}{d_{morton}} \\] Morton numbering then consists in converting these integer coordinates into base \\(2\\) \\((ix_2,iy_2,iz_2)\\), then interleaving the bits of each coordinate to form a single base-\\(2\\) integer \\(M_2\\). Finally, converting back to base \\(10\\), the Morton index \\(M\\) of the point \\((x,y,z)\\) is obtained.\n\nExample 1.1 \\(ix=2\\), \\(iy=5\\), \\(iz=0\\).\nWe then have \\(ix_2=010\\), \\(iy_2=101\\), and \\(iz_2=000\\).\nInterleaving the bits gives \\(M_2=010100010\\), and in base \\(10\\), \\(M=162\\).\n\n\n\n\n\n\n\nNote\n\n\n\nAll of CERF’s indexing tools can be found in the module ~/sources/UTI/mod_zorder.f90.\n\n\n\n\n\n1.1.3.3 Refinement criterion\nIn general, a refinement criterion is defined for each physical model considered. This criterion is often derived from physical considerations based, for example, on the gradient of one of the primitive variables, or, more relevantly and more elaborately, on an a posteriori error estimate or the construction of an error indicator. In the case of hyperbolic problems, there exists a more general criterion derived from the work of Croisille (1990). This criterion is based on the numerical production of entropy. It was used by Golay (2009), Ersoy, Golay, and Yushchenko (2013), and studied at length by Puppo (2004). The entropy inequality of Lax (1971) allows us to state that for any hyperbolic system of conservation laws of the form Equation 1.1, there exists a mathematical entropy \\(s(\\overrightarrow{w})\\) and an entropy flux \\(\\overrightarrow{\\psi}_s(\\overrightarrow{w})\\) such that: \\[\\frac{\\partial s(\\overrightarrow{w})}{\\partial t} + \\nabla \\cdot \\overrightarrow{\\psi}_s(\\overrightarrow{w}) \\leq 0 \\tag{1.4}\\] where the entropy flux \\(\\overrightarrow{\\psi}_s(\\overrightarrow{w})\\) must satisfy a compatibility equation \\(\\nabla_w \\psi_s = \\nabla_w \\overline{\\overline{F}} \\cdot \\nabla_w s\\). For smooth solutions, Equation 1.4 is an equality. In the case of a shock, it is a strict inequality. However, in numerical simulations, the numerical production of entropy is observed not to be zero! It turns out that the numerical production of entropy is a very good error indicator, except at shocks, where it tends to infinity. We therefore define \\(p^s_k\\), the numerical entropy production over cell \\(C_k\\), as:\n\\[p^s_k = \\frac{1}{|C_k|} \\int_{C_k} \\left( \\frac{\\partial s(\\overrightarrow{w})}{\\partial t} + \\nabla \\cdot \\overrightarrow{\\psi}_s(\\overrightarrow{w}) \\right) d\\Omega\\].\nThe refinement criterion \\(C_b\\) is defined per block as the ratio of the numerical entropy production in the block to the total numerical entropy production over the domain. \\[ C_b = \\frac{\\int_{Bloc} \\left( \\frac{\\partial s(\\overrightarrow{w})}{\\partial t} + \\nabla \\cdot \\overrightarrow{\\psi}_s(\\overrightarrow{w}) \\right) d\\Omega}{\\int_{\\Omega} \\left( \\frac{\\partial s(\\overrightarrow{w})}{\\partial t} + \\nabla \\cdot \\overrightarrow{\\psi}_s(\\overrightarrow{w}) \\right) d\\Omega}\\] In CERF, the refinement criterion is defined per block. If the criterion exceeds a threshold value \\(\\alpha_{max}\\), the block is refined. If the criterion is below a threshold value \\(\\alpha_{min}\\), the block is coarsened. The threshold values are defined by the user.",
    "crumbs": [
      "<span class='chapter-number'>1</span>  <span class='chapter-title'>Modeling</span>"
    ]
  },
  {
    "objectID": "chapitre1.html#sec-sv",
    "href": "chapitre1.html#sec-sv",
    "title": "1  Modeling",
    "section": "1.2 Saint-Venant model",
    "text": "1.2 Saint-Venant model\n\n1.2.1 Governing equations\nThe shallow-water flow model, or Saint-Venant model, is a simplified flow model. It is derived from the Navier-Stokes equations by neglecting the viscosity terms, assuming a hydrostatic pressure distribution, and assuming that the water depth is much smaller than the characteristic length of the flow. The original Barré de Saint-Venant model is a 1D model Saint-Venant (1871); the derivation of the model is explained, for example, by Marche (2007) or Poussel (2024). Development of the model in CERF was initiated by Pons (2018). In a frame of reference \\((O,x,y,z)\\), we consider a shallow-water flow over a domain where the bathymetry is defined by \\(z_b(x,y)\\), with a water depth \\(h(x,y)\\), so that the free-surface level is defined by \\(\\zeta=h+z_b\\), as illustrated in Figure 1.8.\n\n\n\n\n\n\n\nFigure 1.8: Notation of the shallow-water flow model\n\n\n\n\nThe Saint-Venant equations are as follows: \\[ \\frac{\\partial h}{\\partial t} + \\nabla \\cdot (h \\overrightarrow{u}) = 0  \\tag{1.5}\\] \\[ \\frac{\\partial h\\overrightarrow{u}}{\\partial t} + \\nabla \\cdot (h \\overrightarrow{u} \\otimes \\overrightarrow{u} + g\\frac{h^2}{2} \\overline{\\overline{I}} ) = -gh \\nabla z_b -gh\\overrightarrow{S_f}  \\tag{1.6}\\]\nwhere \\(\\overrightarrow{u}=&lt;u,v&gt;^T\\) represents the flow velocity, \\(g\\) the acceleration due to gravity, \\(S_f\\) the friction term, and where \\(\\overline{\\overline{I}}\\) denotes the identity matrix, \\(\\nabla\\) the gradient operator, \\(\\nabla \\cdot\\) the divergence operator, and \\(\\otimes\\) the tensor product. We define the conservative variables \\(\\overrightarrow{w}=&lt;h,hu,hv&gt;^T\\) and the primitive variables \\(\\overrightarrow{w_p}=&lt;h,u,v&gt;^T\\). This corresponds to the general form Equation 1.1 with: \\[\\overline{\\overline{F}}(\\overrightarrow{w}) = \\left[  \\begin{array}{cc} hu & hv \\\\ hu^2 + g\\frac{h^2}{2} & h u v \\\\ h u v &    hv^2 + g\\frac{h^2}{2}   \\end{array}\\right]\\]\nand\n\\[\\overrightarrow{S} = -gh \\nabla z_b -gh\\overrightarrow{S_f} \\]\n\n\n\n1.2.2 Handling bathymetry\nHandling bathymetry in the finite-volume resolution scheme requires particular care. Indeed, the bathymetry jump inherent to the fact that bathymetry is considered constant per cell can generate spurious waves. The reconstruction method proposed by Audusse et al. (2004) is therefore implemented in CERF. Let \\(\\Delta z_{k/a}\\) be the bathymetry jump at the interface between cells. At the interface, we then distinguish the flux \\(\\overrightarrow{F}_k\\) contributing to cell \\(C_k\\) and the flux \\(\\overrightarrow{F}_a\\) contributing to cell \\(C_a\\), such that:\n\\[\\overrightarrow{F}_k \\left( \\overrightarrow{w}_k(t),\\overrightarrow{w}_a(t),\\overrightarrow{n}_{k/a}, \\Delta z_{k/a}\\right)=\\overrightarrow{F} \\left(\\overrightarrow{w}_k^*(t),\\overrightarrow{w}_a^*(t),\\overrightarrow{n}_{k/a} \\right) + \\left\\{ \\begin{array}{cc} 0 \\\\  \\frac{g}{2} \\left(h_k^2 - {h_k^*}^2 \\right)     \\overrightarrow{n}_{k/a} \\end{array} \\right\\} \\] \\[\\overrightarrow{F}_a \\left( \\overrightarrow{w}_k(t),\\overrightarrow{w}_a(t),\\overrightarrow{n}_{k/a}, \\Delta z_{k/a}\\right)=\\overrightarrow{F} \\left(\\overrightarrow{w}_k^*(t),\\overrightarrow{w}_a^*(t),\\overrightarrow{n}_{k/a} \\right) + \\left\\{ \\begin{array}{cc} 0 \\\\  \\frac{g}{2} \\left(h_a^2 - {h_a^*}^2 \\right)     \\overrightarrow{n}_{k/a} \\end{array} \\right\\} \\] where \\(\\overrightarrow{w}_k^*=&lt;h_k^*,\\overrightarrow{u}_k&gt;^T\\) and \\(\\overrightarrow{w}_a^*=&lt;h_a^*,\\overrightarrow{u}_a&gt;^T\\), with: \\[h_k^* = max \\left(0, h_k - max \\left( 0, \\Delta z_{k/a}\\right) \\right)\\] \\[h_a^* = max \\left(0, h_a - max \\left( 0, \\Delta z_{k/a}\\right) \\right)\\]\n\n\n\n1.2.3 Refinement criterion\nFor the Saint-Venant model, to define the numerical entropy production Equation 1.4, we use the entropy defined by: \\[ s(\\overrightarrow{w}) = \\frac{1}{2} \\left( h ||\\overrightarrow{u}||^2 + gh^2 + 2 g h z_b \\right) \\] and the entropy flux \\[ \\overrightarrow{\\psi}_s(\\overrightarrow{w}) = \\left( s + \\frac{g h^2}{2} \\right) \\overrightarrow{u} .\\]\n\n\n\n1.2.4 Handling friction\nIn the momentum conservation equation Equation 1.6, the friction term is taken into account. In CERF, the friction term is modeled using either the Manning-Strickler term or the Darcy-Weisbach term. The Manning-Strickler term is defined by: \\[ \\overrightarrow{S_f} = C_f  \\frac{||\\overrightarrow{u}|| \\overrightarrow{u}}{h^{4/3}}\\] where \\(C_f = n^2\\) is the Manning-Strickler friction coefficient.\nThe Darcy-Weisbach term is defined by: \\[ \\overrightarrow{S_f} = C_f  \\frac{||\\overrightarrow{u}|| \\overrightarrow{u}}{h}\\]\nThe friction term is handled using a “splitting” method. We first solve problem Equation 1.1 without the friction term, \\[ \\frac{\\partial \\overrightarrow{w}}{\\partial t} + \\nabla \\cdot \\overline{\\overline{F}}(\\overrightarrow{w}) = -gh \\nabla z_b  \\tag{1.7}\\]\nthen we solve \\[ \\frac{\\partial \\overrightarrow{w}}{\\partial t}  = -gh\\overrightarrow{S_f}  \\tag{1.8}\\]\nWe then assume that friction modifies neither the water depth nor the direction of the flow, but only the magnitude of the flow velocity. The friction problem for the velocity can then be solved using a first-order implicit scheme over a time step \\(\\delta t_n\\) and the velocity “penalized” as in Varra et al. (2024).\nLet \\(\\overrightarrow{u}^*\\) be the solution of problem Equation 1.7; problem Equation 1.8 can then be solved and the velocity penalized by: \\[ \\overrightarrow{u} = \\frac{2\\overrightarrow{u}^*}{1+\\sqrt{1 + 4\\beta}}  \\]\nwhere\n\\[\\beta = \\delta t_n g C_f h^{-q}\\]\nwith \\(q=4/3\\) for the Manning-Strickler model and \\(q=1\\) for the Darcy-Weisbach model.\n\n\n\n1.2.5 Riemann solver for the Saint-Venant model\n\n1.2.5.1 Overview\nNumerical fluxes in CERF are computed using a Riemann solver. There are many references concerning this solver, such as Toro (1997). This section describes the Riemann solver used for the Saint-Venant model, based on the course by Gloria Faccanoni.\nWe therefore seek to solve the Saint-Venant problem in one dimension. We consider a shallow-water flow with depth \\(h(x,t)\\) and horizontal velocity \\(u(x,t)\\) over a flat bottom, described by the following problem: \\[\n\\left \\{\n\\begin{array}{ll}\n\\displaystyle\\frac{\\partial h}{\\partial t} + \\frac{\\partial (hu)}{\\partial x} = 0 \\\\\n\\displaystyle\\frac{\\partial (hu)}{\\partial t} + \\frac{\\partial (hu^2 + \\frac{1}{2}gh^2)}{\\partial x} = 0\n\\end{array}\n\\right.\n\\tag{1.9}\\] With initial conditions such that: \\[\n\\overrightarrow{w_p}(x,0) =\\left\\{\n\\begin{array}{ll}\n\\overrightarrow{w_p}_L &\\mbox{ if } x&lt;0 \\\\\n\\overrightarrow{w_p}_R &\\mbox{ if } x&gt;0\n\\end{array}\n\\right.\n\\] where \\(\\overrightarrow{w_p}_L=&lt;h_L,u_L&gt;^T\\) and \\(\\overrightarrow{w_p}_R=&lt;h_R,u_R&gt;^T\\) are two constant states.\nBy expanding Equation 1.9, the system can be rewritten in the following form: \\[\n\\frac{\\partial}{\\partial t} \\left\\{ \\begin{array}{c} h \\\\ u \\end{array} \\right\\} + \\left[ \\begin{array}{cc} u & h \\\\ g & u \\end{array} \\right] \\frac{\\partial}{\\partial t} \\left\\{ \\begin{array}{c} h \\\\ u \\end{array} \\right\\} = \\left\\{ \\begin{array}{c} 0 \\\\ 0 \\end{array} \\right\\}\n\\]\nWe then seek the eigenvalues \\(\\lambda_1\\) and \\(\\lambda_2\\) as well as the eigenvectors \\(\\overrightarrow{r}_1\\) and \\(\\overrightarrow{r}_2\\) of the matrix thus obtained. A quick calculation gives: \\[\n\\lambda_1=u-\\sqrt{gh} , \\lambda_2=u+\\sqrt{gh}\n\\] with \\(\\lambda_1 \\leq \\lambda_2\\) and \\[\\overrightarrow{r}_1=\\left\\{ \\begin{array}{c} -\\sqrt{h} \\\\ \\sqrt{g} \\end{array} \\right\\} , \\overrightarrow{r}_2=\\left\\{ \\begin{array}{c} \\sqrt{h} \\\\ \\sqrt{g} \\end{array} \\right\\}\\].\n\n\n\n\n\n\nNote 1.1\n\n\n\nThe eigenvalues are used to determine the wave propagation speed within the cell Important 1.1 as: \\(\\lambda_k = Max\\left( |u-\\sqrt{gh}| , |u+\\sqrt{gh}| \\right)\\)\n\n\nTo determine the nature of the two characteristic fields, we compute \\(\\nabla_{w_p}\\lambda_1 \\cdot \\overrightarrow{r}_1 = 3/2 \\neq 0\\) and \\(\\nabla_{w_p}\\lambda_2 \\cdot \\overrightarrow{r}_2 = 3/2 \\neq 0\\). The characteristic fields are therefore genuinely nonlinear. There will therefore be no contact discontinuity, and the two waves are either shocks or rarefaction waves.\nSince there are \\(2\\) waves, \\(1\\) Riemann invariant \\(I_{1;1}\\), \\(I_{2;1}\\) can be defined for each wave, such that \\(\\nabla_{w_p}I_{1;1} \\cdot \\overrightarrow{r}_1 = 0\\) and \\(\\nabla_{w_p}I_{2;1} \\cdot \\overrightarrow{r}_2 = 0\\). It is easy to verify that \\(I_{1;1}=u+2\\sqrt{gh}\\) and \\(I_{2;1}=u-2\\sqrt{gh}\\) are Riemann invariants.\nThe entropy weak solution of the Riemann problem consists of at most three constant states \\(\\overrightarrow{w_p}_L\\), \\(\\overrightarrow{w_p}_*\\), \\(\\overrightarrow{w_p}_R\\), separated by two discontinuities. These discontinuities are either shocks or rarefaction waves. In what follows, we determine the state \\(\\overrightarrow{w_p}_*\\).\n\n\n1.2.5.2 Study of the rarefaction waves\nConsider a state \\(\\overrightarrow{w_p}_g\\) and a state \\(\\overrightarrow{w_p}_d\\), connected by a k-rarefaction wave. The following conditions must then be satisfied:\n\n\\(\\lambda_k(\\overrightarrow{w_p}_g)&lt;\\lambda_k(\\overrightarrow{w_p}_d)\\)\nThe Riemann invariant is conserved along the rarefaction wave \\(I_{k;1}(\\overrightarrow{w_p}_g)=I_{k;1}(\\overrightarrow{w_p}_d)\\)\n\nLet us apply these conditions to the Saint-Venant problem.\nConsider first a 1-rarefaction connecting a state \\(\\overrightarrow{w_p}_L\\) to a state \\(\\overrightarrow{w_p}_*\\), separated by a rarefaction wave. The Riemann invariants give \\[ I_{k;1}(\\overrightarrow{w_p}_L)=I_{k;1}(\\overrightarrow{w_p}_*) \\rightarrow u_L+2\\sqrt{g\\ h_L}=u_*+2\\sqrt{g\\ h_*} \\] that is \\[ u_* = u_L + 2\\sqrt{g} \\left( \\sqrt{h_L} - \\sqrt{h_*}\\right)  \\tag{1.10}\\] Furthermore, \\[\\lambda_k(\\overrightarrow{w_p}_L)&lt;\\lambda_k(\\overrightarrow{w_p}_*) \\rightarrow u_L-\\sqrt{g\\ h_L} &lt; u_*-\\sqrt{g\\ h_*} \\rightarrow h_* &lt; h_L\\] and \\[ u_* &gt; u_L\\]\nConsider finally a 2-rarefaction connecting a state \\(\\overrightarrow{w_p}_*\\) to a state \\(\\overrightarrow{w_p}_R\\), separated by a rarefaction wave. The Riemann invariants give \\[ I_{k;1}(\\overrightarrow{w_p}_*)=I_{k;1}(\\overrightarrow{w_p}_R) \\rightarrow u_*-2\\sqrt{g\\ h_*}=u_R-2\\sqrt{g\\ h_R} \\] that is \\[ u_* = u_R + 2\\sqrt{g} \\left(\\sqrt{h_*} - \\sqrt{h_R}\\right)  \\tag{1.11}\\] Furthermore, \\[\\lambda_k(\\overrightarrow{w_p}_*)&lt;\\lambda_k(\\overrightarrow{w_p}_R) \\rightarrow u_*+\\sqrt{g\\ h_*} &lt; u_R+\\sqrt{g\\ h_R} \\rightarrow h_* &lt; h_R\\] and \\[ u_* &lt; u_R\\]\n\n\n1.2.5.3 Study of the shocks\nConsider a state \\(\\overrightarrow{w_p}_g\\) and a state \\(\\overrightarrow{w_p}_d\\), connected by a \\(k\\) shock wave. The following conditions must then be satisfied:\n\nThe Rankine-Hugoniot conditions, where \\(\\dot\\sigma_k\\) denotes the shock speed: \\(\\displaystyle \\dot\\sigma_k=\\frac{\\left[\\overrightarrow{F}(\\overrightarrow{w})\\right]}{\\left[\\overrightarrow{w}\\right]}\\)\nThe Lax entropy conditions for a genuinely nonlinear field \\(k\\): \\[ \\left\\{ \\begin{array}{l} \\lambda_{k-1}(\\overrightarrow{w_p}_g) &lt; \\dot\\sigma_k  \\\\ \\lambda_{k}(\\overrightarrow{w_p}_d) &lt; \\dot\\sigma_k &lt; \\lambda_{k}(\\overrightarrow{w_p}_g) \\\\ \\dot\\sigma_k &lt; \\lambda_{k+1}(\\overrightarrow{w_p}_d) \\end{array} \\right.  \\]\n\nApplied to the Saint-Venant problem, the Rankine-Hugoniot conditions give: \\[ \\dot\\sigma_k = \\frac{h_d u_d - h_g u_g}{h_d - h_g} =  \\frac{h_d u_d^2 + \\frac{g}{2}h_d^2 - h_g u_g^2 + \\frac{g}{2}h_g^2}{h_d u_d - h_g u_g} \\] in particular, we can define \\(j_k\\) such that\n\\[ j_k = h_g u_g - \\dot\\sigma_k h_g = h_d u_d - \\dot\\sigma_k h_d .\\]\nBy eliminating the shock speed, we obtain: \\[\\left( u_g - u_d \\right)^2 = \\frac{g}{2 h_d h_g} \\left( h_g - h_d \\right)^2 \\left( h_g + h_d \\right)\\] and therefore \\[ j_k =  \\frac{g}{2} \\frac{h_g - h_d}{u_d - u_g} \\]\nConsider first a 1-shock connecting a state \\(\\overrightarrow{w_p}_L\\) to a state \\(\\overrightarrow{w_p}_*\\), separated by a shock wave. The Lax conditions give: \\[ \\displaystyle\\left\\{ \\begin{array}{l} \\lambda_1(\\overrightarrow{w_p}_*) &lt; \\dot\\sigma_1 &lt; \\lambda_1(\\overrightarrow{w_p}_L) \\\\ \\dot\\sigma_1 &lt; \\lambda_2(\\overrightarrow{w_p}_*) \\end{array} \\right. \\rightarrow\n\\left\\{ \\begin{array}{l} u_*-\\sqrt{g h_*} &lt; \\dot\\sigma_1 &lt; u_L-\\sqrt{g h_L} \\\\ \\dot\\sigma_1 &lt; u_*+\\sqrt{g h_*} \\end{array} \\right.\\]\nThe Rankine-Hugoniot conditions give: \\[ \\displaystyle\\left\\{ \\begin{array}{l} j_1 = h_*\\left( u_*-\\dot\\sigma_1\\right) = h_L\\left( u_L-\\dot\\sigma_1\\right) \\\\  \\displaystyle\\dot\\sigma_1 =  \\frac{h_* u_*^2 + \\frac{g}{2}h_*^2 - h_L u_L^2 - \\frac{g}{2}h_L^2}{h_* u_* - h_L u_L} \\end{array} \\right.\\]\nSo, from the Lax conditions and using the first Rankine-Hugoniot condition, we obtain: \\[ \\displaystyle\\left\\{ \\begin{array}{l} u_*-\\sqrt{g h_*} &lt; \\dot\\sigma_1 \\\\ u_L-\\sqrt{g h_L} &gt;  \\dot\\sigma_1 \\\\  u_*+\\sqrt{g h_*} &gt; \\dot\\sigma_1 \\end{array} \\right. \\rightarrow \\left\\{ \\begin{array}{l} h_*\\sqrt{g h_*} &gt; j_1 \\\\ h_L\\sqrt{g h_L} &lt;  j_1 \\\\  - h_*\\sqrt{g h_*} &lt; j_1 \\end{array} \\right.\\] That is \\[ h_* &gt; h_L \\] and since \\(j_1=\\frac{g}{2} \\frac{h_L - h_*}{u_* - u_L} &gt; 0\\), we have \\[ u_* &lt; u_L \\] and using the expression for \\(\\dot\\sigma_1\\), we find the velocity \\[ u_*=u_L+\\left(h_L-h_*\\right)\\sqrt{\\frac{g}{2}\\frac{h_L + h_*}{h_L h_*}} \\tag{1.12}\\]\nConsider finally a 2-shock connecting a state \\(\\overrightarrow{w_p}_*\\) to a state \\(\\overrightarrow{w_p}_R\\), separated by a shock wave. The Lax conditions give: \\[ \\displaystyle\\left\\{ \\begin{array}{l} \\lambda_1(\\overrightarrow{w_p}_*) &lt; \\dot\\sigma_2 \\\\ \\lambda_2(\\overrightarrow{w_p}_R) &lt; \\dot\\sigma_2 &lt; \\lambda_2(\\overrightarrow{w_p}_*) \\end{array} \\right. \\rightarrow\n\\left\\{ \\begin{array}{l} u_*-\\sqrt{g h_*} &lt; \\dot\\sigma_2  \\\\ &lt; u_R+\\sqrt{g h_R} &lt; \\dot\\sigma_2 &lt; u_*+\\sqrt{g h_*} \\end{array} \\right.\\]\nThe Rankine-Hugoniot conditions give: \\[ j_2 = h_*\\left( u_*-\\dot\\sigma_2 \\right) = h_R \\left( u_R-\\dot\\sigma_2 \\right)= \\frac{g}{2} \\frac{h_* - h_R}{u_R - u_*}\\]\nSo, from the Lax conditions and using the first Rankine-Hugoniot condition, we obtain: \\[ \\displaystyle\\left\\{ \\begin{array}{l} u_*-\\sqrt{g h_*} &lt; \\dot\\sigma_2 \\\\ u_R+\\sqrt{g h_R} &lt;  \\dot\\sigma_2 \\\\  u_*+\\sqrt{g h_*} &gt; \\dot\\sigma_2 \\end{array} \\right. \\rightarrow \\left\\{ \\begin{array}{l} h_*\\sqrt{g h_*} &gt; j_2 \\\\ -h_R\\sqrt{g h_R} &gt;  j_2 \\\\  - h_*\\sqrt{g h_*} &gt; j_2 \\end{array} \\right.\\]\nThat is \\[ h_* &gt; h_R \\] and since \\(j_2 &lt; 0\\), we have \\[ u_* &gt; u_L \\] and using the expression for \\(\\dot\\sigma_2\\), we find the velocity \\[ u_* = u_R + \\left( h_* - h_R \\right) \\sqrt{\\frac{g}{2} \\frac{h_R + h_*}{h_R h_*}}  \\tag{1.13}\\]\n\n\n1.2.5.4 Constructing the solution\nTo determine the solution of the Riemann problem, the state \\(\\overrightarrow{w_p}_*\\) must be determined as a function of the states \\(\\overrightarrow{w_p}_L\\) and \\(\\overrightarrow{w_p}_R\\). The intermediate state is connected to the left state by a 1-wave, using Equation 1.12 and Equation 1.10: \\[ u_* = v^L(h_*) = \\left\\{ \\begin{array}{l} u_L - \\displaystyle\\left( \\sqrt{h_*} - \\sqrt{h_L}\\right) \\frac{2\\sqrt{g}}{\\sqrt{h_L} + \\sqrt{h_*}} \\mbox{ if } h_*  \\leq h_L \\mbox{ (1-rarefaction)}\\\\ u_L - \\displaystyle\\left( \\sqrt{h_*} - \\sqrt{h_L}\\right) \\sqrt{\\frac{g}{2} \\frac{h_* + h_L}{h_* h_L}} \\mbox{ if } h_*  &gt; h_L  \\mbox{ (1-shock)}\\end{array} \\right.   \\] and the intermediate state is connected to the right state by a 2-wave, using Equation 1.11 and Equation 1.13: \\[ u_* = v^R(h_*) = \\left\\{ \\begin{array}{l} u_R + \\displaystyle\\left( \\sqrt{h_*} - \\sqrt{h_R}\\right) \\frac{2\\sqrt{g}}{\\sqrt{h_R} + \\sqrt{h_*}} \\mbox{ if } h_*  \\leq h_R \\mbox{ (2-rarefaction)}\\\\ u_R + \\left( \\sqrt{h_*} - \\sqrt{h_R}\\right) \\displaystyle\\sqrt{\\frac{g}{2} \\frac{h_* + h_R}{h_* h_R}} \\mbox{ if } h_*  &gt; h_R  \\mbox{ (2-shock)} \\end{array} \\right.   \\]\nSo to find \\(h_*\\), the equation \\(v^L(h_*)-v^R(h_*)=0\\) must be solved using a Newton-Raphson scheme.\nBy introducing the function \\(\\zeta(h,\\psi)=\\left\\{ \\begin{array}{l}  \\displaystyle\\frac{2\\sqrt{g}}{\\sqrt{h} + \\sqrt{\\psi}} \\mbox{ if } h  \\leq \\psi \\\\ \\displaystyle\\sqrt{\\frac{g}{2} \\frac{h + \\psi}{h \\psi}} \\mbox{ if } h  &gt; \\psi \\end{array} \\right.\\), we will simply write \\[u_* = u_R + \\left( \\sqrt{h_*} - \\sqrt{h_R}\\right) \\zeta(h_*,h_R)\\]\nThe previous calculations allow us to define the speeds of the shock waves \\(\\dot\\sigma_1 = u_L - h_* \\displaystyle\\sqrt{\\frac{g}{2} \\frac{h_L + h_*}{h_L h_*}}\\) and \\(\\dot\\sigma_2 = u_R + h_*  \\displaystyle\\sqrt{\\frac{g}{2} \\frac{h_R + h_*}{h_R h_*}}\\) .\nand the speeds of the rarefaction waves \\(\\lambda_1 = u_L - \\sqrt{g h_L}\\) , \\(\\lambda_1^* = u_* - \\sqrt{g h_*}\\), \\(\\lambda_2^* = u_* + \\sqrt{g h_*}\\) and \\(\\lambda_2 = u_R + \\sqrt{g h_R}\\).\nThe solution in a rarefaction zone is given by: \\[\\overrightarrow{w_p}_{1det}(x,t) = \\left\\{ \\begin{array}{l}  \\displaystyle \\frac{1}{9g} \\left( u_L+2 \\sqrt{g h_L} - \\frac{x}{t} \\right)^2 \\\\   \\displaystyle \\frac{1}{3} \\left( u_L + 2 \\sqrt{g h_L} + 2 \\frac{x}{t} \\right) \\end{array} \\right\\}\\] and \\[\\overrightarrow{w_p}_{2det}(x,t) = \\left\\{ \\begin{array}{l}  \\displaystyle \\frac{1}{9g} \\left( -u_R + 2 \\sqrt{g h_R} + \\frac{x}{t} \\right)^2 \\\\   \\displaystyle \\frac{1}{3} \\left( u_R - 2 \\sqrt{g h_R} + 2 \\frac{x}{t} \\right)  \\end{array} \\right\\} .\\]\nWe can now write the complete solution of the Riemann problem.\nCase 1: \\(h_* &gt; h_L\\) and \\(h_* &gt; h_R\\), the solution consists of a 1-shock and a 2-shock:\n\\[\\overrightarrow{w_p}(x,t) = \\left\\{ \\begin{array}{ll} \\overrightarrow{w_p}_L & \\mbox{ if } x &lt; \\dot\\sigma_1 t \\\\ \\overrightarrow{w_p}_* & \\mbox{ if } \\dot\\sigma_1 t &lt; x &lt; \\dot\\sigma_2 t \\\\ \\overrightarrow{w_p}_R  & \\mbox{ if } x &gt; \\dot\\sigma_2 t \\end{array} \\right.\\]\nCase 2: \\(h_L &lt; h_* &lt; h_R\\), the solution consists of a 1-shock and a 2-rarefaction: \\[\\overrightarrow{w_p}(x,t) = \\left\\{ \\begin{array}{ll} \\overrightarrow{w_p}_L & \\mbox{ if } x &lt; \\dot\\sigma_1 t \\\\ \\overrightarrow{w_p}_* & \\mbox{ if } \\dot\\sigma_1 t &lt; x &lt; \\lambda_2^* t \\\\ \\overrightarrow{w_p}_{2det}(x,t)  & \\mbox{ if } \\lambda_2^* t &lt; x &lt; \\lambda_2^R t \\\\  \\overrightarrow{w_p}_R  & \\mbox{ if } x &gt; \\lambda_2^R t \\end{array} \\right.\\]\nCase 3: \\(h_R &lt; h_* &lt; h_L\\), the solution consists of a 1-rarefaction and a 2-shock: \\[\\overrightarrow{w_p}(x,t) = \\left\\{ \\begin{array}{ll} \\overrightarrow{w_p}_L & \\mbox{ if } x &lt; \\lambda_1^L t \\\\ \\overrightarrow{w_p}_{1det}(x,t) & \\mbox{ if } \\lambda_1^L t &lt; x &lt; \\lambda_1^* t \\\\ \\overrightarrow{w_p}_*  & \\mbox{ if } \\lambda_1^* t &lt; x &lt; \\dot\\sigma_2 t \\\\  \\overrightarrow{w_p}_R  & \\mbox{ if } x &gt; \\dot\\sigma_2 t \\end{array} \\right.\\]\nCase 4: \\(h_* &lt; h_L\\) and \\(h_* &lt; h_R\\), the solution consists of a 1-rarefaction and a 2-rarefaction: \\[\\overrightarrow{w_p}(x,t) = \\left\\{ \\begin{array}{ll} \\overrightarrow{w_p}_L & \\mbox{ if } x &lt; \\lambda_1^L t \\\\ \\overrightarrow{w_p}_{1det}(x,t) & \\mbox{ if } \\lambda_1^L t &lt; x &lt; \\lambda_1^* t \\\\ \\overrightarrow{w_p}_* & \\mbox{ if } \\lambda_1^* t &lt; x &lt; \\lambda_2^* t \\\\ \\overrightarrow{w_p}_{2det}(x,t) & \\mbox{ if } \\lambda_2^* t &lt; x &lt; \\lambda_2^R t\\\\ \\overrightarrow{w_p}_R  & \\mbox{ if } x &gt; \\lambda_2^R t \\end{array} \\right.\\]\n\n\n\n\n\n\nNote\n\n\n\nIn the case where dry zones appear, the algorithm is adapted as explained, for example, in this course by Gloria Faccanoni or in the thesis of Pons (2018).\n\n\n\n\n\n\n\n\nNote\n\n\n\nThis Riemann solver is used in CERF. It is implemented in Fortran90 in ~/sources/PHY/svriemann.F90.\n\n\n\n\n\n\n\n\nImportant 1.2\n\n\n\nIn practice, it is necessary to define a threshold value defining the “dry” zone. In CERF, we therefore define \\(h_{crit}\\) such that if \\(h&lt;h_{crit}\\), the zone is considered dry and \\(h =0\\) and \\(\\overrightarrow{u} = \\overrightarrow{0}\\) are imposed. This value can be modified by the user.",
    "crumbs": [
      "<span class='chapter-number'>1</span>  <span class='chapter-title'>Modeling</span>"
    ]
  },
  {
    "objectID": "chapitre1.html#sec-bf",
    "href": "chapitre1.html#sec-bf",
    "title": "1  Modeling",
    "section": "1.3 Two-fluid flow model",
    "text": "1.3 Two-fluid flow model\n\n1.3.1 Governing equations\nThis model was originally developed to simulate flows in marine hydrodynamics, wave propagation and wave breaking. For breaking, the most physically relevant approach is to consider the two fluids, air and water, with behavior governed by an incompressible Navier-Stokes model, with or without turbulence. Handling the interface is, moreover, a topic in its own right! Whatever discretization approach is used, these simulations are very costly in terms of computation time. Moreover, due to numerical dissipation, it is difficult to simulate flows at very large scales. This is why, to simulate wave propagation phenomena, models of the Saint-Venant, Serre-Green-Naghdi, or inviscid irrotational flow type are often preferred. These models, however, model breaking phenomena poorly, or not at all. The literature on these topics is very extensive, and we do not claim to provide a state of the art here. We will simply offer a partial comparative study in Helluy, Philippe et al. (2005).\nThe two-fluid model presented here stems from the desire to find a compromise between the relevance of Navier-Stokes-type models and the simplicity of Saint-Venant-type models. We therefore consider a mixture of two fluids, so that no interface “tracking” is required and no surface tension can be taken into account. We further assume, which holds true in most cases studied in hydrodynamics, that viscosity can be neglected, which, from a numerical standpoint, advantageously frees us from modeling the Laplacian operator. Finally, as in Chorin (1967), we consider that the two fluids are very weakly compressible. This frees us from the need to impose incompressibility, which is costly in computation time; however, it becomes necessary to add an equation of state for the mixture. After investigation, according to Golay and Helluy (2007), it appears that an isothermal equation of state is sufficient for most of the cases studied, so that we can also dispense with solving the energy equation. We thus obtain a Eulerian model of the form:\n\\[ \\frac{\\partial \\rho}{\\partial t} + \\nabla \\cdot (\\rho \\overrightarrow{u}) = 0  \\tag{1.14}\\] \\[ \\frac{\\partial \\rho\\overrightarrow{u}}{\\partial t} + \\nabla \\cdot \\left(\\rho \\overrightarrow{u} \\otimes \\overrightarrow{u} + p \\overline{\\overline{I}} \\right) = \\rho\\overrightarrow{g}  \\tag{1.15}\\]\nwhere \\(\\rho\\) is the density of the mixture, \\(\\overrightarrow{u}=&lt;u,v,w&gt;^T\\) is the velocity of the mixture, \\(\\overrightarrow{g}\\) is the gravity vector, and where \\(\\overline{\\overline{I}}\\) denotes the identity matrix, \\(\\nabla \\cdot\\) the divergence operator, and \\(\\otimes\\) the tensor product. Equation 1.14 represents conservation of mass and Equation 1.15 represents conservation of momentum. \\(p\\) is the pressure of the mixture, determined by the equation of state:\n\\[ p = p_0 + c_0^2 \\left(\\rho - \\left( \\phi \\rho_a + (1 - \\phi) \\rho_w  \\right)  \\right)  \\tag{1.16}\\]\nwhere \\(\\rho_a\\) represents the density of air, \\(\\rho_w\\) the density of water, \\(p_0\\) the reference pressure, and \\(c_0\\) the speed of sound in the mixture. \\(\\phi\\) represents the volume fraction of air in the mixture. \\(\\phi\\) satisfies a transport equation, which we take in non-conservative form (Golay and Helluy (2007)), as follows:\n\\[ \\frac{\\partial \\phi}{\\partial t} + \\overrightarrow{u} \\cdot \\nabla \\phi  = 0  \\tag{1.17}\\]\nTypically, we take \\(p_0 = 10^5 Pa\\), \\(\\rho_a = 1. kg/m^3\\), \\(\\rho_w = 1000. kg/m^3\\) and \\(c_0 = 20. m/s\\).\nIn hydrodynamics, flows are on the order of a few \\(m/s\\), so that a speed of sound on the order of \\(20 m/s\\) leads to a Mach number \\(M = \\displaystyle \\frac{||\\overrightarrow{u}||}{c_0} &lt; 0.3\\). The mixture can then be considered nearly incompressible. This speed of sound is entirely artificial and does not correspond to any physical reality. It is chosen for numerical reasons. Indeed, the choice of \\(c_0\\) influences the stability of the explicit numerical scheme used to solve the equations, through the CFL condition (Important 1.1).\nHowever, as it happens, according to Wood’s equation (Figure 1.9), while the speed of sound in water is on the order of 1481 \\(m/s\\) and in air on the order of 343 \\(m/s\\), the speed of sound in an air-water mixture is on the order of 23 \\(m/s\\). This is in no way a justification for the model used here, but it is interesting to note that the speed of sound in an air-water mixture is of the same order as the artificial speed of sound used in the two-fluid model.\n\n\n\n\n\n\n\nFigure 1.9: Speed of sound in an air-water mixture according to Wood’s equation\n\n\n\n\nEquations Equation 1.14, Equation 1.15 and Equation 1.17, together with Equation 1.16, constitute the hyperbolic problem solved by CERF to model two-fluid flows. We define the conservative variables \\(\\overrightarrow{w}=&lt;\\rho,\\rho u,\\rho v,\\rho w,\\phi&gt;^T\\) and the primitive variables \\(\\overrightarrow{w_p}=&lt;\\rho, u, v, w, p&gt;^T\\). This corresponds to the general form Equation 1.1.\n\n\n\n1.3.2 Refinement criterion\nFor the two-fluid model, to define the numerical entropy production Equation 1.4, we use the entropy defined by: \\[ s(\\overrightarrow{w}) = \\frac{1}{2}\\rho ||\\overrightarrow{u}||^2 + c_0^2 \\rho ln(\\rho) - c_0^2 \\left( \\rho_w - \\rho_a \\right) \\phi\\] and the entropy flux \\[ \\overrightarrow{\\psi}_s(\\overrightarrow{w}) = \\left( s + c_0^2 \\rho + c_0^2 \\left( \\rho_w - \\rho_a \\right) \\right) \\overrightarrow{u} = \\left( \\frac{1}{2}\\rho ||\\overrightarrow{u}||^2 + c_0^2 \\rho \\left( ln(\\rho) + 1 \\right) \\right) \\overrightarrow{u} .\\]\n\n\n\n1.3.3 Free-surface sharpening\nThe finite-volume resolution scheme for the hyperbolic system Equation 1.14, Equation 1.15 and Equation 1.17, even at second order, is known to be dissipative. As a result, the air-water interface diffuses and widens. To counteract this phenomenon, an interface sharpening model can be used. The idea is to add a source term to the transport equation for \\(\\phi\\) Equation 1.17, which penalizes mixing zones using a splitting method. Over a fictitious time step, we therefore compute: \\[ \\frac{\\partial \\phi}{\\partial \\tau} = \\phi^2 (1 - \\phi)^2(\\phi - c)\\] where \\(c\\) is determined so as to conserve the volume fraction, i.e.:\n\\[ c = \\frac{\\int _{\\Omega} \\phi^3 (1 - \\phi)^2 d\\Omega}{\\int _{\\Omega} \\phi^2 (1 - \\phi)^2 d\\Omega}\\]\nNote that for \\(\\phi=0\\) or \\(\\phi=1\\), this source term has no effect. Finally, since \\(\\phi\\) is modified, for the sake of conservation we also modify the density and momentum of the mixture:\n\\[ \\frac{\\partial \\rho}{\\partial \\tau}  = \\phi^2 (1 - \\phi)^2(\\phi - c)(\\rho_a - \\rho_w)\\] \\[ \\frac{\\partial \\left(\\rho \\overrightarrow{u} \\right)}{\\partial \\tau} = \\phi^2 (1 - \\phi)^2(\\phi - c) \\overrightarrow{u}(\\rho_a - \\rho_w)\\]\nThe fictitious time step is calibrated to represent a fraction of the flux. This fraction is determined by a sharpening coefficient \\(\\alpha_{sharp}\\). In practice, this coefficient is on the order of \\(10^{-3} - 10^{-2}\\).\n\n\n\n1.3.4 The “rigid fluid” concept\nFollowing the fictitious domain principle, part of the cells can be penalized to simulate an obstacle by imposing a zero velocity. Following the work of Coquerelle and Cottet (2008), part of the cells can also be penalized to obtain “rigid solid”-type behavior. The extension of the work of Coquerelle and Cottet (2008) to a two-fluid model was carried out by Altazin (2017). Consider a “solid” occupying the domain \\(\\Omega_s\\), of uniform density \\(\\rho_s\\), center of gravity \\(G\\), mass \\(M_s\\), and inertia matrix \\(\\overline{\\overline{J_s}}\\). The velocity of any point \\(P\\) of the solid is given by:\n\\[ \\overrightarrow{v}_s = \\overrightarrow{v}_G + \\overrightarrow{\\omega}_G \\wedge \\overrightarrow{GP}\\]\nwhere \\(\\overrightarrow{v}_G\\) is the velocity of the center of gravity, \\(\\overrightarrow{\\omega}_G\\) is the angular velocity of the solid, and \\(\\overrightarrow{GP}\\) is the vector connecting the center of gravity to the point \\(P\\). By minimizing \\(\\int _{\\Omega_s} || \\rho\\overrightarrow{u} - \\rho_s\\overrightarrow{v}_s||^2 d\\Omega_s\\), we obtain, for a uniform solid density: \\[\\overrightarrow{v}_G=\\frac{\\int _{\\Omega_s} \\rho\\overrightarrow{u} d\\Omega_s}{M_s}\\]\nand\n\\[\\overrightarrow{\\omega}_G = \\overline{\\overline{J_s}}^{-1} \\int _{\\Omega_s} \\overrightarrow{GP} \\wedge \\overrightarrow{u} .\\]\nAt the end of each time step, the position and velocity of the solid’s center of gravity and the solid’s rotational velocity are computed. \\(\\rho_s\\) and \\(\\overrightarrow{v}_s\\) are then assigned to all the cells of the “solid”. In CERF, the solid is determined by its envelope. This envelope is defined by a triangular mesh. The cells inside the envelope are determined using a “ray casting” technique. In CERF, a solid can be considered “fixed”, “with imposed velocity”, or “free”.\n\n\n\n1.3.5 Riemann solver for the two-fluid model\n\n1.3.5.1 Overview\nNumerical fluxes in CERF are computed using a Riemann solver. There are many references concerning this solver, such as Toro (1997). This section describes the Riemann solver used for the two-fluid flow model, based on Golay and Helluy (2007).\nWe therefore seek to solve the isothermal two-fluid problem in one dimension. We consider a flow with density \\(\\rho(x,t)\\), horizontal velocity \\(u(x,t)\\), pressure \\(p(x,t)\\), and volume fraction \\(\\phi(x,t)\\), described by the following problem: \\[\n\\left \\{\n\\begin{array}{ll}\n\\displaystyle\\frac{\\partial \\rho}{\\partial t} + \\frac{\\partial (\\rho u)}{\\partial x} = 0 \\\\\n\\displaystyle\\frac{\\partial (\\rho u)}{\\partial t} + \\frac{\\partial (\\rho u^2 + p)}{\\partial x} = 0 \\\\\n\\displaystyle\\frac{\\partial \\rho \\phi}{\\partial t} + \\frac{\\partial \\rho \\phi u}{\\partial x} = 0 \\\\\n\\end{array}\n\\right.\n\\tag{1.18}\\] With an equation of state of the form: \\[ p = p_0 + c_0^2 \\left(\\rho - \\left( \\phi \\rho_a + (1 - \\phi) \\rho_w  \\right)  \\right) \\] and initial conditions such that: \\[\n\\overrightarrow{w_p}(x,0) =\\left\\{\n\\begin{array}{ll}\n\\overrightarrow{w_p}_L &\\mbox{ if } x&lt;0 \\\\\n\\overrightarrow{w_p}_R &\\mbox{ if } x&gt;0\n\\end{array}\n\\right.\n\\] where \\(\\overrightarrow{w_p}_L=&lt;\\rho_L, u_L, \\phi_L&gt;^T\\) and \\(\\overrightarrow{w_p}_R=&lt;\\rho_R, u_R, \\phi_R&gt;^T\\) are two constant states.\nBy expanding Equation 1.18, the system can be rewritten in the following form: \\[\n\\frac{\\partial}{\\partial t} \\left\\{ \\begin{array}{c} \\rho \\\\ u \\\\ \\phi \\end{array} \\right\\} + \\left[ \\begin{array}{ccc} u & \\rho & 0 \\\\ \\displaystyle\\frac{c_0^2}{\\rho} & u & \\displaystyle\\frac{c_0^2\\left( \\rho_w - \\rho_a \\right)}{\\rho} \\\\ 0 & 0 & u \\end{array} \\right] \\frac{\\partial}{\\partial x} \\left\\{ \\begin{array}{c} \\rho \\\\ u \\\\ \\phi \\end{array} \\right\\} = \\left\\{ \\begin{array}{c} 0 \\\\ 0 \\\\ 0\\end{array} \\right\\} .\n\\]\nWe then seek the eigenvalues \\(\\lambda_1\\), \\(\\lambda_2\\), and \\(\\lambda_3\\) as well as the eigenvectors \\(\\overrightarrow{r}_1\\), \\(\\overrightarrow{r}_2\\), and \\(\\overrightarrow{r}_3\\) of the matrix thus obtained. A quick calculation gives: \\[\n\\lambda_1=u - c_0 , \\lambda_2=u, \\lambda_3=u + c_0\n\\] with \\(\\lambda_1 \\leq \\lambda_2 \\leq \\lambda_3\\) and \\[\\overrightarrow{r}_1=\\left\\{ \\begin{array}{c} -\\rho/c_0 \\\\ 1 \\\\0 \\end{array} \\right\\} , \\overrightarrow{r}_2=\\left\\{ \\begin{array}{c} \\left( \\rho_w - \\rho_a \\right) \\\\ 0 \\\\ -1 \\end{array} \\right\\}, \\overrightarrow{r}_3=\\left\\{ \\begin{array}{c} \\rho/c_0 \\\\ 1 \\\\0 \\end{array} \\right\\} .\\]\n\n\n\n\n\n\nNote 1.2\n\n\n\nThe eigenvalues are used to determine the wave propagation speed within the cell Important 1.1 as: \\(\\lambda_k = Max\\left( |u - c_0| , |u|, |u + c_0| \\right)\\)\n\n\nTo determine the nature of the three characteristic fields, we compute \\(\\nabla_{w_p}\\lambda_1 \\cdot \\overrightarrow{r}_1 = 1 \\neq 0\\), \\(\\nabla_{w_p}\\lambda_2 \\cdot \\overrightarrow{r}_2 = 0\\) and \\(\\nabla_{w_p}\\lambda_3 \\cdot \\overrightarrow{r}_3 = 1 \\neq 0\\). Characteristic fields 1 and 3 are therefore genuinely nonlinear, while field 2 is linearly degenerate. There will therefore be one contact discontinuity and two waves, which are either shocks or rarefaction waves.\nSince there are \\(3\\) waves, \\(2\\) Riemann invariants can be defined for each wave. The Riemann invariants of wave \\(k\\) are such that \\(\\nabla_{w_p}I_{k;1} \\cdot \\overrightarrow{r}_k = 0\\) and \\(\\nabla_{w_p}I_{k;2} \\cdot \\overrightarrow{r}_k = 0\\). It is easy to verify that:\n\\(I_{1;1} = \\phi\\) , \\(I_{1;2} = u + c_0 ln(\\rho)\\),\n\\(I_{2;1} = u\\) , \\(I_{2;2} = p\\),\n\\(I_{3;1} = \\phi\\) , \\(I_{3;2} = u - c_0 ln(\\rho)\\).\nThe entropy weak solution of the Riemann problem consists of at most four constant states \\(\\overrightarrow{w_p}_L\\), \\(\\overrightarrow{w_p}_I\\), \\(\\overrightarrow{w_p}_{II}\\), \\(\\overrightarrow{w_p}_R\\), separated by a shock or a rarefaction wave, a contact discontinuity, then a shock or a rarefaction wave (Figure 1.10). In what follows, we determine the states \\(\\overrightarrow{w_p}_I\\) and \\(\\overrightarrow{w_p}_{II}.\\)\n\n\n\n\n\n\n\nFigure 1.10: Diagram of the solutions to the Riemann problem with the isothermal two-fluid model\n\n\n\n\n\n\n1.3.5.2 Study of the rarefaction waves\nConsider a state \\(\\overrightarrow{w_p}_g\\) and a state \\(\\overrightarrow{w_p}_d\\), connected by a k-rarefaction wave. The following conditions must then be satisfied:\n\n\\(\\lambda_k(\\overrightarrow{w_p}_g)&lt;\\lambda_k(\\overrightarrow{w_p}_d)\\)\nThe Riemann invariants are conserved along the rarefaction wave \\(I_{k;i}(\\overrightarrow{w_p}_g)=I_{k;i}(\\overrightarrow{w_p}_d)\\) for \\(i=1,2\\)\n\nLet us apply these conditions to the isothermal two-fluid problem.\nConsider first a 1-rarefaction connecting a state \\(\\overrightarrow{w_p}_L\\) to a state \\(\\overrightarrow{w_p}_I\\), separated by a rarefaction wave. The Riemann invariants give:\n\\[\\left\\{ \\begin{array}{l}I_{1;1}(\\overrightarrow{w_p}_L) = I_{1;1}(\\overrightarrow{w_p}_I) \\\\ I_{1;2}(\\overrightarrow{w_p}_L) = I_{1;2}(\\overrightarrow{w_p}_I)\\end{array} \\right. \\rightarrow \\left\\{ \\begin{array}{l}\\phi_L=\\phi_I \\\\ u_L + c_0 ln(\\rho_L) = u_I + c_0 ln(\\rho_I)\\end{array} \\right.\\] that is \\[ u_I = u_L + c_0 ln\\left( \\frac{\\rho_L}{\\rho_I}\\right)  \\tag{1.19}\\] Furthermore, \\[\\lambda_1(\\overrightarrow{w_p}_L)&lt;\\lambda_1(\\overrightarrow{w_p}_I) \\rightarrow u_L &lt; u_I \\] and consequently, from Equation 1.19, \\(\\rho_L &gt; \\rho_I\\).\nConsider finally a 3-rarefaction connecting a state \\(\\overrightarrow{w_p}_{II}\\) to a state \\(\\overrightarrow{w_p}_R\\), separated by a rarefaction wave. The condition on the velocities is written: \\[\\lambda_3(\\overrightarrow{w_p}_{II})&lt;\\lambda_3(\\overrightarrow{w_p}_R) \\rightarrow u_{II} &lt; u_R .\\]\nThe Riemann invariants give \\[\\left\\{ \\begin{array}{l}I_{3;1}(\\overrightarrow{w_p}_{II}) = I_{3;1}(\\overrightarrow{w_p}_R) \\\\ I_{3;2}(\\overrightarrow{w_p}_{II}) = I_{3;2}(\\overrightarrow{w_p}_R)\\end{array} \\right. \\rightarrow \\left\\{ \\begin{array}{l}\\phi_{II}=\\phi_R \\\\ u_{II} - c_0 ln(\\rho_{II}) = u_R - c_0 ln(\\rho_R)\\end{array} \\right.\\] that is \\[ u_{II} = u_R - c_0 ln\\left( \\frac{\\rho_R}{\\rho_{II}}\\right)  \\tag{1.20}\\]\nso since \\(u_{II} &lt; u_R\\), we have \\(\\rho_{II} &lt; \\rho_R\\).\n\n\n1.3.5.3 Study of the shocks\nConsider a state \\(\\overrightarrow{w_p}_g\\) and a state \\(\\overrightarrow{w_p}_d\\), connected by a \\(k\\) shock wave. The following conditions must then be satisfied:\n\nThe Rankine-Hugoniot conditions, where \\(\\dot\\sigma_k\\) denotes the shock speed: \\(\\displaystyle \\dot\\sigma_k=\\frac{\\left[\\overrightarrow{F}(\\overrightarrow{w})\\right]}{\\left[\\overrightarrow{w}\\right]}\\)\nThe Lax entropy conditions for a genuinely nonlinear field \\(k\\): \\[ \\left\\{ \\begin{array}{l} \\lambda_{k-1}(\\overrightarrow{w_p}_g) &lt; \\dot\\sigma_k  &lt; \\lambda_{k}(\\overrightarrow{w_p}_g) \\\\ \\lambda_{k}(\\overrightarrow{w_p}_d) &lt; \\dot\\sigma_k &lt; \\lambda_{k+1}(\\overrightarrow{w_p}_d) \\end{array} \\right.  \\]\n\nApplied to the two-fluid problem, the Rankine-Hugoniot conditions give: \\[ \\dot\\sigma_k = \\frac{\\rho_d u_d - \\rho_g u_g}{\\rho_d - \\rho_g} =  \\frac{\\rho_d u_d^2 + p_d - \\rho_g u_g^2 - p_g}{\\rho_d u_d - \\rho_g u_g} = \\frac{\\rho_d u_d \\phi_d - \\rho_g u_g \\phi_g}{\\rho_d \\phi_d - \\rho_g \\phi_g} \\] \\[ \\left\\{ \\begin{array}{l} \\rho_d \\dot\\sigma_k - \\rho_g \\dot\\sigma_k =  \\rho_d u_d - \\rho_g u_g \\\\ \\rho_d u_d \\dot\\sigma_k - \\rho_g u_g \\dot\\sigma_k = \\rho_d u_d^2 + p_d - \\rho_g u_g^2 - p_g \\\\ \\rho_d \\phi_d \\dot\\sigma_k - \\rho_g \\phi_g \\dot\\sigma_k =  \\rho_d u_d \\phi_d - \\rho_g u_g \\phi_g\\end{array} \\right.\\] \\[ \\left\\{ \\begin{array}{l}  \\rho_g (u_g - \\dot\\sigma_k) =  \\rho_d (u_d - \\dot\\sigma_k) \\\\ u_d \\rho_d (u_d - \\dot\\sigma_k) -u_g \\rho_g (u_g - \\dot\\sigma_k) + p_d -  p_g = 0 \\\\ \\phi_g \\rho_g (u_g - \\dot\\sigma_k) =  \\phi_d \\rho_d (u_d - \\dot\\sigma_k) \\end{array} \\right.\\] So, denoting \\(j_k = \\rho_g (u_g - \\dot\\sigma_k) =  \\rho_d (u_d - \\dot\\sigma_k)\\), we obtain: \\[\\phi_g = \\phi_d\\]\nand\n\\[j_k = \\frac{p_d - p_g}{u_g - u_d}\\]\nBy eliminating the shock speed, we also have: \\[\\left( \\rho_d u_d - \\rho_g u_g \\right)^2 = \\left( \\rho_d - \\rho_g \\right) \\left( \\rho_d u_d^2 - \\rho_g u_g^2 + p_d - p_g \\right)\\] or equivalently\n\\[\\rho_d \\rho_g \\left( u_d - u_g \\right)^2 = \\left( \\rho_d - \\rho_g \\right) \\left( p_d - p_g \\right)\\]\nConsider first a 1-shock connecting a state \\(\\overrightarrow{w_p}_L\\) to a state \\(\\overrightarrow{w_p}_I\\), separated by a shock wave. The Lax admissibility conditions are written: \\[ \\displaystyle\\left\\{ \\begin{array}{l} \\dot\\sigma_1 &lt; \\lambda_1(\\overrightarrow{w_p}_L) \\\\ \\lambda_1(\\overrightarrow{w_p}_I) &lt;\\dot\\sigma_1 &lt; \\lambda_2(\\overrightarrow{w_p}_I) \\end{array} \\right. \\rightarrow\n\\left\\{ \\begin{array}{l} \\dot\\sigma_1  &lt; u_L - c_0 \\\\ u_I - c_0 &lt; \\dot\\sigma_1 \\\\ \\dot\\sigma_1 &lt; u_I \\end{array} \\right. \\rightarrow \\left\\{ \\begin{array}{l} \\rho_L c_0 &lt; j_1 \\\\ j_1 &lt; \\rho_I c_0 \\\\ \\dot\\sigma_1 &lt; u_I \\end{array} \\right.\\]\nSo in particular \\(u_I &lt; u_L\\) and \\(\\rho_L &lt; \\rho_I\\). The Rankine-Hugoniot conditions give: \\[ \\displaystyle\\left\\{ \\begin{array}{l} \\phi_L =  \\phi_I \\\\  \\left( u_I - u_L  \\right)^2 = \\displaystyle\\frac{\\rho_I - \\rho_L}{\\rho_I \\rho_L}\\left( p_I - p_L  \\right) \\end{array} \\right.\\] from which, since \\(u_I &lt; u_L\\) and \\(\\rho_L &lt; \\rho_I\\): \\[ u_I=u_L -\\sqrt{\\frac{\\rho_I - \\rho_L}{\\rho_I \\rho_L}\\left( p_I - p_L  \\right)} = u_L + c_0 \\frac{\\rho_L - \\rho_I}{\\sqrt{\\rho_I \\rho_L}} \\tag{1.21}\\]\nConsider finally a 3-shock connecting a state \\(\\overrightarrow{w_p}_{II}\\) to a state \\(\\overrightarrow{w_p}_R\\), separated by a shock wave. The Lax conditions give: \\[ \\displaystyle\\left\\{ \\begin{array}{l} \\lambda_2(\\overrightarrow{w_p}_{II}) &lt; \\dot\\sigma_3 &lt; \\lambda_3(\\overrightarrow{w_p}_{II}) \\\\ \\lambda_3(\\overrightarrow{w_p}_R) &lt; \\dot\\sigma_3  \\end{array} \\right. \\rightarrow\n\\left\\{ \\begin{array}{l} u_{II} &lt; \\dot\\sigma_3 \\\\   \\dot\\sigma_3 &lt; u_{II} + c_0 \\\\ u_R + c_0 &lt; \\dot\\sigma_3 \\end{array} \\right. \\rightarrow\n\\left\\{ \\begin{array}{l} u_{II} &lt; \\dot\\sigma_3 \\\\   - c_0 \\rho_{II} &lt; j_3 \\\\ j_3 &lt; - c_0 \\rho_R \\end{array} \\right.\\]\nSo in particular \\(u_R &lt; u_{II}\\) and \\(\\rho_R &lt; \\rho_{II}\\). The Rankine-Hugoniot conditions are written: \\[ \\displaystyle\\left\\{ \\begin{array}{l} \\phi_{II} =  \\phi_R \\\\  \\left( u_R - u_{II}  \\right)^2 = \\displaystyle\\frac{\\rho_R - \\rho_{II}}{\\rho_R \\rho_{II}}\\left( p_R - p_{II}  \\right) \\end{array} \\right.\\] from which, since \\(u_R &lt; u_{II}\\) and \\(\\rho_R &lt; \\rho_{II}\\): \\[ u_{II}=u_R + \\sqrt{\\frac{\\rho_R - \\rho_{II}}{\\rho_R \\rho_{II}}\\left( p_R - p_{II}  \\right)} = u_R - c_0 \\frac{\\rho_R - \\rho_{II}}{\\sqrt{\\rho_R \\rho_{II}}} \\tag{1.22}\\]\n\n\n1.3.5.4 Study of the contact discontinuity\nConsider a state \\(\\overrightarrow{w_p}_g\\) and a state \\(\\overrightarrow{w_p}_d\\), connected by a contact discontinuity. Across a k-contact discontinuity, entropy and the Riemann invariants are conserved. Since the discontinuity moves at speed \\(\\dot\\sigma_k\\), we have: \\[\\displaystyle \\dot\\sigma_k = \\lambda_{k}(\\overrightarrow{w_p}_g) = \\lambda_{k}(\\overrightarrow{w_p}_d) = \\frac{\\left[\\overrightarrow{F}(\\overrightarrow{w})\\right]}{\\left[\\overrightarrow{w}\\right]}\\]\nConsider then a 2-contact discontinuity connecting the state \\(\\overrightarrow{w_p}_I\\) to the state \\(\\overrightarrow{w_p}_{II}\\). We then have \\[ \\dot\\sigma_2 = u_I = u_{II} = \\frac{\\rho_{II} u_{II} - \\rho_I u_I}{\\rho_{II} - \\rho_I} =  \\frac{\\rho_{II} u_{II}^2 + p_{II} - \\rho_I u_I^2 - p_I}{\\rho_{II} u_{II} - \\rho_I u_I} = \\frac{\\rho_{II} u_{II} \\phi_{II} - \\rho_I u_I \\phi_I}{\\rho_{II} \\phi_{II} - \\rho_I \\phi_I} \\]\nFrom which we deduce: \\[ u_I = u_{II} \\mbox{  and  } p_I = p_{II} \\tag{1.23}\\]\n\n\n1.3.5.5 Constructing the solution\nTo simplify notation, let us define a function \\(H(a,b)\\) such that \\(H(a,b) = \\displaystyle\\left\\{ \\begin{array}{ll} c_0 \\displaystyle\\frac{\\rho_a - \\rho_b}{\\sqrt{\\rho_a \\rho_b}} & \\mbox{  if  } \\rho_a &lt; \\rho_b \\\\  c_0 ln \\left(\\displaystyle\\frac{\\rho_a}{\\rho_b} \\right) & \\mbox{  if  } \\rho_a &gt; \\rho_b  \\end{array} \\right.\\).\nFor the first wave, using Equation 1.19, Equation 1.21, we can then write \\(u_I = u_L + H(\\rho_L, \\rho_I)\\), and for the third wave, using Equation 1.20, Equation 1.22, \\(u_{II} = u_R - H(\\rho_R, \\rho_{II})\\).\nWe know from the study of the contact discontinuity that: \\[u_I - u_{II} = 0\\] so\n\\[u_L + H(\\rho_L, \\rho_I) - u_R + H(\\rho_R, \\rho_{II}) = 0\\] Now, \\(\\rho_I=\\rho_I(p_I,\\phi_I)\\) and \\(\\rho_{II}=\\rho_{II}(p_{II},\\phi_{II})\\). We also know that \\(p_I = p_{II}\\); we will therefore denote the intermediate pressure \\(p^*\\) such that \\(p^* = p_I = p_{II}\\). Moreover, we know \\(\\phi_I = \\phi_L\\) and \\(\\phi_{II} = \\phi_R\\). We can then write: \\[u_L - u_R + H(\\rho_L, \\rho(p^*,\\phi_L)) + H(\\rho_R, \\rho(p^*,\\phi_R)) = 0\\] This is therefore a nonlinear equation in \\(p^*\\) that can be solved using a Newton-Raphson method. Once \\(p^*\\) is known, the states \\(\\overrightarrow{w_p}_I\\) and \\(\\overrightarrow{w_p}_{II}\\) can be determined.\n\n\n\n\n\n\nNote\n\n\n\nThis Riemann solver is used in CERF. It is implemented in Fortran90 in ~/sources/PHY/riemisot.F90.\n\n\n\n\n\n\n\nAltazin, Thomas. 2017. “Un Modèle d’interaction Fluide-Structure En Régime Compressible Faible Mach.” PhD thesis. http://www.theses.fr/2017TOUL0021/document.\n\n\nAudusse, Emmanuel, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoı⁁t Perthame. 2004. “A Fast and Stable Well-Balanced Scheme with Hydrostatic Reconstruction for Shallow Water Flows.” SIAM Journal on Scientific Computing 25 (6): 2050–65. https://doi.org/10.1137/S1064827503431090.\n\n\nChorin, Alexandre Joel. 1967. “A Numerical Method for Solving Incompressible Viscous Flow Problems.” Journal of Computational Physics 2 (1): 12–26. https://doi.org/https://doi.org/10.1016/0021-9991(67)90037-X.\n\n\nCoquerelle, M., and G.-H. Cottet. 2008. “A Vortex Level Set Method for the Two-Way Coupling of an Incompressible Fluid with Colliding Rigid Bodies.” Journal of Computational Physics 227 (21): 9121–37. https://doi.org/https://doi.org/10.1016/j.jcp.2008.03.041.\n\n\nCroisille, Jean-Pierre. 1990. “Contribution à l’étude Théorique Et à l’approximation Par Éléments Finis Du Système Hyperbolique de La Dynamique Des Gaz Multidimensionnelle Et Multiespèces.” PhD thesis. http://www.theses.fr/1990PA066095.\n\n\nErsoy, Mehmet, Frédéric Golay, and Lyudmyla Yushchenko. 2013. Open Mathematics 11 (8): 1392–1415. https://doi.org/doi:10.2478/s11533-013-0252-6.\n\n\nGolay, Frédéric. 2009. “Numerical Entropy Production and Error Indicator for Compressible Flows.” Comptes Rendus. Mécanique 337 (4): 233–37. https://doi.org/10.1016/j.crme.2009.04.004.\n\n\nGolay, Frédéric, and Philippe Helluy. 2007. “Numerical Schemes for Low Mach Wave Breaking.” International Journal of Computational Fluid Dynamics 21 (2): 69–86. https://doi.org/10.1080/10618560701343382.\n\n\nHelluy, Philippe, Golay, Frédéric, Caltagirone, Jean-Paul, Lubin, Pierre, Vincent, Stéphane, Drevard, Deborah, Marcer, Richard, et al. 2005. “Numerical Simulations of Wave Breaking.” ESAIM: M2AN 39 (3): 591–607. https://doi.org/10.1051/m2an:2005024.\n\n\nLax, Peter. 1971. “Shock Waves and Entropy.” In Contributions to Nonlinear Functional Analysis, edited by Eduardo H. Zarantonello, 603–34. Academic Press. https://doi.org/https://doi.org/10.1016/B978-0-12-775850-3.50018-2.\n\n\nMarche, Fabien. 2007. “Derivation of a New Two-Dimensional Viscous Shallow Water Model with Varying Topography, Bottom Friction and Capillary Effects.” European Journal of Mechanics - B/Fluids 26 (1): 49–63. https://doi.org/https://doi.org/10.1016/j.euromechflu.2006.04.007.\n\n\nPons, Kevin. 2018. “Dynamique Des Écoulements à Surface Libre Dans Les Plages Sableuses.” Thèse de doctorat, IMATH, Université de Toulon: Université de Toulon.\n\n\nPoussel, Camille. 2024. “Dynamique Des Écoulements à Surface Libre Dans Les Plages Sableuses.” Thèse de doctorat, IMATH, Université de Toulon: Université de Toulon.\n\n\nPuppo, Gabriella. 2004. “Numerical Entropy Production for Central Schemes.” SIAM Journal on Scientific Computing 25 (4): 1382–1415. https://doi.org/10.1137/S1064827502386712.\n\n\nSaint-Venant, Barré de. 1871. “AJC (1871b).\" Théorie Et Equations générales Du Mouvement Non Permanent Des Eaux, Avec Application Aux Crues Des Rivières Et à l’introduction Des Marées Dans Leur Lit (2ème Note).\".” Comptes Rendus Des séances de l’Académie Des Sciences, 237–40.\n\n\nThomas Altazin, Frédéric Golay, Mehmet Ersoy, and Lyudmyla Yushchenko. 2016. “Numerical Investigation of BB-AMR Scheme Using Entropy Production as Refinement Criterion.” International Journal of Computational Fluid Dynamics 30 (3): 256–71. https://doi.org/10.1080/10618562.2016.1194977.\n\n\nToro, E. F. 1997. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Springer. https://books.google.fr/books?id=6QFAAQAAIAAJ.\n\n\nVarra, Giada, Veronica Pepe, Renata Della Morte, and Luca Cozzolino. 2024. “A Novel Efficient and Robust Treatment of the Friction Source Term in 2D Shallow Water Inundation Models.” Journal of Hydrology 634: 131045. https://doi.org/10.1016/j.jhydrol.2024.131045.",
    "crumbs": [
      "<span class='chapter-number'>1</span>  <span class='chapter-title'>Modeling</span>"
    ]
  },
  {
    "objectID": "chapitre2.html",
    "href": "chapitre2.html",
    "title": "2  User Guide",
    "section": "",
    "text": "2.1 Software installation\nThis chapter describes the installation of the CERF software, as well as its use.",
    "crumbs": [
      "<span class='chapter-number'>2</span>  <span class='chapter-title'>User Guide</span>"
    ]
  },
  {
    "objectID": "chapitre2.html#software-installation",
    "href": "chapitre2.html#software-installation",
    "title": "2  User Guide",
    "section": "",
    "text": "2.1.1 Prerequisites\nCERF is software written in Fortran 90 using the MPI library for communication between the different processes.\n\n\n\n\n\n\nImportant\n\n\n\nIt is therefore necessary to have:\n- a Linux operating system or emulator (Ubuntu, CentOS, WSL on Windows 10 or 11, Mac OS, ….)\n- a Fortran 90 compiler (gfortran, ifort, …)\n- the MPI library (OpenMPI)\n\n\nThe output files of CERF are in Tecplot format (VTK to come). It is therefore recommended to have visualization software such as ParaView or Visit.\n\n\n2.1.2 Installation\nTo install CERF, after downloading the archive ‘cerf_v1.0.tar’, simply unpack it in a directory of your choice.\ntar xvf cerf_v1.0.tar .\nThe unpacked archive will then present a directory tree as illustrated in Figure 2.1.\n\n\n\n\n\n\n\nFigure 2.1: CERF directory tree\n\n\n\n\nThe main CERF directory contains 4 subdirectories: doc, sources, lib and exec.\n\nthe doc directory contains the software documentation\n\nthe sources directory contains the software sources\n\nthe lib directory contains the libraries created during compilation, required to build the executables\n\nthe exec directory contains all the software’s executables. This is the working directory\n\nTo install CERF, simply go to the exec directory and run the Makefile.\ncd CERF/exec\nmake\n\n\n\n\n\n\nNote\n\n\n\nBy default, the Makefile uses the gfortran compiler and the OpenMPI library. However, these settings can be changed in the Makefile.\n\nParameter FC = “gfortran” or “ifort” or ….. for a sequential computation\n\nParameter FC = “mpif90” or ….. for a parallel computation\n\nParameter MPI = “_MPI” or “NONE_MPI” depending on whether the MPI library is installed or not\n\nParameter FFLAGS to redefine the compilation settings\n\n\n\n\n\n\n\n\n\n\nNote\n\n\n\nThe Makefile is configured to recompile only newly created or modified source files. In the event of a full recompilation, for example following a change of machine, it is advisable to clean the object files and executables before rerunning the Makefile, using the command:\nmake clean\n\n\nOnce the Makefile has finished running, the CERF executables are created in the exec directory: cerf, cerf_amr, cerf_input and cerf2tec.\n\n\n\n\n\n\nNote\n\n\n\nDepending on the operating system used, the executables may or may not have the “.exe” extension.\n\n\n\n\n\n2.1.3 Running a computation\n\n\n\n\n\n\n\nFigure 2.2: Steps of a CERF computation\n\n\n\n\nAs illustrated in Figure 2.2, running a CERF computation involves several steps. A pre-processing phase is used to generate the mesh, the initial conditions of the computation, and the computation parameters. This phase is carried out by the executable cerf_input. Once the mesh has been generated, the computation is performed by the executable cerf, and the mesh can be modified by the executable cerf_amr. Finally, the executable cerf2tec converts the output files to Tecplot format for a post-processing phase.\n\ncerf_input\nAs input, cerf_input reads a configuration file .inp in a format specific to CERF, in which the physical and numerical parameters of the model to be processed, the mesh, the initial conditions, the boundary conditions, etc., are defined. If the mesh is defined by a file in GMSH V2.2 format, it will also read the .msh file. If rigid fluid penalizations are defined, it will also read the objects in the .obj files.\nAs output, cerf_input generates a binary file cerfbin_000 containing the definition of the master mesh, and as many binary files cerfbin_* as required domains, containing all the domain information: parameters, mesh, degrees of freedom, boundary conditions, etc.\ncerf\nAs input, cerf reads the binary files cerfbin_* generated by cerf_input or modified by cerf_amr. It then performs a computation in parallel (via MPI processes) on the defined domains over the required time period. During this step, the mesh is fixed. As output, once the computation is complete, cerf updates the binary files cerfbin_* with the new values of the computed fields.\ncerf_amr\nAs input, cerf_amr reads the binary files cerfbin_* generated by cerf. It then refines the mesh according to criteria defined in the configuration file .inp. As output, cerf_amr updates the binary files cerfbin_* with the new mesh. In addition, cerf_amr generates a file cerfamr.dat containing information related to the master mesh: refinement criterion, refinement level, etc. If rigid fluid penalizations are defined, cerf_amr also generates solide_* files containing the new position of the objects.\ncerf2tec\nAs input, cerf2tec reads the binary files cerfbin_. It then converts the binary files into a single ASCII file, named cerfout.dat*, in Tecplot format for visualizing the results.\n\nRunning a computation using only command lines can become tedious. This is why scripts are often used. Below is an example of a script that lets you run a complete CERF computation with a single command. Simply modify it according to your needs.\n\n\n\n\n\n\nTipExample shell script\n\n\n\nBelow is an example shell script for running a computation from the file dam.inp. This script runs the computation on 4 processors, then modifies the mesh 20 times. The output files are renamed according to the iteration so that the progress of the computation can be visualized.\n#!/bin/bash\n./cerf_input dam.inp\n./cerf2tec\nmv cerfout.dat d_0.dat\nmv cerfamr.dat d_mame_0.dat\nfor i in $(seq 1 20 )\ndo\nmpirun -np 4 ./cerf\n./cerf2tec\nmv cerfout.dat d_$i.dat\n./cerf_amr\nmv cerfamr.dat d_mame_$i.dat\ndone",
    "crumbs": [
      "<span class='chapter-number'>2</span>  <span class='chapter-title'>User Guide</span>"
    ]
  },
  {
    "objectID": "chapitre2.html#the-configuration-file-.inp",
    "href": "chapitre2.html#the-configuration-file-.inp",
    "title": "2  User Guide",
    "section": "2.2 The configuration file *.inp",
    "text": "2.2 The configuration file *.inp\n\n2.2.1 General overview\nThe configuration file is a text file used to define all the parameters and entities required for the simulation. It is divided into several sections, each beginning with a 4-character header. The headers, to be provided in this order, are as follows: PHYS, MAME, BATH, INIT, COND, MESH, OBST, NUME. Each header is followed on the same line by 1 or 2 integers defining the code and, if applicable, the argument. All lines starting with “!” are comments and are ignored by the program. Comments may be placed between sections, but not within a section. Each section is then completed with lines of real or integer data, or keywords, that define the simulation parameters.\n\n\n\n\n\n\nTipExample line with header, code, and argument\n\n\n\n! here is a line with header, code, and argument\nMAME    1    1\nThe first line is a comment line. The second line is a line with a header. It starts with the header MAME, followed by 2 integers, 1 and 1 defining the code and the argument.\n\n\n\n\n2.2.2 Defining the physical parameters: PHYS\nThe PHYS header is used to define the physical parameters of the model. The model is applied to the entire computational domain. The available models are the Saint-Venant model Section 1.2 and the two-fluid model Section 1.3.\nThe model used is selected by setting the keyword MODE, which is mandatory and takes the value \\(2\\) for the two-fluid model and \\(1\\) for the Saint-Venant model.\nFor the isothermal two-fluid model, the following real-valued parameters can be defined as needed:\n\nGPES: Orientation of the gravity field along \\(1/x\\), \\(2/y\\), \\(3/z\\), or none \\(0\\) (default value \\(2.\\))\nPREF: Reference pressure \\(p_0\\) of the equation of state Equation 1.16 (default value \\(10^5\\))\nCSSM: Artificial speed of sound \\(c_0\\) of the equation of state Equation 1.16 (default value \\(20.\\))\nREFA: Air density \\(\\rho_a\\) of the equation of state Equation 1.16 (default value \\(1.\\))\nREFW: Water density \\(\\rho_w\\) of the equation of state Equation 1.16 (default value \\(1000.\\))\nSHAR: Interface sharpening coefficient Section 1.3.3 (default value \\(0.\\))\n\nFor the Saint-Venant model, the following parameters can be defined as needed:\n\nWLTR: Threshold value defining the dry zone (default value \\(10^{-4}\\))\nFRMO: Inclusion of friction. \\(0\\) no friction, \\(1\\) Manning-Strickler, 2/Darcy-Weisbach (default value \\(0\\)). The value of the friction coefficient, which may vary spatially, is defined later in the BATH section.\nBAGR: The bathymetry gradient is either given (by a user-defined function or by bilinear interpolation) BAGR=0, or computed using the Barth method BAGR=1 (default value \\(0\\)).\n\n\n\n\n\n\n\nTipExample parameter definition for the two-fluid model\n\n\n\nPHYS 0\nMODE 2.\nSHAR 0.01\nThis defines a two-fluid model for an air-water flow with gravity along \\(y\\) and an interface sharpening coefficient of \\(0.01\\).\n\n\n\n\n\n\n\n\nTipExample parameter definition for the Saint-Venant model\n\n\n\nPHYS 0\nMODE 1.\nFRMO 1.\nThis defines a shallow-water flow model with Manning-Strickler-type friction.\n\n\n\n\n\n2.2.3 Defining the master mesh: MAME\nThe MAME header is used to define the master mesh. As a reminder, this mesh must be conformal! CERF has a rudimentary mesh generator (code \\(1\\)) and reads ASCII mesh files in GMSH V2.2 format (code \\(2\\)).\nThe mesh generated by CERF is a Cartesian mesh; the domain boundaries must be defined on one line, and the discretization in each direction on another line.\n\nFor a “1D” mesh: argument \\(0\\)\n\n1 line consisting of 2 reals: \\(x_{min}\\), \\(x_{max}\\)\n1 line consisting of 1 integer: \\(n_x\\)\n\nFor a “2D” mesh: argument \\(1\\)\n\n1 line consisting of 4 reals: \\(x_{min}\\), \\(x_{max}\\), \\(y_{min}\\), \\(y_{max}\\)\n1 line consisting of 2 integers: \\(n_x\\), \\(n_y\\)\n\nFor a “3D” mesh: argument \\(2\\)\n\n1 line consisting of 6 reals: \\(x_{min}\\), \\(x_{max}\\), \\(y_{min}\\), \\(y_{max}\\), \\(z_{min}\\), \\(z_{max}\\)\n1 line consisting of 3 integers: \\(n_x\\), \\(n_y\\), \\(n_z\\)\n\n\nIn 1D, the dimensions along \\(y\\) and \\(z\\) are automatically determined by CERF, such that \\(y_{min} = z_{min} = - \\frac{x_{max}-x_{min}}{2 n_x}\\) and \\(y_{max} = z_{max} = \\frac{x_{max}-x_{min}}{2 n_x}\\).\nIn 2D, the dimension along \\(z\\) is automatically determined by CERF, such that \\(z_{min} = - Min \\left(  \\frac{x_{max} - x_{min}}{2 n_x} , \\frac{y_{max} - y_{min}}{2 n_y} \\right)\\) and \\(z_{max} =  Min \\left(  \\frac{x_{max} - x_{min}}{2 n_x} , \\frac{y_{max} - y_{min}}{2 n_y} \\right)\\).\n\n\n\n\n\n\nTipExample of a master mesh definition via CERF\n\n\n\nMAME    1    1\n0. 1.  0.  1.\n2 2\nThis defines a “2D” mesh consisting of 4 hexahedral-shaped blocks, 2 along \\(x\\), 2 along \\(y\\), and 1 along \\(z\\). The blocks have size \\(0.5 \\times 0.5\\) along \\(x\\) and \\(y\\). The dimension along \\(z\\) is automatically determined by CERF.\n\n\nFor more complex meshes, a mesh generated by an external mesh generator, such as GMSH, can be used. In this case, the mesh must be defined in a file in GMSH V2.2 format. The MAME header must then be given code \\(2\\), and the next line contains the name of the file, which must be located in the exec directory.\n\n\n\n\n\n\nTipExample of a master mesh definition via GMSH\n\n\n\nMAME    2    0\ntoto.msh\nThis defines a “2D” mesh from the file toto.msh in GMSH V2.2 format.\n\n\n\n\n\n2.2.4 Defining the initial conditions: INIT\nThe INIT header is used to define the initial conditions of the problem to be solved. The initial conditions can be defined for the shock tube case (code \\(0\\)), via a user-defined function (code \\(1\\)), or by zones (code \\(2\\)).\nAs a reminder, the primitive variables of the Saint-Venant model number \\(n_{dof}=3\\): \\(h\\), \\(u\\), \\(v\\).\nAs a reminder, the primitive variables of the two-fluid model number \\(n_{dof}=5\\): \\(\\rho\\), \\(u\\), \\(v\\), \\(w\\), \\(p\\).\n\nFor a “shock tube” initialization: code \\(0\\)\nThe argument is \\(0\\). Then:\n\n1 line consisting of \\(2\\) reals for the Saint-Venant model (\\(h\\), \\(u\\)) or \\(3\\) reals for the two-fluid model (\\(\\rho\\), \\(u\\), \\(p\\)) at the initial time on the “left” side (\\(x&lt;0\\))\n1 line consisting of \\(2\\) reals for the Saint-Venant model (\\(h\\), \\(u\\)) or \\(3\\) reals for the two-fluid model (\\(\\rho\\), \\(u\\), \\(p\\)) at the initial time on the “right” side (\\(x&gt;0\\))\n\nFor an initialization via a user-defined function: code \\(1\\)\nThe argument defines the number of the user-defined function, written in the routine ~/CERF/sources/PHY/USER/fct_iniw.f90.\nFor an initialization by zones: code \\(2\\)\nThe argument defines the number of zones to create. The order in which zones are defined is important. It allows zones to be “nested”. It is therefore recommended to define the first zone so that it broadly covers the study domain. Then, for each zone, 2 lines are used:\n\n1 line consisting of 6 reals defining the zone: \\(x_{min}\\), \\(x_{max}\\), \\(y_{min}\\), \\(y_{max}\\), \\(z_{min}\\), \\(z_{max}\\)\n1 line consisting of \\(n_{dof}\\) reals defining the primitive variables at the initial time\n\n\n\n\n\n\n\n\nTip 2.1: Example of initialization by zones\n\n\n\nINIT    2    2\n-1. 4.1 -1. 3.1 -1. 1.\n1. 0. 0. 0.  1.e5\n-1. 1. -1. 2. -1. 1.\n1000. 0. 0. 0.  1.e5\nIn this example, we aim to initialize a “dam-break” problem. In a “2D” domain of \\(4 \\times 3\\) filled with air, a column of water of height \\(2\\) and width \\(1\\) is initialized on the left. We therefore define 2 zones. The first zone broadly covers the entire domain, and the second covers the domain occupied by the water. Throughout the domain, we initialize \\(\\rho=1\\), \\(\\overrightarrow{u}=\\overrightarrow{0}\\), and \\(p_0=10.^5\\). With the second zone, in the part occupied by the water, we initialize \\(\\rho=1000\\), \\(\\overrightarrow{u}=\\overrightarrow{0}\\), and \\(p_0=10.^5\\).\n\n\n\n\n\n\n\n\nImportant\n\n\n\nIf a user-defined function is used for initialization, it is necessary to rerun compilation and linking using the make command!\n\n\n\n\n\n2.2.5 Defining bathymetry and friction: BATH\nThe BATH header is used to define the bathymetry and friction parameters, naturally in the case where the Saint-Venant model is selected.\n\nIn the case of a flat bottom with no friction: code \\(0\\)\n\nIn the case of bathymetry and friction defined by a user-defined function: code \\(1\\)\nThe argument defines the number of the user-defined function, written in the routine ~/CERF/sources/PHY/USER/fct_bathy.f90.\nIn the case of bathymetry and friction defined bilinearly per block: code \\(2\\)\nArgument \\(0\\) reads the file, and argument \\(-1\\) creates it. The next line specifies the name of the ASCII file containing the bathymetry and friction information. This file must be located in the exec directory. When creating the file (argument \\(-1\\)), it consists of 2 columns containing the \\(x\\) and \\(y\\) coordinates of the four vertices of each block of the master mesh. The user must then complete the file by adding two columns of reals giving, for each coordinate, the bathymetry and friction values.\n\n\n\n\n2.2.6 Defining the boundary conditions: COND\nThe COND header is used to define the boundary conditions of the problem. Boundary conditions are defined by a number (written on one line) and, if applicable, real values (written on another line). The available boundary conditions are:\n\n0: free outflow condition (copied through)\n1: mirror condition\n2: Dirichlet condition on the primitive variables. The values of the \\(n_{dof}\\) primitive variables must then be defined on the next line. If the value \\(10^{20}\\) is used, the variable is considered free.\n3: user-defined condition. The user must then define the boundary condition function in the routine ~/CERF/sources/PHY/USER/fct_cond.f90. The next line defines the number of the user-defined function.\n4: Dirichlet condition on the conservative variables. The values of the \\(n_{dof}\\) conservative variables must then be defined on the next line. If the value \\(10^{20}\\) is used, the variable is considered free.\n999: No boundary condition. Fluxes are not computed.\n\nThe definition of the boundary conditions at the domain edges then depends on the type of master mesh used:\n\nIn the case of a “shock tube”: code \\(0\\)\nThis is the case where the MAME section has code \\(1\\), argument \\(0\\), and the INIT section has code \\(0\\). In this case, no definition is necessary.\nIn the case of a “2D” mesh: code \\(1\\)\nThis is the case where the MAME section has code \\(1\\), argument \\(1\\). The boundary conditions are then defined on the domain edges in the order \\(x_{min}\\), \\(x_{max}\\), \\(y_{min}\\), \\(y_{max}\\). For each edge, 1 or 2 lines are used to define the boundary condition.\nIn the case of a “3D” mesh: code \\(2\\)\nThis is the case where the MAME section has code \\(1\\), argument \\(2\\). The boundary conditions are then defined on the domain edges in the order \\(x_{min}\\), \\(x_{max}\\), \\(y_{min}\\), \\(y_{max}\\), \\(z_{min}\\), \\(z_{max}\\). For each edge, 1 or 2 lines are used to define the boundary condition.\nIn the case of a “GMSH” mesh: code \\(3\\)\nThis is the case where the MAME section has code \\(2\\). The “zone” numbers defined by GMSH are then used. The argument defines the number of zones to process. These zones correspond to the domain edges. For each zone, 1 or 2 lines are used to define the boundary condition.\n\n\n\n\n\n\n\nTipExample boundary conditions\n\n\n\nCOND    1    0\n1\n1\n1\n1\nIn this example, “mirror” conditions are applied on the 4 edges of a rectangular “2D” domain.\n\n\n\n\n\n2.2.7 Defining the mesh: MESH\nThe MESH header is used to define the domain mesh, i.e. the discretization of the blocks defined in the MAME section. In addition, the mesh (de)refinement parameters are defined here. The default code is \\(0\\).\nThe following keywords can then be defined/modified:\n\nNBDS: Number of domains to create, i.e. number of MPI processes (default value \\(1.\\))\nMARL: Maximum refinement level (“MAximal Refinement Level”) (default value \\(0.\\))\nCOPA: Threshold value of the coarsening criterion (“mesh COarsening PArameter”) (default value \\(0.002\\))\nREPA: Threshold value of the refinement criterion (“mesh REfinement PArameter”) (default value \\(0.02\\))\nFDRA: Number of the user-defined function specifying the refinement to apply. This function is defined in the routine ~/CERF/sources/UTI/USER/fct_mesh.f90\nNZRA: Number of zones defining the refinement to apply (default 1.). Followed, for each zone, by \\(1\\) integer and \\(6\\) reals defining the refinement level to apply to each block in the zone and the zone boundaries \\(x_{min}\\), \\(x_{max}\\), \\(y_{min}\\), \\(y_{max}\\), \\(z_{min}\\), \\(z_{max}\\).\n\n\n\n\n\n\n\nTip 2.2: Example of refinement by zones\n\n\n\nMESH 0\nNBDS 4\nMARL 3\nCOPA 0.8\nREPA 1.\nNZRA 3\n0 -1. 4.1 -1. 3.1 -1. 1.\n3 -1. 1.25 -1. 2.25 -1. 1.\n0 -1. .75 -1. 1.75 -1. 1.\nWe return to the “dam-break” example Tip 2.2. In a “2D” domain of \\(4 \\times 3\\) filled with air, we consider a column of water on the left of height \\(2\\) and width \\(1\\). The objective is to refine the mesh around the air-water interface. We therefore define 3 zones. The first zone broadly covers the entire domain and is assigned refinement level \\(0\\). The second covers the domain occupied by the water with a margin of \\(0.25\\), to which the maximum refinement level \\(3\\) is assigned. And finally, the third covers the domain occupied by the water with a margin of \\(-0.25\\), to which the minimum refinement level \\(0\\) is assigned, in order to minimize the number of cells. CERF will automatically (i.e. without user intervention) adjust the mesh level between levels \\(0\\) and \\(3\\) in order to satisfy the criteria imposed by the mesh generator. In addition, the domain is decomposed into 4 domains, and the mesh refinement criterion management parameters are modified Section 1.1.3.3.\n\n\n\n\n\n2.2.8 Defining “obstacles”: OBST\nThe OBST header is used to define a set of cells to which a “rigid fluid” penalization is applied, naturally in the case where the two-fluid model is selected. The default code is \\(0\\), and the argument defines the number of obstacles to process.\nFor each obstacle, the following is written on one line: a file name, the obstacle density, and the penalization type.\n\nfile name: this is a character string (max 20 characters) defining the name of the ASCII file containing the mesh of the obstacle envelope, with triangular elements in .OBJ format. This file must be located in the exec directory.\nobstacle density: this is a real number defining the density of the obstacle. This density must satisfy \\(\\rho_a \\leq \\rho_{obs} \\leq \\rho_w\\).\npenalization type: this is an integer defining the type of penalization. The available penalization types are:\n\n-1: “fixed” type penalization. A zero velocity is imposed.\n0: “free” type penalization. The obstacle behaves as a free rigid solid.\n&gt;0: “imposed” type penalization. The obstacle behaves as a rigid solid for which the velocity of the center of gravity and the rotational velocity are imposed. The integer value defines the number of the user-defined function describing the motion in the routine ~/CERF/sources/PHY/USER/fct_mvsol.f90.\n\n\n\n\n\n\n\n\nTipExample obstacle file in .obj format\n\n\n\n\n\nv 0. 0. 0.\nv 1. 0. 0.\nv 1. 1. 0.\nv 0. 1. 0.\nv 0. 0. 1.\nv 1. 0. 1.\nv 1. 1. 1.\nv 0. 1. 1.\nf 1 2 3\nf 1 3 4 \nf 2 6 7\nf 2 7 3\nf 6 5 7\nf 5 8 7 \nf 5 1 4\nf 5 4 8\nf 3 4 7\nf 4 7 8\nf 2 5 6\nf 2 1 5\nIn this example, we define a cube-shaped obstacle with side length \\(1\\), centered at \\((0.5,0.5,0.5)\\). The vertices of the obstacle are defined, followed by its faces. Vertices are defined by the letter v followed by the vertex coordinates. Faces are defined by the letter f followed by the indices of the face’s vertices. The figure below illustrates the shape of the obstacle thus defined. \n\n\n\n\n\n\n2.2.9 Defining the numerical parameters: NUME\nThe NUME header is used to define the numerical parameters of the problem to be solved. The default code is \\(0\\).\nThe following keywords can then be defined/modified:\n\nTINT: Computation time interval (Time INTerval) (default value \\(0.\\))\nCCFL: CFL coefficient Important 1.1 (default value \\(0.9\\))\nTDOR: Time discretization order (Time Discretization ORder) 1 or 2 (default value \\(1.\\))\nSDOR: Space discretization order (Space Discretization ORder) 1 or 2 (default value \\(1.\\))\nLIMI: Limiter type 0/Barth, 1/Cartesian (default value \\(0.\\))\nNPRO: Number of probes &lt; 10 (default value \\(0.\\)). Then, for each probe, one line containing 3 reals (coordinates \\(x\\), \\(y\\), \\(z\\) of the probe)\nSAVE: Save type 0/No, 1/Yes. The units digit refers to saving in cerfbin-xxx files. The tens digit refers to saving in 1D gnuplot format. (default value \\(01.\\))",
    "crumbs": [
      "<span class='chapter-number'>2</span>  <span class='chapter-title'>User Guide</span>"
    ]
  },
  {
    "objectID": "chapitre2.html#the-cerf-source-files",
    "href": "chapitre2.html#the-cerf-source-files",
    "title": "2  User Guide",
    "section": "2.3 The CERF source files",
    "text": "2.3 The CERF source files\n\n2.3.1 Structure and modules\nAs illustrated in Figure 2.1, the CERF source code is organized into several directories:\n\n./sources/PHY contains all the lowest-level routines for handling physical variables and entities. The global physical variables and the data types (in the Fortran 90 sense), such as the degrees of freedom, their gradients, the fluxes, or the physical characteristics, are defined in the module phy_typ (./sources/PHY/mod_phy_typ.f90). The module phy (./sources/PHY/mod_phy.f90), in turn, contains the module phy_typ as well as all the interfaces of the routines contained in the PHY directory. During compilation, the compiled result of the routines in this directory is archived in the library ./lib/lib_PHY.a.\n./sources/GEO contains all the routines for handling geometric variables and entities: computing interface fluxes, computing source terms, etc. The global geometric variables and the data types (in the Fortran 90 sense), such as the list of faces, the list of cells, the mesh, etc., are defined in the module geo_typ (./sources/GEO/mod_geo_typ.f90). The module geo_typ contains the module phy. They are nested. The module geo (./sources/GEO/mod_geo.f90), in turn, contains the module geo_typ as well as all the interfaces of the routines contained in the GEO directory. During compilation, the compiled result of the routines in this directory is archived in the library ./lib/lib_GEO.a.\n./sources/NUM contains all the higher-level routines for handling numerical variables and entities: the rk2 algorithm, saving, MPI communication, etc. The global numerical variables and the data types (in the Fortran 90 sense), such as the “computation”, are defined in the module num_typ (./sources/NUM/mod_num_typ.f90). The module num_typ contains the module geo. They are nested. The module num (./sources/NUM/mod_num.f90), in turn, contains the module num_typ as well as all the interfaces of the routines contained in the NUM directory. During compilation, the compiled result of the routines in this directory is archived in the library ./lib/lib_NUM.a.\n./sources/UTI contains the programs and all the routines required for pre- and post-processing, and therefore in particular for managing the AMR mesh. The module uti (./sources/UTI/mod_uti.f90) contains the module num as well as all the interfaces of the routines contained in the UTI directory. During compilation, the compiled result of the routines in this directory is archived in the library ./lib/lib_UTI.a.\n\nThe modules phy, geo, num, and uti are nested and interdependent. The data types defined in them are the basic tools of the finite-volume method. As part of managing AMR meshes, it was necessary to develop new tools, which have been grouped into two modules.\n\nzorder (./sources/UTI/mod_zorder.f90): this module contains the tools for managing the “z-order” data structure used to manage the AMR mesh. It contains all the routines required for handling Morton codes.\nmsh (./sources/UTI/mod_msh.f90): this module contains the tools for managing the “mesh” data structure used to manage the AMR mesh. It contains all the routines required for managing the AMR mesh. The data structures of the finite-volume code are stored in “fixed” arrays for performance reasons. Managing mesh graphs, on the other hand, is more complex and requires dynamic data structures with linked lists. The msh module contains all the tools required for managing these data structures.\n\n\n\n2.3.2 Dependency graph of the main program\nBelow is the dependency graph of the main program cerf (excluding MPI routines and display routines). Routines depending on the phy module are outlined in green, routines depending on the geo module are outlined in blue, and routines depending on the num module are outlined in red. The HTML documentation provides direct access to a routine’s source code by clicking on its name in the graph Figure 2.3.\n\n\n\n\n\n\nflowchart LR\n    A[cerf] --&gt; B(init_calc)\n    B --&gt; AA(empty_mesh)\n    A --&gt; E(rk2)\n    A --&gt; C(save_1d)\n    A --&gt; D(save_calc)\n    E --&gt; F(dtloc)\n    E --&gt; G(synchro_mpi)\n    E --&gt; H(gradprim)\n    E --&gt; I(gradlim)\n    E --&gt; J(gradlimxy)\n    E --&gt; K(synchro_w)\n    E --&gt; L(entropie)\n    E --&gt; M(calc_flux)\n    E --&gt; N(calc_source)\n    E --&gt; O(friction)\n    E --&gt; P(sharp)\n    E --&gt; Q(u_penalisation)\n    F --&gt; R(velcf)\n    R --&gt; T(vitson)\n    M --&gt; U(wext)\n    M --&gt; V(prim2bal)\n    M --&gt; W(bal2prim)\n    M --&gt; X(fluxnum)\n    M --&gt; Y(fluxnumnoncons)\n    Q --&gt; Z(raycasting)\n    U --&gt; W\n    U --&gt; V\n    C --&gt; W\n    C --&gt; AB(riemisot)\n    C --&gt; AC(svriemann)\n    H --&gt; U\n    I --&gt; U\n    J --&gt; U\n    B --&gt; S(get_calc)\n    N --&gt; AD(sources)\n    O --&gt; V\n    P --&gt; W\n    Q --&gt; V\n    Q --&gt; AE(fct_mvsol)\n\n    classDef num  color:#000000,stroke:#FF0000,stroke-width:4px\n    classDef geo  color:#000000,stroke:#0000FF,stroke-width:4px\n    classDef phy  color:#000000,stroke:#00FF00,stroke-width:4px\n    classDef uti  color:#000000,stroke:#800080,stroke-width:4px\n    class A,B,C,D,E,F,G,K,M,N,O,P,Q,S,Z num\n    class H,I,J,AA geo\n    class L,U,V,W,X,Y,AB,AC,R,T,AD,AE phy\n    %%class G uti\n\n  click A href \"../robodoccerf/NUM/cerf_f90.html\" _blank\n  click B href \"../robodoccerf/NUM/init_calc_f90.html\" _blank\n  click E href \"../robodoccerf/NUM/rk2_f90.html\" _blank\n  click C href \"../robodoccerf/NUM/save_1d_f90.html\" _blank\n  click D href \"../robodoccerf/NUM/save_calc_f90.html\" _blank\n  click F href \"../robodoccerf/NUM/dtloc_f90.html\" _blank\n  click G href \"../robodoccerf/GEO/synchro_mpi_f90.html\" _blank\n  click H href \"../robodoccerf/GEO/gradprim_f90.html\" _blank\n  click I href \"../robodoccerf/GEO/gradlim_f90.html\" _blank\n  click J href \"../robodoccerf/GEO/gradlimxy_f90.html\" _blank\n  click K href \"../robodoccerf/NUM/synchro_w_f90.html\" _blank\n  click L href \"../robodoccerf/PHY/entropie_f90.html\" _blank\n  click M href \"../robodoccerf/NUM/calc_flux_f90.html\" _blank\n  click N href \"../robodoccerf/NUM/calc_source_f90.html\" _blank\n  click O href \"../robodoccerf/NUM/friction_f90.html\" _blank\n  click P href \"../robodoccerf/NUM/sharp_f90.html\" _blank\n  click Q href \"../robodoccerf/NUM/u_penalisation_f90.html\" _blank\n  click R href \"../robodoccerf/PHY/velcf_f90.html\" _blank\n  click S href \"../robodoccerf/NUM/get_calc_f90.html\" _blank\n  click T href \"../robodoccerf/PHY/vitson_f90.html\" _blank\n  click U href \"../robodoccerf/PHY/wext_f90.html\" _blank\n  click V href \"../robodoccerf/PHY/prim2bal_f90.html\" _blank\n  click W href \"../robodoccerf/PHY/bal2prim_f90.html\" _blank\n  click X href \"../robodoccerf/PHY/fluxnum_f90.html\" _blank\n  click Y href \"../robodoccerf/PHY/fluxnumnoncons_f90.html\" _blank\n  click Z href \"../robodoccerf/NUM/raycasting_f90.html\" _blank\n  click AA href \"../robodoccerf/GEO/empty_mesh_f90.html\" _blank\n  click AB href \"../robodoccerf/PHY/riemisot_f90.html\" _blank\n  click AC href \"../robodoccerf/PHY/svriemann_f90.html\" _blank\n  click AD href \"../robodoccerf/PHY/sources_f90.html\" _blank\n  click AE href \"../robodoccerf/PHY/USER/fct_mvsol_f90.html\" _blank\n\n\n\n\n\n\nFigure 2.3: Dependency graph of the cerf routines\n\n\n\n\n\n\n\n2.3.3 Dependency graph of the pre-processing program\nBelow is the dependency graph of the pre-processing program cerf_input (excluding MPI routines and display routines). Routines depending on the phy module are outlined in green, routines depending on the geo module are outlined in blue, routines depending on the num module are outlined in red, and routines depending on the uti module are outlined in purple. The HTML documentation provides direct access to a routine’s source code by clicking on its name in the graph Figure 2.4.\n\n\n\n\n\n\nflowchart LR\n    A[cerf_input] --&gt; B(read_phys)\n    A --&gt; C(read_cond)\n    A --&gt; D(read_init)\n    A --&gt; E(read_bath)\n    A --&gt; F(read_obst)\n    A --&gt; G(read_mame)\n    A --&gt; H(read_nume)\n    A --&gt; I(read_mesh)\n    A --&gt; J(block2mesh)\n    A --&gt; K(save_calc)\n    A --&gt; L(amr2tec)\n    J --&gt; P(init_def_bloc)\n    J --&gt; R(splitdom_Z)\n    J --&gt; S(inter)\n    J --&gt; T(empty_mesh)\n    J --&gt; U(initialise)\n    J --&gt; V(mesh1d)\n    J --&gt; W(mesh2d)\n    J --&gt; X(mesh3d)\n    J --&gt; Y(init_mesh)\n    J --&gt; Z(trivertex)\n    J --&gt; AA(ini_w)\n    J --&gt; K\n    J --&gt; AC(copy_sol)\n    C --&gt; AD(prim2bal)\n    C --&gt; AE(bal2prim)\n    D --&gt; AD\n    F --&gt; M(read_obj)\n    G --&gt; N(mame2bloc)\n    G --&gt; O(read_gmsh)\n    P --&gt; Q(fct_mesh)\n    V --&gt; AF(psr)\n    V --&gt; AG(pni)\n    V --&gt; AH(volcel)\n    W --&gt; AF\n    W --&gt; AG\n    W --&gt; AH\n    X --&gt; AF\n    X --&gt; AG\n    X --&gt; AH\n    Y --&gt; AI(triangle)\n    AA --&gt; AJ(fct_iniw)\n    AA --&gt; AK(fct_bathy)\n    AC --&gt; AL(empty_sol)\n    N --&gt; ZZ((Module zorder))\n    O --&gt; ZZ\n    O --&gt; MM((Module msh))\n    Z --&gt; ZZ\n    classDef num  color:#000000,stroke:#FF0000,stroke-width:4px\n    classDef geo  color:#000000,stroke:#0000FF,stroke-width:4px\n    classDef phy  color:#000000,stroke:#00FF00,stroke-width:4px\n    classDef uti  color:#000000,stroke:#800080,stroke-width:4px\n    class B,C,D,E,H,I,K,L,N,R,S,AA,AB,AC num\n    class M,T,U,V,W,X,Y,Z,AF,AG,AH,AL geo\n    class Q,AD,AE,AJ,AK phy\n    class A,F,G,J,O,P uti\n\n  click A href \"../robodoccerf/UTI/PRE/cerf_input_f90.html\" _blank\n  click B href \"../robodoccerf/UTI/PRE/read_phy_f90.html\" _blank\n  click C href \"../robodoccerf/UTI/PRE/read_cond_f90.html\" _blank\n  click D href \"../robodoccerf/UTI/PRE/read_init_f90.html\" _blank\n  click E href \"../robodoccerf/UTI/PRE/read_bath_f90.html\" _blank\n  click F href \"../robodoccerf/UTI/PRE/read_obst_f90.html\" _blank\n  click G href \"../robodoccerf/UTI/PRE/read_mame_f90.html\" _blank\n  click H href \"../robodoccerf/UTI/PRE/read_nume.html\" _blank\n  click I href \"../robodoccerf/UTI/PRE/read_mesh_f90.html\" _blank\n  click J href \"../robodoccerf/UTI/PRE/block2mesh_f90.html\" _blank\n  click K href \"../robodoccerf/NUM/save_calc_f90.html\" _blank\n  click L href \"../robodoccerf/UTI/AMR/amr2tec_f90.html\" _blank\n  click M href \"../robodoccerf/UTI/PRE/read_obj_f90.html\" _blank\n  click N href \"../robodoccerf/UTI/PRE/mame2bloc_f90.html\" _blank\n  click O href \"../robodoccerf/UTI/PRE/read_gmsh_f90.html\" _blank\n  click P href \"../robodoccerf/UTI/PRE/init_def_bloc_f90.html\" _blank\n  click Q href \"../robodoccerf/UTI/USER/fct_mesh_f90.html\" _blank\n  click R href \"../robodoccerf/UTI/PRE/splitdom_Z_f90.html\" _blank\n  click S href \"../robodoccerf/UTI/PRE/inter_f90.html\" _blank\n  click T href \"../robodoccerf/GEO/empty_mesh_f90.html\" _blank\n  click U href \"../robodoccerf/UTI/PRE/initialise_f90.html\" _blank\n  click V href \"../robodoccerf/GEO/mesh1d_f90.html\" _blank\n  click W href \"../robodoccerf/GEO/mesh2d_f90.html\" _blank\n  click X href \"../robodoccerf/GEO/mesh3d_f90.html\" _blank\n  click Y href \"../robodoccerf/UTI/PRE/init_mesh_f90.html\" _blank\n  click Z href \"../robodoccerf/UTI/PRE/trivertex_f90.html\" _blank\n  click AA href \"../robodoccerf/UTI/PRE/ini_w_f90.html\" _blank\n  click AC href \"../robodoccerf/UTI/PRE/copy_sol_f90.html\" _blank\n  click AD href \"../robodoccerf/PHY/prim2bal_f90.html\" _blank\n  click AE href \"../robodoccerf/PHY/bal2prim_f90.html\" _blank\n  click AF href \"../robodoccerf/GEO/fonctions_f90.html\" _blank\n  click AG href \"../robodoccerf/GEO/fonctions_f90.html\" _blank\n  click AH href \"../robodoccerf/GEO/fonctions_f90.html\" _blank\n  click AI href \"../robodoccerf/GEO/triangle_f90.html\" _blank\n  click AJ href \"../robodoccerf/PHY/USER/fct_iniw_f90.html\" _blank\n  click AK href \"../robodoccerf/PHY/USER/fct_bathy_f90.html\" _blank\n  click AL href \"../robodoccerf/UTI/PRE/empty_sol_f90.html\" _blank\n  click MM href \"../robodoccerf/UTI/mod_msh_f90.html\" _blank\n  click ZZ href \"../robodoccerf/UTI/mod_zorder_f90.html\" _blank\n\n\n\n\n\n\nFigure 2.4: Dependency graph of the cerf_input routines\n\n\n\n\n\n\n\n2.3.4 Dependency graph of the mesh-modification program\nBelow is the dependency graph of the mesh-modification program cerf_amr (excluding MPI routines and display routines). Routines depending on the phy module are outlined in green, routines depending on the geo module are outlined in blue, routines depending on the num module are outlined in red, and routines depending on the uti module are outlined in purple. The HTML documentation provides direct access to a routine’s source code by clicking on its name in the graph Figure 2.5.\n\n\n\n\n\n\nflowchart LR\n    A[cerf_input] --&gt; B(get_calc)\n    A --&gt; T(empty_mesh)\n    A --&gt; AL(empty_sol)\n    A --&gt; C(critraf)\n    C --&gt; D(critphy)\n    A --&gt; K(save_calc)\n    A --&gt; L(amr2tec)\n    A --&gt; R(splitdom_Z)\n    A --&gt; S(inter)\n    A --&gt; AC(copy_sol)\n    A --&gt; V(mesh1d)\n    A --&gt; W(mesh2d)\n    A --&gt; X(mesh3d)\n    A --&gt; Y(init_mesh)\n    Y --&gt; AI(triangle)\n    A --&gt; Z(trivertex)\n    Z --&gt; ZZ((Module zorder))\n    A --&gt; E(proddl)\n    E --&gt; AK(fct_bathy)\n    E --&gt; AD(prim2bal)\n    C --&gt; AD(prim2bal)\n    V --&gt; AF(psr)\n    V --&gt; AG(pni)\n    V --&gt; AH(volcel)\n    W --&gt; AF\n    W --&gt; AG\n    W --&gt; AH\n    X --&gt; AF\n    X --&gt; AG\n    X --&gt; AH\n    A --&gt; F(raycasting)\n    A --&gt; G(save_obj)\n    classDef num  color:#000000,stroke:#FF0000,stroke-width:4px\n    classDef geo  color:#000000,stroke:#0000FF,stroke-width:4px\n    classDef phy  color:#000000,stroke:#00FF00,stroke-width:4px\n    classDef uti  color:#000000,stroke:#800080,stroke-width:4px\n    class B,C,K,L,R,S,AC,E,F num\n    class T,AL,V,W,X,AF,AG,AH,Y,Z,G geo\n    class D,AK,AD phy\n    class A uti\n\n  click A href \"../robodoccerf/UTI/AMR/cerf_amr_f90.html\" _blank\n  click T href \"../robodoccerf/GEO/empty_mesh_f90.html\" _blank\n  click B href \"../robodoccerf/NUM/get_calc_f90.html\" _blank\n  click AL href \"../robodoccerf/UTI/PRE/empty_sol_f90.html\" _blank\n  click C href \"../robodoccerf/UTI/AMR/critraf_f90.html\" _blank\n  click D href \"../robodoccerf/PHY/critphy_f90.html\" _blank\n  click K href \"../robodoccerf/NUM/save_calc_f90.html\" _blank\n  click L href \"../robodoccerf/UTI/AMR/amr2tec_f90.html\" _blank\n  click R href \"../robodoccerf/UTI/PRE/splitdom_Z_f90.html\" _blank\n  click S href \"../robodoccerf/UTI/PRE/inter_f90.html\" _blank\n  click AC href \"../robodoccerf/UTI/PRE/copy_sol_f90.html\" _blank\n  click V href \"../robodoccerf/GEO/mesh1d_f90.html\" _blank\n  click W href \"../robodoccerf/GEO/mesh2d_f90.html\" _blank\n  click X href \"../robodoccerf/GEO/mesh3d_f90.html\" _blank\n  click AF href \"../robodoccerf/GEO/fonctions_f90.html\" _blank\n  click AG href \"../robodoccerf/GEO/fonctions_f90.html\" _blank\n  click AH href \"../robodoccerf/GEO/fonctions_f90.html\" _blank\n  click Y href \"../robodoccerf/UTI/PRE/init_mesh_f90.html\" _blank\n  click AI href \"../robodoccerf/GEO/triangle_f90.html\" _blank\n  click Z href \"../robodoccerf/UTI/PRE/trivertex_f90.html\" _blank\n  click E href \"../robodoccerf/NUM/proddl_f90.html\" _blank\n  click AK href \"../robodoccerf/PHY/USER/fct_bathy_f90.html\" _blank\n  click AD href \"../robodoccerf/PHY/prim2bal_f90.html\" _blank\n  click F href \"../robodoccerf/NUM/raycasting_f90.html\" _blank\n  click G href \"../robodoccerf/NUM/save_obj_f90.html\" _blank\n  click ZZ href \"../robodoccerf/UTI/mod_zorder_f90.html\" _blank\n\n\n\n\n\n\nFigure 2.5: Dependency graph of the cerf_amr routines",
    "crumbs": [
      "<span class='chapter-number'>2</span>  <span class='chapter-title'>User Guide</span>"
    ]
  },
  {
    "objectID": "chapitre3.html",
    "href": "chapitre3.html",
    "title": "3  Documented examples",
    "section": "",
    "text": "3.1 Examples with the Saint-Venant model\nIn this chapter, we present the implementation of a few classic flow test cases for the Saint-Venant model. The cases covered are as follows:",
    "crumbs": [
      "<span class='chapter-number'>3</span>  <span class='chapter-title'>Documented examples</span>"
    ]
  },
  {
    "objectID": "chapitre3.html#examples-with-the-saint-venant-model",
    "href": "chapitre3.html#examples-with-the-saint-venant-model",
    "title": "3  Documented examples",
    "section": "",
    "text": "Solving a Riemann problem\nHydraulic jump\nFriction, McDonald test case\n…\n\n\n3.1.1 Riemann problem\nWe consider the classic 1D dam-break test case for the Saint-Venant model on a flat bottom. We use the test case from Faccanoni and Mangeney (2013) (test case 5). The domain is defined by \\(x\\in[-2.,2.]\\) and \\(t\\in[0,0.5]\\). The initial conditions are given by: \\[\n\\overrightarrow{w_p}(x,0) =\\left\\{\n\\begin{array}{ll}\n\\overrightarrow{w_p}_L = &lt;h_L,u_L&gt;^T = &lt;0.0046,0.&gt;^T &\\mbox{ if } x&lt;0 \\\\\n\\overrightarrow{w_p}_R = &lt;h_R,u_R&gt;^T = &lt;0.1446,0.&gt;^T &\\mbox{ if } x&gt;0\n\\end{array}\n\\right.\n\\]\nThis is therefore a dam-break case without a dry zone, with a shock and a rarefaction wave, for which the exact solution can be computed using the Riemann solver.\nFirst, we consider a uniform, fixed 1D mesh of 400 cells over a single domain. The data file corresponding to this test case is provided below:\n\n\n\n\n\n\nTip 3.1: Data file ex1_sv_Riemann.inp\n\n\n\n\n\n!-----------------------------------------------------------------------\n! Data file in CERF format\n! Each block of data is characterized in the first line by a header (4 characters)\n! then a code (1 integer) and an argument (1 integer)\n! PHYS MAME INIT BATH COND MESH OBST NUME\n!-----------------------------------------------------------------------\n\n!-----------------------------------------------------------------------\n! Reading physical parameters \n! Header PHYS\n! 1 keyword (4 characters) and a code (1 integer)\n! for now we have the same values for all domains\n!-----------------------------------------------------------------------\n! MODE : physical model      \n!        1/ Shallow water\n!        2/ Euler isothermal two-phases               (default 0.)\n!---- for isothermal Euler model\n! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)\n! PREF : isothermal pressure of reference             (default 1.d5)\n! CSSM : Mixture sound velocity (isothermal)          (default 20.)\n! REFA : Isothermal reference density for Air         (default 1.)\n! REFW : Isothermal reference density for Water       (default 1000.)\n! SHAR : isothermal interface sharpening parameter \n!        recommended value 0.01                       (defaut 0.)\n!---- for Shallow water model\n! WLTR : Value to truncate the Water level            (default 0.0001)\n! FRMO : Friction Model\n! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)\n! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation\n!-----------------------------------------------------------------------  \nPHYS 0\nMODE 1.\n\n!-----------------------------------------------------------------------\n! Reading the MAster MEsh\n! Header  MAME \n!-----------------------------------------------------------------------\n! code 1 -&gt; creation of a mesh\n!\n!           argument: 0 type \"1D shock tube\"\n! next line: xmin,xmax\n! next line: nx\n!\n!           argument: 1 type \"2d box\"\n! next line: xmin,xmax,ymin,ymax\n! next line: nx,ny\n!\n!           argument: 2 type \"3d box\"\n! next line: xmin,xmax,ymin,ymax,zmin,zmax\n! next line: nx,ny,nz\n!\n! code 2 -&gt; Reading an ascii file format GMSH V2.2\n! next line: file name\n!-----------------------------------------------------------------------\nMAME 1 0\n-2. 2.  \n400\n\n!-----------------------------------------------------------------------\n! Reading Initial conditions\n! Header INIT \n!-----------------------------------------------------------------------\n! code 0 -&gt; Definition for a shock tube\n!           primitive variables on left  \n!           primitive variables on right\n!           argument: 0\n! code 1 -&gt; user function \n!           argument: Function's Number \n! code 2 -&gt; Definition using areas\n!           argument: number of areas\n!           and for each area 2 lignes\n!           - xmin,xmax,ymin,ymax,zmin,zmax\n!           - nvar values ( primitive variables)\n! Shallow water model: water level, x velocity, y velocity\n! Air-water flow model: density, x velocity, y velocity,z velocity, pressure\n!-----------------------------------------------------------------------\nINIT 0 0\n0.0046 0. 0. \n0.1446 0. 0.  \n\n!-----------------------------------------------------------------------\n! Reading bathymetry and friction\n! Header BATH \n!-----------------------------------------------------------------------\n! code 0 -&gt; Default : flat bottom and no friction\n\n! code 1 -&gt; bathymetry and friction defined by user function (defined in sources/PHY/fct_bathy.f90)\n!           argument: Function's Number \n! code 2 -&gt; bathymetry and friction defined linearly per block \n!           argument: 0  reading the file \n!           argument:-1  creation of a file containing the coordinates of the vertices of the block \n! next line: file name\n!-----------------------------------------------------------------------\nBATH    0   0\n\n!-----------------------------------------------------------------------\n! Reading boundary condition\n! (maximum 10 for the moment)\n! Header COND \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; Description in the case  \"shock tube\" : Nothing\n!\n! code 1 -&gt; Description in the case  \"2d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax\n!\n! code 2 -&gt; Description in the case  \"3d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax\n!\n! code 3 -&gt; Definition per zone for  GMSH mesh\n!           argument: number of zones defining a boundary condition\n!           number of zone, kind of boundary condition \n!\n! kind of boundary condition:\n! 0 outlet (we copy)\n! 1 miror\n! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values\n! 3 user's function , then next line function number\n! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),\n!   then next line nvar values\n! 999 Non used\n!-----------------------------------------------------------------------\nCOND    0    0\n\n!-----------------------------------------------------------------------\n! Reading mesh parameters in order to define the mesh according to\n! the Master mesh is the one represented by the domain \"0\".\n! Header MESH \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; default value\n!\n! 1 keyword (4 characters) and a real\n! NBDS : NumBer of DomainS                         (default 001.)\n! MARL : MAximal Refinement Level                  (default 000.)\n! COPA : Mesh COarsening PArameter          0&lt;..&lt;1 (default 0.002)\n! REPA : Mesh REfinement PArameter          0&lt;..&lt;1 (default 0.02)\n! FDRA : Function defining the mesh refinement to be applied \n! NZRA : Number of zones where initial blocks are refined (default 1.)\n!        then for each zone\n!        nrb,xmin,xmax,ymin,ymax,zmin,zmax\n!-----------------------------------------------------------------------\nMESH    0 \nNBDS 1\nMARL 0\nNZRA 1\n0 -2.  2. -100. 100. -100. 100.\n\n!-----------------------------------------------------------------------\n! Reading numerical parameter \n! Header  NUME\n!-----------------------------------------------------------------------\n! 1 keyword (4 characters) then a real\n! code 0 -&gt; default value\n!\n! TINT : Time interval to compute                    (default 0.)\n! CCFL : cfl condition                               (default 0.9)\n! TDOR : Time discretization order 1 ou 2         (defaut 1.)\n! SDOR : Space discretization order 1 ou 2        (defaut 1.)\n! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)\n! SAVE : kind of backup, three-digit number          (default 001.)\n!        decade  : shock tube backup 0/no 1/yes \n!        unit    : binary backup     0/no 1/yes \n! NPRO : number of probe  &lt; 10                      (defaut 000.)\n!        then coordinates x,y,z of the probe (1 probe per line)\n!-----------------------------------------------------------------------\nNUME 0\nTINT .5\nCCFL .9\nSDOR 2.\nTDOR 2.\nSAVE 11.\n\n\n\n\nTo run the computation, since there is only one domain to process over a single time interval:\n./cerf_input ex1_sv_Riemann.inp\n./cerf\nSince the data file ex1_sv_Riemann.inp requests output in 1D format (parameter SAVE=11, section NUME), 1D output files have been generated. Files for the numerical solution in water depth and velocity ‘n_h_001’ and ‘n_u_001’, and files for the exact solution in water depth and velocity: ‘e_h_001’ and ‘e_u_001’. These files can be visualized, for example, with gnuplot:\ngnuplot\ngnuplot&gt; set term x11\ngnuplot&gt; plot 'e_h_001' w l, 'n_h_001' w l\ngnuplot&gt; plot 'e_u_001' w l, 'n_u_001' w l\ngnuplot&gt; quit\nThe numerical solution can thus be compared to the exact solution.\n\n\n\n\nNumerical/exact comparison of water depth at \\(t=0.5s\\), uniform mesh of 400 cells\n\n\n\n\n\n\n\nNumerical/exact comparison of velocity at \\(t=0.5s\\), uniform mesh of 400 cells\n\n\n\nTo test the code with a “dry” zone (Faccanoni and Mangeney (2013), test case 1), the same problem can for example be revisited with a different initial condition by modifying the INIT section of the data file ex1_sv_Riemann.inp:\nINIT    0   0\n0. 0. 0. \n0.1446 0. 0.  \nThis then gives:\n\n\n\n\nNumerical/exact comparison of water depth at \\(t=0.5s\\), uniform mesh of 400 cells with a dry zone\n\n\n\n\n\n\n\nNumerical/exact comparison of velocity at \\(t=0.5s\\), uniform mesh of 400 cells with a dry zone\n\n\n\n\n\n3.1.2 Hydraulic jump\nWe consider here a steady-state test case, without friction, with a non-flat bottom. The analytical solution is taken from Delestre et al. (2013). On a domain \\(x\\in[0.,25.]\\), we consider a topography defined by:\n\\[z_b(x) = \\left\\{ \\begin{array}{ll} 0.2 - 0.05 \\left( x - 10 \\right)^2 & \\mbox{ if } 8 &lt; x &lt; 12 \\\\ 0. & \\mbox{ otherwise }  \\end{array}\\right. \\]\nWe consider the case of transcritical flow with a shock, characterized by the following boundary conditions:\n\nOn the left, \\(h = 0.33 m\\) is imposed and the discharge is left free.\nOn the right, \\(hu = 0.18 m^2/s\\) is imposed and the water depth is left free.\n\nAs recommended by Delestre et al. (2013), the initial conditions are given by: \\(h + z_b =0\\) and \\(u=0\\).\nThe computation is carried out on a fixed, regular mesh of 2500 cells distributed over 8 domains. In order to reach a steady state, the computation is run up to \\(T=100s\\). The data file corresponding to this test case is provided below Tip 3.2:\n\n\n\n\n\n\nTip 3.2: Data file ex2_sv_Bump.inp\n\n\n\n\n\n!-----------------------------------------------------------------------\n! Data file in CERF format\n! Each block of data is characterized in the first line by a header (4 characters)\n! then a code (1 integer) and an argument (1 integer)\n! PHYS MAME INIT BATH COND MESH OBST NUME\n!-----------------------------------------------------------------------\n\n!-----------------------------------------------------------------------\n! Reading physical parameters \n! Header PHYS\n! 1 keyword (4 characters) and a code (1 integer)\n! for now we have the same values for all domains\n!-----------------------------------------------------------------------\n! MODE : physical model      \n!        1/ Shallow water\n!        2/ Euler isothermal two-phases               (default 0.)\n!---- for isothermal Euler model\n! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)\n! PREF : isothermal pressure of reference             (default 1.d5)\n! CSSM : Mixture sound velocity (isothermal)          (default 20.)\n! REFA : Isothermal reference density for Air         (default 1.)\n! REFW : Isothermal reference density for Water       (default 1000.)\n! SHAR : isothermal interface sharpening parameter \n!        recommended value 0.01                       (defaut 0.)\n!---- for Shallow water model\n! WLTR : Value to truncate the Water level            (default 0.0001)\n! FRMO : Friction Model\n! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)\n! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation\n!-----------------------------------------------------------------------  \nPHYS 0\nMODE 1.\nFRMO 0.\nBAGR 0.\n!-----------------------------------------------------------------------\n! Reading the MAster MEsh\n! Header  MAME \n!-----------------------------------------------------------------------\n! code 1 -&gt; creation of a mesh\n!\n!           argument: 0 type \"1D shock tube\"\n! next line: xmin,xmax\n! next line: nx\n!\n!           argument: 1 type \"2d box\"\n! next line: xmin,xmax,ymin,ymax\n! next line: nx,ny\n!\n!           argument: 2 type \"3d box\"\n! next line: xmin,xmax,ymin,ymax,zmin,zmax\n! next line: nx,ny,nz\n!\n! code 2 -&gt; Reading an ascii file format GMSH V2.2\n! next line: file name\n!-----------------------------------------------------------------------\nMAME    1    0\n0. 25.  \n2500\n!-----------------------------------------------------------------------\n! Reading Initial conditions\n! Header INIT \n!-----------------------------------------------------------------------\n! code 0 -&gt; Definition for a shock tube\n!           primitive variables on left  \n!           primitive variables on right\n!           argument: 0\n! code 1 -&gt; user function \n!           argument: Function's Number \n! code 2 -&gt; Definition using areas\n!           argument: number of areas\n!           and for each area 2 lignes\n!           - xmin,xmax,ymin,ymax,zmin,zmax\n!           - nvar values ( primitive variables)\n! Shallow water model: water level, x velocity, y velocity\n! Air-water flow model: density, x velocity, y velocity,z velocity, pressure\n!-----------------------------------------------------------------------\nINIT    1   1\n\n!-----------------------------------------------------------------------\n! Reading bathymetry and friction\n! Header BATH \n!-----------------------------------------------------------------------\n! code 0 -&gt; Default : flat bottom and no friction\n\n! code 1 -&gt; bathymetry and friction defined by user function (defined in sources/PHY/fct_bathy.f90)\n!           argument: Function's Number \n! code 2 -&gt; bathymetry and friction defined linearly per block \n!           argument: 0  reading the file \n!           argument:-1  creation of a file containing the coordinates of the vertices of the block \n! next line: file name\n!-----------------------------------------------------------------------\nBATH    1   1\n\n!-----------------------------------------------------------------------\n! Reading boundary condition\n! (maximum 10 for the moment)\n! Header COND \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; Description in the case  \"shock tube\" : Nothing\n!\n! code 1 -&gt; Description in the case  \"2d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax\n!\n! code 2 -&gt; Description in the case  \"3d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax\n!\n! code 3 -&gt; Definition per zone for  GMSH mesh\n!           argument: number of zones defining a boundary condition\n!           number of zone, kind of boundary condition \n!\n! kind of boundary condition:\n! 0 outlet (we copy)\n! 1 miror\n! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values\n! 3 user's function , then next line function number\n! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),\n!   then next line nvar values\n! 999 Non used\n!-----------------------------------------------------------------------\nCOND    2    0\n4\n1.e20 0.18 0. \n2\n0.33 1.e20 0. \n999\n999\n999\n999\n\n!-----------------------------------------------------------------------\n! Reading mesh parameters in order to define the mesh according to\n! the Master mesh is the one represented by the domain \"0\".\n! Header MESH \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; default value\n!\n! 1 keyword (4 characters) and a real\n! NBDS : NumBer of DomainS                         (default 001.)\n! MARL : MAximal Refinement Level                  (default 000.)\n! COPA : Mesh COarsening PArameter          0&lt;..&lt;1 (default 0.002)\n! REPA : Mesh REfinement PArameter          0&lt;..&lt;1 (default 0.02)\n! FDRA : Function defining the mesh refinement to be applied \n! NZRA : Number of zones where initial blocks are refined (default 1.)\n!        then for each zone\n!        nrb,xmin,xmax,ymin,ymax,zmin,zmax\n!-----------------------------------------------------------------------\nMESH 0\nNBDS 8\nNZRA 1\n0 0.  25. -100. 100. -100. 100.\n\n!-----------------------------------------------------------------------\n! Reading numerical parameter \n! Header  NUME\n!-----------------------------------------------------------------------\n! 1 keyword (4 characters) then a real\n! code 0 -&gt; default value\n!\n! TINT : Time interval to compute                    (default 0.)\n! CCFL : cfl condition                               (default 0.9)\n! TDOR : Time discretization order 1 ou 2         (defaut 1.)\n! SDOR : Space discretization order 1 ou 2        (defaut 1.)\n! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)\n! SAVE : kind of backup, three-digit number          (default 001.)\n!        decade  : shock tube backup 0/no 1/yes \n!        unit    : binary backup     0/no 1/yes \n! NPRO : number of probe  &lt; 10                      (defaut 000.)\n!        then coordinates x,y,z of the probe (1 probe per line)\n!-----------------------------------------------------------------------\nNUME 0\nTINT 100.\nCCFL .9\nSDOR 2.\nTDOR 2.\nSAVE 01.\n\n\n\n\n\nTo run the computation, since there is only one domain to process over a single time interval, but distributed over 8 domains, the following commands are used:\n./cerf_input ex2_sv_Bump.inp\nmpirun -np 8 ./cerf\n./cerf2tec\nThe analytical solution to this problem can be obtained from the Swashes website. In the doc directory of CERF, the analytical solution of \\(\\eta = h + z_b\\) on 2500 points in Tecplot format is provided for reference. The numerical solution can then be compared to the exact solution in the figure below (Figure 3.1).\n\n\n\n\n\n\n\nFigure 3.1: Numerical (red) / exact (black) comparison of water level at \\(t=100.s\\), uniform mesh of 2500 cells\n\n\n\n\nThe numerical solution can also be visualized with Visit or Paraview:\n\n\n\n\n\n\n\nFigure 3.2: Water level isovalues at \\(t=100.s\\), uniform mesh of 2500 cells\n\n\n\n\nIn figure Figure 3.2, the isovalues of the water level are shown at \\(t=100s\\). Note that the mesh has been modified for post-processing purposes. The coordinates of the nodes at the base of the hexahedra correspond to \\(z_b\\), while the nodes on the upper part of the hexahedra correspond to \\(z_b+h\\). This makes it possible to visualize the topography and the water level at the same time.\nThis computation can also be performed using automatic mesh refinement. To do this, the data file is modified accordingly ex2_sv_Bump.inp. A mesh consisting of only 250 blocks is then created, and mesh refinement is allowed only up to level 2. All other computation parameters are kept the same as for the fixed mesh.\n\n\n\n\n\n\nTip 3.3: Data file ex2_sv_Bump_amr.inp\n\n\n\n\n\n!-----------------------------------------------------------------------\n! Data file in CERF format\n! Each block of data is characterized in the first line by a header (4 characters)\n! then a code (1 integer) and an argument (1 integer)\n! PHYS MAME INIT BATH COND MESH OBST NUME\n!-----------------------------------------------------------------------\n\n!-----------------------------------------------------------------------\n! Reading physical parameters \n! Header PHYS\n! 1 keyword (4 characters) and a code (1 integer)\n! for now we have the same values for all domains\n!-----------------------------------------------------------------------\n! MODE : physical model      \n!        1/ Shallow water\n!        2/ Euler isothermal two-phases               (default 0.)\n!---- for isothermal Euler model\n! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)\n! PREF : isothermal pressure of reference             (default 1.d5)\n! CSSM : Mixture sound velocity (isothermal)          (default 20.)\n! REFA : Isothermal reference density for Air         (default 1.)\n! REFW : Isothermal reference density for Water       (default 1000.)\n! SHAR : isothermal interface sharpening parameter \n!        recommended value 0.01                       (defaut 0.)\n!---- for Shallow water model\n! WLTR : Value to truncate the Water level            (default 0.0001)\n! FRMO : Friction Model\n! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)\n! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation\n!----------------------------------------------------------------------- \nPHYS 0\nMODE 1.\n\n!-----------------------------------------------------------------------\n! Reading the MAster MEsh\n! Header  MAME \n!-----------------------------------------------------------------------\n! code 1 -&gt; creation of a mesh\n!\n!           argument: 0 type \"1D shock tube\"\n! next line: xmin,xmax\n! next line: nx\n!\n!           argument: 1 type \"2d box\"\n! next line: xmin,xmax,ymin,ymax\n! next line: nx,ny\n!\n!           argument: 2 type \"3d box\"\n! next line: xmin,xmax,ymin,ymax,zmin,zmax\n! next line: nx,ny,nz\n!\n! code 2 -&gt; Reading an ascii file format GMSH V2.2\n! next line: file name\n!-----------------------------------------------------------------------\nMAME    1    0\n0. 25.  \n2500\n!-----------------------------------------------------------------------\n! Reading Initial conditions\n! Header INIT \n!-----------------------------------------------------------------------\n! code 0 -&gt; Definition for a shock tube\n!           primitive variables on left  \n!           primitive variables on right\n!           argument: 0\n! code 1 -&gt; user function \n!           argument: Function's Number \n! code 2 -&gt; Definition using areas\n!           argument: number of areas\n!           and for each area 2 lignes\n!           - xmin,xmax,ymin,ymax,zmin,zmax\n!           - nvar values ( primitive variables)\n! Shallow water model: water level, x velocity, y velocity\n! Air-water flow model: density, x velocity, y velocity,z velocity, pressure\n!-----------------------------------------------------------------------\nINIT    1   1\n\n!-----------------------------------------------------------------------\n! Reading bathymetry and friction\n! Header BATH \n!-----------------------------------------------------------------------\n! code 0 -&gt; Default : flat bottom and no friction\n\n! code 1 -&gt; bathymetry and friction defined by user function (defined in sources/PHY/fct_bathy.f90)\n!           argument: Function's Number \n! code 2 -&gt; bathymetry and friction defined linearly per block \n!           argument: 0 \n! next line: file name\n! \n! code 3 -&gt; bathymetry and friction defined  per block at the highest mesh refinement level \n!           argument: 0 \n! next line: file name\n!-----------------------------------------------------------------------\nBATH    1   1\n\n!-----------------------------------------------------------------------\n! Reading boundary condition\n! (maximum 10 for the moment)\n! Header COND \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; Description in the case  \"shock tube\" : Nothing\n!\n! code 1 -&gt; Description in the case  \"2d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax\n!\n! code 2 -&gt; Description in the case  \"3d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax\n!\n! code 3 -&gt; Definition per zone for  GMSH mesh\n!           argument: number of zones defining a boundary condition\n!           number of zone, kind of boundary condition \n!\n! kind of boundary condition:\n! 0 outlet (we copy)\n! 1 miror\n! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values\n! 3 user's function , then next line function number\n! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),\n!   then next line nvar values\n! 999 Non used\n!-----------------------------------------------------------------------\nCOND    2    0\n4\n1.e20 0.18 0. 0. 0. 0.\n2\n0.33 1.e20 0. 0. 0. 0.\n999\n999\n999\n999\n\n!-----------------------------------------------------------------------\n! Reading mesh parameters in order to define the mesh according to\n! the Master mesh is the one represented by the domain \"0\".\n! Header MESH \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; default value\n!\n! 1 keyword (4 characters) and a real\n! NBDS : NumBer of DomainS                         (default 001.)\n! MARL : MAximal Refinement Level                  (default 000.)\n! COPA : Mesh COarsening PArameter          0&lt;..&lt;1 (default 0.002)\n! REPA : Mesh REfinement PArameter          0&lt;..&lt;1 (default 0.02)\n! FDRA : Function defining the mesh refinement to be applied \n! NZRA : Number of zones where initial blocks are refined (default 1.)\n!        then for each zone\n!        nrb,xmin,xmax,ymin,ymax,zmin,zmax\n!-----------------------------------------------------------------------\nMESH 0\nNBDS 8\nMARL 2\nNZRA 2\n0 0.  25. -100. 100. -100. 100.\n2 7.8  12.1 -100. 100. -100. 100.\n\n\n!-----------------------------------------------------------------------\n! Reading numerical parameter \n! Header  NUME\n!-----------------------------------------------------------------------\n! 1 keyword (4 characters) then a real\n! code 0 -&gt; default value\n!\n! TINT : Time interval to compute                    (default 0.)\n! CCFL : cfl condition                               (default 0.9)\n! TDOR : Time discretization order 1 ou 2         (defaut 1.)\n! SDOR : Space discretization order 1 ou 2        (defaut 1.)\n! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)\n! SAVE : kind of backup, three-digit number          (default 001.)\n!        decade  : shock tube backup 0/no 1/yes \n!        unit    : binary backup     0/no 1/yes \n! NPRO : number of probe  &lt; 10                      (defaut 000.)\n!        then coordinates x,y,z of the probe (1 probe per line)\n!-----------------------------------------------------------------------\nNUME 0\nTINT 10.\nCCFL .9\nSDOR 2.\nTDOR 2.\nSAVE 01.\n\n\n\n\n\nTo run the computation with mesh refinement, the following script is used:\n#!/bin/bash\n./cerf_input ex2_sv_Bump2.inp\n./cerf2tec\nmv cerfout.dat b_0.dat\nmv cerfamr.dat b_mame_0.dat\nfor i in $(seq 1 20 )\ndo\nmpirun -np 8 ./cerf\n./cerf2tec\nmv cerfout.dat b_$i.dat\n./cerf_amr\nmv cerfamr.dat b_mame_$i.dat\ndone\nAs before, the results can be compared with the exact solution.\n\n\n\n\n\n\n\nFigure 3.3: Numerical comparison on a fixed mesh (red) / numerical on a dynamic mesh (green) / exact (black) water level at \\(t=100.s\\)\n\n\n\n\nThe numerical solution can also be visualized with Visit or Paraview:\n\n\n\n\n\n\n\nFigure 3.4: Water level isovalues, dynamic mesh of 344 cells at \\(t=100.s\\)\n\n\n\n\nIn both Figure 3.3 and Figure 3.4, it can be seen that mesh refinement made it possible to reduce the number of cells while maintaining good accuracy in the numerical solution.\n\n\n3.1.3 Friction, McDonald test case\nWe consider here a steady-state test case, with friction, on a non-flat bottom. The semi-analytical solution is taken from Delestre et al. (2013) (case 3.2.1).\nOn a domain \\(x\\in[0.,1000.]\\), we consider a subcritical flow with a constant discharge at the inlet of \\(q = h u = 2 m^2/s\\) and a prescribed water depth at the outlet \\(h_{e}(1000)\\). The domain is initially dry, with \\(h = 0\\) and \\(u = 0\\). A Manning friction with \\(n = 0.033\\) is imposed.\nThe water depth is given by: \\[ h_{ex}(x) = \\left( \\frac{4}{g} \\right)^{1/3} \\left( 1 + \\frac{1}{2} exp \\left( - 16 \\left( \\frac{x}{1000} - \\frac{1}{2} \\right) \\right) \\right)\\]\nUsing Equation 1.6, in the steady-state case, the topography can then be recovered by integrating: \\[\\frac{\\partial z_b}{\\partial x} = \\left( \\frac{q^2}{g h^3} - 1 \\right) \\frac{\\partial h}{\\partial x} - S_f\\]\nThe computation is carried out on a fixed, regular mesh of 1000 cells. In order to reach a steady state, the computation is run up to \\(T=1500s\\). The data file corresponding to this test case is provided below Tip 3.4:\n\n\n\n\n\n\nTip 3.4: Data file ex3_sv_McDo.inp\n\n\n\n\n\n!-----------------------------------------------------------------------\n! Data file in CERF format\n! Each block of data is characterized in the first line by a header (4 characters)\n! then a code (1 integer) and an argument (1 integer)\n! PHYS MAME INIT BATH COND MESH OBST NUME\n!-----------------------------------------------------------------------\n\n!-----------------------------------------------------------------------\n! Reading physical parameters \n! Header PHYS\n! 1 keyword (4 characters) and a code (1 integer)\n! for now we have the same values for all domains\n!-----------------------------------------------------------------------\n! MODE : physical model      \n!        1/ Shallow water\n!        2/ Euler isothermal two-phases               (default 0.)\n!---- for isothermal Euler model\n! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)\n! PREF : isothermal pressure of reference             (default 1.d5)\n! CSSM : Mixture sound velocity (isothermal)          (default 20.)\n! REFA : Isothermal reference density for Air         (default 1.)\n! REFW : Isothermal reference density for Water       (default 1000.)\n! SHAR : isothermal interface sharpening parameter \n!        recommended value 0.01                       (defaut 0.)\n!---- for Shallow water model\n! WLTR : Value to truncate the Water level            (default 0.0001)\n! FRMO : Friction Model\n! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)\n! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation\n!-----------------------------------------------------------------------  \nPHYS 0\nMODE 1.\nFRMO 1.\n\n!-----------------------------------------------------------------------\n! Reading the MAster MEsh\n! Header  MAME \n!-----------------------------------------------------------------------\n! code 1 -&gt; creation of a mesh\n!\n!           argument: 0 type \"1D shock tube\"\n! next line: xmin,xmax\n! next line: nx\n!\n!           argument: 1 type \"2d box\"\n! next line: xmin,xmax,ymin,ymax\n! next line: nx,ny\n!\n!           argument: 2 type \"3d box\"\n! next line: xmin,xmax,ymin,ymax,zmin,zmax\n! next line: nx,ny,nz\n!\n! code 2 -&gt; Reading an ascii file format GMSH V2.2\n! next line: file name\n!-----------------------------------------------------------------------\nMAME    1    0\n0. 1000.  \n1000\n!-----------------------------------------------------------------------\n! Reading Initial conditions\n! Header INIT \n!-----------------------------------------------------------------------\n! code 0 -&gt; Definition for a shock tube\n!           primitive variables on left  \n!           primitive variables on right\n!           argument: 0\n! code 1 -&gt; user function \n!           argument: Function's Number \n! code 2 -&gt; Definition using areas\n!           argument: number of areas\n!           and for each area 2 lignes\n!           - xmin,xmax,ymin,ymax,zmin,zmax\n!           - nvar values ( primitive variables)\n! Shallow water model: water level, x velocity, y velocity\n! Air-water flow model: density, x velocity, y velocity,z velocity, pressure\n!-----------------------------------------------------------------------\nINIT    2   1\n-1. 1001. -100. 100. -100. 100.\n0. 0. 0.\n\n!-----------------------------------------------------------------------\n! Reading bathymetry and friction\n! Header BATH \n!-----------------------------------------------------------------------\n! code 0 -&gt; Default : flat bottom and no friction\n\n! code 1 -&gt; bathymetry and friction defined by user function (defined in sources/PHY/fct_bathy.f90)\n!           argument: Function's Number \n! code 2 -&gt; bathymetry and friction defined linearly per block \n!           argument: 0  reading the file \n!           argument:-1  creation of a file containing the coordinates of the vertices of the block \n! next line: file name\n!-----------------------------------------------------------------------\n\nBATH    1   2\n!-----------------------------------------------------------------------\n! Reading boundary condition\n! (maximum 10 for the moment)\n! Header COND \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; Description in the case  \"shock tube\" : Nothing\n!\n! code 1 -&gt; Description in the case  \"2d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax\n!\n! code 2 -&gt; Description in the case  \"3d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax\n!\n! code 3 -&gt; Definition per zone for  GMSH mesh\n!           argument: number of zones defining a boundary condition\n!           number of zone, kind of boundary condition \n!\n! kind of boundary condition:\n! 0 outlet (we copy)\n! 1 miror\n! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values\n! 3 user's function , then next line function number\n! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),\n!   then next line nvar values\n! 999 Non used\n!-----------------------------------------------------------------------\nCOND    2    0\n4\n1.e20 2. 0.\n2\n0.74832355831838937 1.e20 0. 0.\n999\n999\n999\n999\n\n!-----------------------------------------------------------------------\n! Reading mesh parameters in order to define the mesh according to\n! the Master mesh is the one represented by the domain \"0\".\n! Header MESH \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; default value\n!\n! 1 keyword (4 characters) and a real\n! NBDS : NumBer of DomainS                         (default 001.)\n! MARL : MAximal Refinement Level                  (default 000.)\n! COPA : Mesh COarsening PArameter          0&lt;..&lt;1 (default 0.002)\n! REPA : Mesh REfinement PArameter          0&lt;..&lt;1 (default 0.02)\n! FDRA : Function defining the mesh refinement to be applied \n! NZRA : Number of zones where initial blocks are refined (default 1.)\n!        then for each zone\n!        nrb,xmin,xmax,ymin,ymax,zmin,zmax\n!-----------------------------------------------------------------------\nMESH 0\nNBDS 1\nMARL 0\nNZRA 1\n0 0.  1000. -100. 100. -100. 100.\n\n!-----------------------------------------------------------------------\n! Reading numerical parameter \n! Header  NUME\n!-----------------------------------------------------------------------\n! 1 keyword (4 characters) then a real\n! code 0 -&gt; default value\n!\n! TINT : Time interval to compute                    (default 0.)\n! CCFL : cfl condition                               (default 0.9)\n! TDOR : Time discretization order 1 ou 2         (defaut 1.)\n! SDOR : Space discretization order 1 ou 2        (defaut 1.)\n! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)\n! SAVE : kind of backup, three-digit number          (default 001.)\n!        decade  : shock tube backup 0/no 1/yes \n!        unit    : binary backup     0/no 1/yes \n! NPRO : number of probe  &lt; 10                      (defaut 000.)\n!        then coordinates x,y,z of the probe (1 probe per line)\n!-----------------------------------------------------------------------\nNUME 0\nTINT 50.\nCCFL .9\nSDOR 1.\nTDOR 1.\nSAVE 01.\n\n\n\n\n\nTo run the computation, in order to visualize the evolution towards the steady state, the computation time is split into \\(30\\) steps using the script below:\n\n\n\n\n\n\nTipScript ex3_sv_McDo.sh\n\n\n\n#!/bin/bash\n./cerf_input ex3_sv_McDo.inp\n./cerf2tec\nmv cerfout1d.dat md_0.dat\nfor i in $(seq 1 30 )\ndo\n./cerf\n./cerf2tec\nmv cerfout1d.dat md_$i.dat\ndone\n\n\n\nThe analytical solution to this problem can be obtained from the Swashes website. In the doc directory of CERF, the analytical solution on 10000 points in Tecplot format is provided for reference. The numerical solution can then be compared to the exact solution in the figure below (Figure 3.5).\n\n\n\n\n\n\n\nFigure 3.5: Numerical (black) / exact (red) comparison of the water level \\(h\\) (top) and velocity \\(u\\) (bottom), with a uniform mesh of 1000 cells at first order in space and time",
    "crumbs": [
      "<span class='chapter-number'>3</span>  <span class='chapter-title'>Documented examples</span>"
    ]
  },
  {
    "objectID": "chapitre3.html#examples-with-the-eulerian-two-fluid-model",
    "href": "chapitre3.html#examples-with-the-eulerian-two-fluid-model",
    "title": "3  Documented examples",
    "section": "3.2 Examples with the Eulerian two-fluid model",
    "text": "3.2 Examples with the Eulerian two-fluid model\nIn this chapter, we present the implementation of a few classic test cases for two-fluid flows. The cases covered are as follows:\n\nSolving a Riemann problem\nDam break\nA cylinder falling onto a body of water \n…\n\n\n3.2.1 Riemann problem\nIn order to validate the implementation of the model, we consider the classic “shock tube” test case with the isothermal two-fluid model. We use the test case from Frédéric Golay and Helluy (2007). The domain is defined by \\(x\\in[-2.,2.]\\) and \\(t\\in[0,0.25]\\). The initial conditions are given by: \\[\n\\overrightarrow{w_p}(x,0) =\\left\\{\n\\begin{array}{ll}\n\\overrightarrow{w_p}_L = &lt;\\rho_L,u_L,p_L&gt;^T = &lt;0.0046,\\;0,\\;1.5 \\;10^5&gt;^T &\\mbox{ if } x&lt;0 \\\\\n\\overrightarrow{w_p}_R = &lt;\\rho_R,u_R,p_R&gt;^T = &lt;0.1446,\\;0.,\\;1. \\;10^5&gt;^T &\\mbox{ if } x&gt;0\n\\end{array}\n\\right.\n\\]\nWe consider a uniform, fixed 1D mesh of \\(400\\) cells over a single domain. The data file corresponding to this test case is provided below:\n\n\n\n\n\n\nTip 3.5: Data file ex1_bf_Riemann.inp\n\n\n\n\n\n!-----------------------------------------------------------------------\n! Data file in CERF format\n! Each block of data is characterized in the first line by a header (4 characters)\n! then a code (1 integer) and an argument (1 integer)\n! PHYS MAME INIT BATH COND MESH OBST NUME\n!-----------------------------------------------------------------------\n\n!-----------------------------------------------------------------------\n! Reading physical parameters \n! Header PHYS\n! 1 keyword (4 characters) and a code (1 integer)\n! for now we have the same values for all domains\n!-----------------------------------------------------------------------\n! MODE : physical model      \n!        1/  Shallow water\n!        2/  Euler isothermal two-phases              (default 0.)\n!---- for isothermal Euler model\n! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)\n! PREF : isothermal pressure of reference             (default 1.d5)\n! CSSM : Mixture sound velocity (isothermal)          (default 20.)\n! REFA : Isothermal reference density for Air         (default 1.)\n! REFW : Isothermal reference density for Water       (default 1000.)\n! SHAR : isothermal interface sharpening parameter \n!        recommended value 0.01                       (defaut 0.)\n!---- for Shallow water model\n! WLTR : Value to truncate the Water level            (default 0.0001)\n! FRMO : Friction Model\n! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)\n! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation\n!-----------------------------------------------------------------------  \nPHYS 0\nMODE 2.\n\n\n!-----------------------------------------------------------------------\n! Reading the MAster MEsh\n! Header  MAME \n!-----------------------------------------------------------------------\n! code 1 -&gt; creation of a mesh\n!\n!           argument: 0 type \"1D shock tube\"\n! next line: xmin,xmax\n! next line: nx\n!\n!           argument: 1 type \"2d box\"\n! next line: xmin,xmax,ymin,ymax\n! next line: nx,ny\n!\n!           argument: 2 type \"3d box\"\n! next line: xmin,xmax,ymin,ymax,zmin,zmax\n! next line: nx,ny,nz\n!\n! code 2 -&gt; Reading an ascii file format GMSH V2.2\n! next line: file name\n!-----------------------------------------------------------------------\nMAME    1    0\n-10. 10.  \n400\n!-----------------------------------------------------------------------\n! Reading Initial conditions\n! Header INIT \n!-----------------------------------------------------------------------\n! code 0 -&gt; Definition for a shock tube\n!           primitive variables on left  \n!           primitive variables on right\n!           argument: 0\n! code 1 -&gt; user function \n!           argument: Function's Number \n! code 2 -&gt; Definition using areas\n!           argument: number of areas\n!           and for each area 2 lignes\n!           - xmin,xmax,ymin,ymax,zmin,zmax\n!           - nvar values ( primitive variables)\n! Shallow water model: water level, x velocity, y velocity\n! Air-water flow model: density, x velocity, y velocity,z velocity, pressure\n!-----------------------------------------------------------------------\nINIT    0   0\n1125. 0. 1.5e5\n1. 0. 1.e5\n\n!-----------------------------------------------------------------------\n! Reading boundary condition\n! (maximum 10 for the moment)\n! Header COND \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; Description in the case  \"shock tube\" : Nothing\n!\n! code 1 -&gt; Description in the case  \"2d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax\n!\n! code 2 -&gt; Description in the case  \"3d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax\n!\n! code 3 -&gt; Definition per zone for  GMSH mesh\n!           argument: number of zones defining a boundary condition\n!           number of zone, kind of boundary condition \n!\n! kind of boundary condition:\n! 0 outlet (we copy)\n! 1 miror\n! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values\n! 3 user's function , then next line function number\n! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),\n!   then next line nvar values\n! 999 Non used\n!-----------------------------------------------------------------------\nCOND    0    0\n\n\n!-----------------------------------------------------------------------\n! Reading mesh parameters in order to define the mesh according to\n! the Master mesh is the one represented by the domain \"0\".\n! Header MESH \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; default value\n!\n! 1 keyword (4 characters) and a real\n! NBDS : NumBer of DomainS                         (default 001.)\n! MARL : MAximal Refinement Level                  (default 000.)\n! COPA : Mesh COarsening PArameter          0&lt;..&lt;1 (default 0.002)\n! REPA : Mesh REfinement PArameter          0&lt;..&lt;1 (default 0.02)\n! FDRA : Function defining the mesh refinement to be applied \n! NZRA : Number of zones where initial blocks are refined (default 1.)\n!        then for each zone\n!        nrb,xmin,xmax,ymin,ymax,zmin,zmax\n!-----------------------------------------------------------------------\nMESH 0\nNBDS 1\nMARL 0\nNZRA 1\n0 -10.  10. -100. 100. -100. 100.\n\n!-----------------------------------------------------------------------\n! Reading numerical parameter \n! Header  NUME\n!-----------------------------------------------------------------------\n! 1 keyword (4 characters) then a real\n! code 0 -&gt; default value\n!\n! TINT : Time interval to compute                    (default 0.)\n! CCFL : cfl condition                               (default 0.9)\n! TDOR : Time discretization order 1 ou 2         (defaut 1.)\n! SDOR : Space discretization order 1 ou 2        (defaut 1.)\n! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)\n! SAVE : kind of backup, three-digit number          (default 001.)\n!        decade  : shock tube backup 0/no 1/yes \n!        unit    : binary backup     0/no 1/yes \n! NPRO : number of probe  &lt; 10                      (defaut 000.)\n!        then coordinates x,y,z of the probe (1 probe per line)\n!-----------------------------------------------------------------------\nNUME 0\nTINT 0.25\nCCFL .9\nSDOR 2.\nTDOR 2.\nSAVE 11.\n\n\n\n\n\nTo run the computation, since there is only one domain to process over a single time interval:\n./cerf_input ex1_bf_Riemann.inp\n./cerf\nSince the data file ex1_bf_Riemann.inp requests output in 1D format (parameter SAVE=11, section NUME), 1D output files have been generated. Files for the numerical solution in density, velocity, and pressure ‘n_r_001’, ‘n_u_001’, and ‘n_p_001’, and files for the exact solution in density, velocity, and pressure: ‘e_r_001’, ‘e_u_001’, and ‘e_p_001’. These files can be visualized, for example, with gnuplot:\ngnuplot\ngnuplot&gt; set term x11\ngnuplot&gt; plot 'e_r_001' w l, 'n_r_001' w l\ngnuplot&gt; plot 'e_u_001' w l, 'n_u_001' w l\ngnuplot&gt; plot 'e_p_001' w l, 'n_p_001' w l\ngnuplot&gt; quit\nThe numerical solution can thus be compared to the exact solution.\n\n\n\n\nNumerical/exact comparison of density at \\(t=0.25s\\), uniform mesh of 400 cells\n\n\n\n\n\n\n\nNumerical/exact comparison of velocity at \\(t=0.25s\\), uniform mesh of 400 cells\n\n\n\n\n\n3.2.2 Dam break\nWe consider here a 2D dam-break test case. This is a widely used and well-documented case. From an experimental standpoint, one can rely for example on the experiments carried out by Martin and Moyce (1952). We consider a computational domain \\(x\\in[0.,4.]\\) and \\(y\\in[0.,3.]\\). The domain boundaries are characterized by “mirror”-type or zero-Neumann boundary conditions. At the initial instant, the domain consists of air, except for a region defined by \\(x\\in[0.,2.]\\) and \\(y\\in[0.,3.]\\) where water is present. The master mesh consists of \\(40 \\times 30\\) blocks, which are refined to level 4 around the air/water interface. To improve the quality of the interface, the computation is penalized using an interface stiffening coefficient of \\(0.02\\).\nThe data file corresponding to this test case is provided below:\n\n\n\n\n\n\nTip 3.6: Data file ex2_bf_Dam.inp\n\n\n\n\n\n!-----------------------------------------------------------------------\n! Data file in CERF format\n! Each block of data is characterized in the first line by a header (4 characters)\n! then a code (1 integer) and an argument (1 integer)\n! PHYS MAME INIT BATH COND MESH OBST NUME\n!-----------------------------------------------------------------------\n\n!-----------------------------------------------------------------------\n! Reading physical parameters \n! Header PHYS\n! 1 keyword (4 characters) and a code (1 integer)\n! for now we have the same values for all domains\n!-----------------------------------------------------------------------\n! MODE : physical model      \n!        1/ Shallow water\n!        2/ Euler isothermal two-phases               (default 0.)\n!---- for isothermal Euler model\n! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)\n! PREF : isothermal pressure of reference             (default 1.d5)\n! CSSM : Mixture sound velocity (isothermal)          (default 20.)\n! REFA : Isothermal reference density for Air         (default 1.)\n! REFW : Isothermal reference density for Water       (default 1000.)\n! SHAR : isothermal interface sharpening parameter \n!        recommended value 0.01                       (defaut 0.)\n!---- for Shallow water model\n! WLTR : Value to truncate the Water level            (default 0.0001)\n! FRMO : Friction Model\n! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)\n! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation\n!-----------------------------------------------------------------------  \nPHYS 0\nMODE 3.\nGPES 2\nSHAR 0.02\n!-----------------------------------------------------------------------\n! Reading the MAster MEsh\n! Header  MAME \n!-----------------------------------------------------------------------\n! code 1 -&gt; creation of a mesh\n!\n!           argument: 0 type \"1D shock tube\"\n! next line: xmin,xmax\n! next line: nx\n!\n!           argument: 1 type \"2d box\"\n! next line: xmin,xmax,ymin,ymax\n! next line: nx,ny\n!\n!           argument: 2 type \"3d box\"\n! next line: xmin,xmax,ymin,ymax,zmin,zmax\n! next line: nx,ny,nz\n!\n! code 2 -&gt; Reading an ascii file format GMSH V2.2\n! next line: file name\n!-----------------------------------------------------------------------\nMAME    1    1\n0. 4.  0. 3.\n40 30\n!-----------------------------------------------------------------------\n! Reading Initial conditions\n! Header INIT \n!-----------------------------------------------------------------------\n! code 0 -&gt; Definition for a shock tube\n!           primitive variables on left  \n!           primitive variables on right\n!           argument: 0\n! code 1 -&gt; user function \n!           argument: Function's Number \n! code 2 -&gt; Definition using areas\n!           argument: number of areas\n!           and for each area 2 lignes\n!           - xmin,xmax,ymin,ymax,zmin,zmax\n!           - nvar values ( primitive variables)\n! Shallow water model: water level, x velocity, y velocity\n! Air-water flow model: density, x velocity, y velocity,z velocity, pressure\n!-----------------------------------------------------------------------\nINIT    2    2\n-1. 4.1 -1. 3.1 -1. 1.\n1. 0. 0. 0.  1.e5 0.\n-1. 1. -1. 2. -1. 1.\n1000. 0. 0. 0.  1.e5 0.\n\n!-----------------------------------------------------------------------\n! Reading boundary condition\n! (maximum 10 for the moment)\n! Header COND \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; Description in the case  \"shock tube\" : Nothing\n!\n! code 1 -&gt; Description in the case  \"2d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax\n!\n! code 2 -&gt; Description in the case  \"3d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax\n!\n! code 3 -&gt; Definition per zone for  GMSH mesh\n!           argument: number of zones defining a boundary condition\n!           number of zone, kind of boundary condition \n!\n! kind of boundary condition:\n! 0 outlet (we copy)\n! 1 miror\n! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values\n! 3 user's function , then next line function number\n! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),\n!   then next line nvar values\n! 999 Non used\n!-----------------------------------------------------------------------\nCOND    1    0\n1\n1\n1\n1\n\n!-----------------------------------------------------------------------\n! Reading mesh parameters in order to define the mesh according to\n! the Master mesh is the one represented by the domain \"0\".\n! Header MESH \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; default value\n!\n! 1 keyword (4 characters) and a real\n! NBDS : NumBer of DomainS                         (default 001.)\n! MARL : MAximal Refinement Level                  (default 000.)\n! COPA : Mesh COarsening PArameter          0&lt;..&lt;1 (default 0.002)\n! REPA : Mesh REfinement PArameter          0&lt;..&lt;1 (default 0.02)\n! FDRA : Function defining the mesh refinement to be applied \n! NZRA : Number of zones where initial blocks are refined (default 1.)\n!        then for each zone\n!        nrb,xmin,xmax,ymin,ymax,zmin,zmax\n!-----------------------------------------------------------------------\nMESH 0\nNBDS 8\nMARL 4\nCOPA 0.9\nREPA 3\nNZRA 3\n0 -1. 4.1 -1. 3.1 -1. 1.\n4 -1. 1.5 -1. 2.25 -1. 1.\n0 -1. .75 -1. 1.75 -1. 1.\n\n!-----------------------------------------------------------------------\n! Reading numerical parameter \n! Header  NUME\n!-----------------------------------------------------------------------\n! 1 keyword (4 characters) then a real\n! code 0 -&gt; default value\n!\n! TINT : Time interval to compute                    (default 0.)\n! CCFL : cfl condition                               (default 0.9)\n! TDOR : Time discretization order 1 ou 2         (defaut 1.)\n! SDOR : Space discretization order 1 ou 2        (defaut 1.)\n! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)\n! SAVE : kind of backup, three-digit number          (default 001.)\n!        decade  : shock tube backup 0/no 1/yes \n!        unit    : binary backup     0/no 1/yes \n! NPRO : number of probe  &lt; 10                      (defaut 000.)\n!        then coordinates x,y,z of the probe (1 probe per line)\n!-----------------------------------------------------------------------\nNUME 0\nTINT 0.01\nCCFL .9\nSDOR 2.\nTDOR 2.\nSAVE 1.\n\n\n\n\n\nThe expected velocities are on the order of a few meters per second. Therefore, since the blocks have a size of \\(0.1 m\\), the time interval between each remeshing can be estimated to be \\(0.01 s\\). To speed up the computation, the mesh is decomposed into 8 domains (or equivalently 8 MPI processes). To run a computation up to \\(T=0.8 s\\), the following script can be used:\n#!/bin/bash\n./cerf_input ex2_bf_Dam.inp\n./cerf2tec\nmv cerfout.dat d_0.dat\nmv cerfamr.dat d_mame_0.dat\nfor i in $(seq 1 80 )\ndo\nmpirun -np 8 ./cerf\n./cerf2tec\nmv cerfout.dat d_$i.dat\n./cerf_amr\nmv cerfamr.dat d_mame_$i.dat\ndone\nIn the animation below Figure 3.6, the density of the mixture is shown (blue for water and red for air). The interface can be seen to be “sharp”, i.e., spread over 3-4 cells.\n\n\n\n\n\n\n\nFigure 3.6: Dam break, density\n\n\n\n\n\n\n3.2.3 A cylinder falling onto a body of water\nWe consider here a 2D test case of a solid object falling onto a water surface. This fluid-structure interaction problem is treated here using CERF’s “rigid-fluid” tool (Frederic Golay (2018)). This is a widely used and well-documented case (see for example H. Sun and Faltinsen (2006), P. Sun, Ming, and Zhang (2015)). From an experimental standpoint, one can rely for example on the experiments carried out by Greenhow and Lin (2015). We consider a computational domain \\(x\\in[0.,2.]\\) and \\(y\\in[0.,2.]\\). The domain boundaries are characterized by “mirror”-type or zero-Neumann boundary conditions. At the initial instant, the domain consists of air, except for a region defined by \\(x\\in[0.,2.]\\) and \\(y\\in[0.,0.65]\\) where water is present. The rigid solid is defined by a mesh in .OBJ format. The cells inside this envelope behave as a free rigid solid, subject to gravity, with a density of \\(500 kg/m^3\\). The master mesh consists of \\(40 \\times 40\\) blocks, which are refined to level \\(4\\) around the air/water interface and around the cylinder, whose center is positioned \\(0.5m\\) from the surface of the water, as illustrated in figure Figure 3.7. To improve the quality of the interface, the computation is penalized using an interface stiffening coefficient of \\(0.01\\).\n\n\n\n\n\n\n\nFigure 3.7: Description of the test case: a cylinder falling onto a water surface\n\n\n\n\nThe data file corresponding to this test case is provided below:\n\n\n\n\n\n\nTip 3.7: Data file ex3_bf_cyl.inp\n\n\n\n\n\n!-----------------------------------------------------------------------\n! Data file in CERF format\n! Each block of data is characterized in the first line by a header (4 characters)\n! then a code (1 integer) and an argument (1 integer)\n! PHYS MAME INIT BATH COND MESH OBST NUME\n!-----------------------------------------------------------------------\n\n!-----------------------------------------------------------------------\n! Reading physical parameters \n! Header PHYS\n! 1 keyword (4 characters) and a code (1 integer)\n! for now we have the same values for all domains\n!-----------------------------------------------------------------------\n! MODE : physical model      \n!        1/ Shallow water\n!        2/ Euler isothermal two-phases               (default 0.)\n!---- for isothermal Euler model\n! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)\n! PREF : isothermal pressure of reference             (default 1.d5)\n! CSSM : Mixture sound velocity (isothermal)          (default 20.)\n! REFA : Isothermal reference density for Air         (default 1.)\n! REFW : Isothermal reference density for Water       (default 1000.)\n! SHAR : isothermal interface sharpening parameter \n!        recommended value 0.01                       (defaut 0.)\n!---- for Shallow water model\n! WLTR : Value to truncate the Water level            (default 0.0001)\n! FRMO : Friction Model\n! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)\n! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation\n!-----------------------------------------------------------------------  \nPHYS 0\nMODE 2.\nGPES 2\nSHAR 0.01\n!-----------------------------------------------------------------------\n! Reading the MAster MEsh\n! Header  MAME \n!-----------------------------------------------------------------------\n! code 1 -&gt; creation of a mesh\n!\n!           argument: 0 type \"1D shock tube\"\n! next line: xmin,xmax\n! next line: nx\n!\n!           argument: 1 type \"2d box\"\n! next line: xmin,xmax,ymin,ymax\n! next line: nx,ny\n!\n!           argument: 2 type \"3d box\"\n! next line: xmin,xmax,ymin,ymax,zmin,zmax\n! next line: nx,ny,nz\n!\n! code 2 -&gt; Reading an ascii file format GMSH V2.2\n! next line: file name\n!-----------------------------------------------------------------------\nMAME    1    1\n0. 2.  0. 2.\n40 40\n!-----------------------------------------------------------------------\n! Reading Initial conditions\n! Header INIT \n!-----------------------------------------------------------------------\n! code 0 -&gt; Definition for a shock tube\n!           primitive variables on left  \n!           primitive variables on right\n!           argument: 0\n! code 1 -&gt; user function \n!           argument: Function's Number \n! code 2 -&gt; Definition using areas\n!           argument: number of areas\n!           and for each area 2 lignes\n!           - xmin,xmax,ymin,ymax,zmin,zmax\n!           - nvar values ( primitive variables)\n! Shallow water model: water level, x velocity, y velocity\n! Air-water flow model: density, x velocity, y velocity,z velocity, pressure\n!-----------------------------------------------------------------------\nINIT    2    2\n-1. 2.1 -1. 2.1 -1. 1.\n1. 0. 0. 0.  1.e5 0.\n-1. 2. -1. 0.65 -1. 1.\n1000. 0. 0. 0.  1.e5 0.\n\n!-----------------------------------------------------------------------\n! Reading boundary condition\n! (maximum 10 for the moment)\n! Header COND \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; Description in the case  \"shock tube\" : Nothing\n!\n! code 1 -&gt; Description in the case  \"2d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax\n!\n! code 2 -&gt; Description in the case  \"3d box\" (at least 1 line per type) :\n!           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax\n!\n! code 3 -&gt; Definition per zone for  GMSH mesh\n!           argument: number of zones defining a boundary condition\n!           number of zone, kind of boundary condition \n!\n! kind of boundary condition:\n! 0 outlet (we copy)\n! 1 miror\n! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values\n! 3 user's function , then next line function number\n! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),\n!   then next line nvar values\n! 999 Non used\n!-----------------------------------------------------------------------\nCOND    1    0\n1\n1\n1\n1\n\n!-----------------------------------------------------------------------\n! Reading mesh parameters in order to define the mesh according to\n! the Master mesh is the one represented by the domain \"0\".\n! Header MESH \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; default value\n!\n! 1 keyword (4 characters) and a real\n! NBDS : NumBer of DomainS                         (default 001.)\n! MARL : MAximal Refinement Level                  (default 000.)\n! COPA : Mesh COarsening PArameter          0&lt;..&lt;1 (default 0.002)\n! REPA : Mesh REfinement PArameter          0&lt;..&lt;1 (default 0.02)\n! FDRA : Function defining the mesh refinement to be applied \n! NZRA : Number of zones where initial blocks are refined (default 1.)\n!        then for each zone\n!        nrb,xmin,xmax,ymin,ymax,zmin,zmax\n!-----------------------------------------------------------------------\nMESH 0\nNBDS 10\nMARL 4\nCOPA 1\nREPA 2\nNZRA 4\n0 -1. 2.1 0. 2.1 -1. 1\n4 -1. 2.1 0.6 0.7 -1. 1\n0 -1. 2.1 0. 0.6 -1. 1\n4 0.9 1.1  1.0 1.3  -1. 1\n\n!-----------------------------------------------------------------------\n! Reading rigid obstacle parameters\n! Header OBST \n!-----------------------------------------------------------------------\n!\n! code 0 -&gt; default value\n!           argument: number of obstacles\n! 1 line per obstacle: file name, density, kind of motion\n! (  -1: fixed, 0: free, &gt;0 : function number imposing the movement\n!\n!-----------------------------------------------------------------------\nOBST    0    1\nex3_bf_cyl.obj 500. 0\n\n!-----------------------------------------------------------------------\n! Reading numerical parameter \n! Header  NUME\n!-----------------------------------------------------------------------\n! 1 keyword (4 characters) then a real\n! code 0 -&gt; default value\n!\n! TINT : Time interval to compute                    (default 0.)\n! CCFL : cfl condition                               (default 0.9)\n! TDOR : Time discretization order 1 ou 2         (defaut 1.)\n! SDOR : Space discretization order 1 ou 2        (defaut 1.)\n! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)\n! SAVE : kind of backup, three-digit number          (default 001.)\n!        decade  : shock tube backup 0/no 1/yes \n!        unit    : binary backup     0/no 1/yes \n! NPRO : number of probe  &lt; 10                      (defaut 000.)\n!        then coordinates x,y,z of the probe (1 probe per line)\n!-----------------------------------------------------------------------\nNUME 0\nTINT 0.01\nCCFL .9\nSDOR 2.\nTDOR 2.\nTIME 1\nSAVE 1.\n\n\n\n\n\nThe expected velocities are on the order of a few meters per second. Therefore, since the blocks have a size of \\(0.05 m\\), the time interval between each remeshing can be estimated to be \\(0.01 s\\). To speed up the computation, the mesh is decomposed into 10 domains (or equivalently 10 MPI processes). To run a computation up to \\(T=2 s\\), the following script can be used:\n#!/bin/bash\n./cerf_input ex3_bf_cyl.inp\n./cerf2tec\nmv cerfout.dat c_000.dat\nmv cerfamr.dat c_mame_000.dat\ncp ex3_bf_cyl.obj cs_000.obj\n\nfor i in $(seq 1 9 )\ndo\nmpirun -np 10 ./cerf\n./cerf2tec\nmv cerfout.dat c_00$i.dat\n./cerf_amr\nmv cerfamr.dat c_mame_00$i.dat\nmv solide_001.obj cs_00$i.obj\ndone\nfor i in $(seq 10 99 )\ndo\nmpirun -np 10 ./cerf\n./cerf2tec\nmv cerfout.dat c_0$i.dat\n./cerf_amr\nmv cerfamr.dat c_mame_0$i.dat\nmv solide_001.obj cs_0$i.obj\ndone\nfor i in $(seq 100 200 )\ndo\nmpirun -np 10 ./cerf\n./cerf2tec\nmv cerfout.dat c_$i.dat\n./cerf_amr\nmv cerfamr.dat c_mame_$i.dat\nmv solide_001.obj cs_$i.obj\ndone\nIn the images below, the density of the mixture is shown (blue for water and red for air). The interface can be seen to be “sharp”, i.e., spread over 3-4 cells, and the dynamic mesh can be seen to “follow” the regions of interest closely. The first figure Figure 3.8, at \\(t=0.3 s\\), corresponds to the instant when the cylinder touches the water surface at a velocity of about \\(3m/s\\), which is consistent with theory (\\(\\sqrt{2gh}\\)) and with experiment. The second figure Figure 3.9, at \\(t=0.45 s\\), illustrates the generation of the wave created as the cylinder enters the water.\n\n\n\n\n\n\n\nFigure 3.8: Cylinder fall, density and mesh at \\(t=0.3s\\)\n\n\n\n\n\n\n\n\n\n\n\nFigure 3.9: Cylinder fall, density and mesh at \\(t=0.45s\\)\n\n\n\n\nIn the animation below Figure 3.10, the density of the mixture is shown (blue for water and red for air). The cylinder is shown in white.\n\n\n\n\n\n\n\nFigure 3.10: Cylinder fall: density\n\n\n\n\nIn the figure below Figure 3.11, the vertical trajectory of the cylinder over time is compared with the experimental results of Greenhow and Lin (2015). Good agreement is observed between the two approaches.\n\n\n\n\n\n\n\nFigure 3.11: Vertical trajectory of the cylinder over time: comparison between the numerical simulation (solid line) and the experiment (dots)\n\n\n\n\n\n\n\n\n\n\nDelestre, Olivier, Carine Lucas, Pierre-Antoine Ksinant, Frédéric Darboux, Christian Laguerre, T.-N.-Tuoi Vo, François James, and Stéphane Cordier. 2013. “SWASHES: A Compilation of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies.” International Journal for Numerical Methods in Fluids 72 (3): 269–300. https://doi.org/https://doi.org/10.1002/fld.3741.\n\n\nFaccanoni, Gloria, and Anne Mangeney. 2013. “Exact Solution for Granular Flows.” International Journal for Numerical and Analytical Methods in Geomechanics 37 (10): 1408–33. https://doi.org/https://doi.org/10.1002/nag.2124.\n\n\nGolay, Frederic. 2018. “Monolithic Fluid Structure Interaction Model, Application to Water Entry Problem.” In Topical Problems of Fluid Mechanics 2018. Prague, Czech Republic: Institute of Thermomechanics AS CR v.v.i. https://doi.org/10.14311/TPFM.2018.016.\n\n\nGolay, Frédéric, and Philippe Helluy. 2007. “Numerical Schemes for Low Mach Wave Breaking.” International Journal of Computational Fluid Dynamics 21 (2): 69–86. https://doi.org/10.1080/10618560701343382.\n\n\nGreenhow, Martin, and and Lin. 2015. “Greenhow m. And Lin w-m.\"nonlinear Free Surface Effects: Experiments and Theory\" MIT Internal Report 83-19 (1983),” June. https://doi.org/10.13140/RG.2.1.4129.2964.\n\n\nMartin, J. C., and W. J. Moyce. 1952. “Part IV. An Experimental Study of the Collapse of Liquid Columns on a Rigid Horizontal Plane.” Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 244 (882): 312–24. http://www.jstor.org/stable/91519.\n\n\nSun, Hui, and Odd M. Faltinsen. 2006. “Water Impact of Horizontal Circular Cylinders and Cylindrical Shells.” Applied Ocean Research 28 (5): 299–311. https://doi.org/https://doi.org/10.1016/j.apor.2007.02.002.\n\n\nSun, Pengnan, Furen Ming, and Aman Zhang. 2015. “Numerical Simulation of Interactions Between Free Surface and Rigid Body Using a Robust SPH Method.” Ocean Engineering 98: 32–49. https://doi.org/https://doi.org/10.1016/j.oceaneng.2015.01.019.",
    "crumbs": [
      "<span class='chapter-number'>3</span>  <span class='chapter-title'>Documented examples</span>"
    ]
  },
  {
    "objectID": "references.html",
    "href": "references.html",
    "title": "Références",
    "section": "",
    "text": "Altazin, Thomas. 2017. “Un Modèle d’interaction Fluide-Structure\nEn Régime Compressible Faible Mach.” PhD thesis. http://www.theses.fr/2017TOUL0021/document.\n\n\nAudusse, Emmanuel, François Bouchut, Marie-Odile Bristeau, Rupert Klein,\nand Benoı⁁t Perthame. 2004. “A Fast and Stable Well-Balanced\nScheme with Hydrostatic Reconstruction for Shallow Water Flows.”\nSIAM Journal on Scientific Computing 25 (6): 2050–65. https://doi.org/10.1137/S1064827503431090.\n\n\nChorin, Alexandre Joel. 1967. “A Numerical Method for Solving\nIncompressible Viscous Flow Problems.” Journal of\nComputational Physics 2 (1): 12–26. https://doi.org/https://doi.org/10.1016/0021-9991(67)90037-X.\n\n\nCoquerelle, M., and G.-H. Cottet. 2008. “A Vortex Level Set Method\nfor the Two-Way Coupling of an Incompressible Fluid with Colliding Rigid\nBodies.” Journal of Computational Physics 227 (21):\n9121–37. https://doi.org/https://doi.org/10.1016/j.jcp.2008.03.041.\n\n\nCroisille, Jean-Pierre. 1990. “Contribution à l’étude Théorique Et\nà l’approximation Par Éléments Finis Du Système Hyperbolique de La\nDynamique Des Gaz Multidimensionnelle Et Multiespèces.” PhD\nthesis.\nhttp://www.theses.fr/1990PA066095.\n\n\nDelestre, Olivier, Carine Lucas, Pierre-Antoine Ksinant, Frédéric\nDarboux, Christian Laguerre, T.-N.-Tuoi Vo, François James, and Stéphane\nCordier. 2013. “SWASHES: A Compilation of Shallow Water Analytic\nSolutions for Hydraulic and Environmental Studies.”\nInternational Journal for Numerical Methods in Fluids 72 (3):\n269–300. https://doi.org/https://doi.org/10.1002/fld.3741.\n\n\nErsoy, Mehmet, Frédéric Golay, and Lyudmyla Yushchenko. 2013. Open\nMathematics 11 (8): 1392–1415. https://doi.org/doi:10.2478/s11533-013-0252-6.\n\n\nFaccanoni, Gloria, and Anne Mangeney. 2013. “Exact Solution for\nGranular Flows.” International Journal for Numerical and\nAnalytical Methods in Geomechanics 37 (10): 1408–33.\nhttps://doi.org/https://doi.org/10.1002/nag.2124.\n\n\nGolay, Frederic. 2018. “Monolithic Fluid\nStructure Interaction Model, Application to Water Entry\nProblem.” In Topical Problems of\nFluid Mechanics 2018. Prague, Czech Republic:\nInstitute of Thermomechanics AS CR v.v.i. https://doi.org/10.14311/TPFM.2018.016.\n\n\nGolay, Frédéric. 2009. “Numerical Entropy Production and Error\nIndicator for Compressible Flows.” Comptes Rendus.\nMécanique 337 (4): 233–37. https://doi.org/10.1016/j.crme.2009.04.004.\n\n\nGolay, Frédéric, and Philippe Helluy. 2007. “Numerical Schemes for\nLow Mach Wave Breaking.” International Journal of\nComputational Fluid Dynamics 21 (2): 69–86. https://doi.org/10.1080/10618560701343382.\n\n\nGreenhow, Martin, and and Lin. 2015. “Greenhow m. And Lin\nw-m.\"nonlinear Free Surface Effects: Experiments and Theory\" MIT\nInternal Report 83-19 (1983),” June. https://doi.org/10.13140/RG.2.1.4129.2964.\n\n\nHelluy, Philippe, Golay, Frédéric, Caltagirone, Jean-Paul, Lubin,\nPierre, Vincent, Stéphane, Drevard, Deborah, Marcer, Richard, et al.\n2005. “Numerical Simulations of Wave Breaking.” ESAIM:\nM2AN 39 (3): 591–607. https://doi.org/10.1051/m2an:2005024.\n\n\nLax, Peter. 1971. “Shock Waves and Entropy.” In\nContributions to Nonlinear Functional Analysis, edited by\nEduardo H. Zarantonello, 603–34. Academic Press. https://doi.org/https://doi.org/10.1016/B978-0-12-775850-3.50018-2.\n\n\nMarche, Fabien. 2007. “Derivation of a New Two-Dimensional Viscous\nShallow Water Model with Varying Topography, Bottom Friction and\nCapillary Effects.” European Journal of Mechanics -\nB/Fluids 26 (1): 49–63. https://doi.org/https://doi.org/10.1016/j.euromechflu.2006.04.007.\n\n\nMartin, J. C., and W. J. Moyce. 1952. “Part IV. An Experimental\nStudy of the Collapse of Liquid Columns on a Rigid Horizontal\nPlane.” Philosophical Transactions of the Royal Society of\nLondon. Series A, Mathematical and Physical Sciences 244 (882):\n312–24. http://www.jstor.org/stable/91519.\n\n\nPons, Kevin. 2018. “Dynamique Des Écoulements à Surface Libre Dans\nLes Plages Sableuses.” Thèse de doctorat, IMATH, Université de\nToulon: Université de Toulon.\n\n\nPoussel, Camille. 2024. “Dynamique Des Écoulements à Surface Libre\nDans Les Plages Sableuses.” Thèse de doctorat, IMATH, Université\nde Toulon: Université de Toulon.\n\n\nPuppo, Gabriella. 2004. “Numerical Entropy Production for Central\nSchemes.” SIAM Journal on Scientific Computing 25 (4):\n1382–1415. https://doi.org/10.1137/S1064827502386712.\n\n\nSaint-Venant, Barré de. 1871. “AJC (1871b).\" Théorie\nEt Equations générales Du Mouvement Non\nPermanent Des Eaux, Avec Application Aux Crues Des Rivières\nEt à l’introduction Des Marées\nDans Leur Lit (2ème Note).\".” Comptes Rendus Des\nséances de l’Académie Des Sciences,\n237–40.\n\n\nSun, Hui, and Odd M. Faltinsen. 2006. “Water Impact of Horizontal\nCircular Cylinders and Cylindrical Shells.” Applied Ocean\nResearch 28 (5): 299–311. https://doi.org/https://doi.org/10.1016/j.apor.2007.02.002.\n\n\nSun, Pengnan, Furen Ming, and Aman Zhang. 2015. “Numerical\nSimulation of Interactions Between Free Surface and Rigid Body Using a\nRobust SPH Method.” Ocean Engineering 98: 32–49.\nhttps://doi.org/https://doi.org/10.1016/j.oceaneng.2015.01.019.\n\n\nThomas Altazin, Frédéric Golay, Mehmet Ersoy, and Lyudmyla Yushchenko.\n2016. “Numerical Investigation of BB-AMR Scheme Using Entropy\nProduction as Refinement Criterion.” International Journal of\nComputational Fluid Dynamics 30 (3): 256–71. https://doi.org/10.1080/10618562.2016.1194977.\n\n\nToro, E. F. 1997. Riemann Solvers and Numerical Methods for Fluid\nDynamics: A Practical Introduction. Springer. https://books.google.fr/books?id=6QFAAQAAIAAJ.\n\n\nVarra, Giada, Veronica Pepe, Renata Della Morte, and Luca Cozzolino.\n2024. “A Novel Efficient and Robust Treatment of the Friction\nSource Term in 2D Shallow Water Inundation Models.” Journal\nof Hydrology 634: 131045. https://doi.org/10.1016/j.jhydrol.2024.131045.",
    "crumbs": [
      "Références"
    ]
  }
]