Some documentation on EB and Poiseuille V&V (#146)

* Add Poiseuille V&V.

* Add some EB documentation.

* Minor fixes.

* Add plot python for Poiseuille (lifted from PeleC).

* Add an input file for the Poiseuille test.

* Update AMReX-Hydro & PP

Docs/source/LMeXControls.rst

@@ -184,6 +184,52 @@ Chemistry integrator
 
 Note that the last four parameters belong to the Reactor class of PelePhysics but are specified here for completeness. In particular, CVODE is the adequate choice of integrator to tackle PeleLMeX large time step sizes. Several linear solvers are available depending on whether or not GPU are employed: on CPU, `dense_direct` is a finite-difference direct solver, `denseAJ_direct` is an analytical-jacobian direct solver (preferred choice), `sparse_direct` is an analytical-jacobian sparse direct solver based on the KLU library and `GMRES` is a matrix-free iterative solver; on GPU `GMRES` is a matrix-free iterative solver (available on all the platforms), `sparse_direct` is a batched block-sparse direct solve based on NVIDIA's cuSparse (only with CUDA), `magma_direct` is a batched block-dense direct solve based on the MAGMA library (available with CUDA and HIP.
 
+Embedded Geometry
+-----------------
+
+`PeleLMeX` geomtry relies on AMReX implementation of the EB method. Simple geometrical objects 
+can thus be constructed using `AMReX internal parser<https://amrex-codes.github.io/amrex/docs_html/EB.html>`_.
+For instance, setting up a sphere of radius 5 mm can be achieved:
+
+::
+
+    eb2.geom_type = sphere
+    eb2.sphere_radius = 0.005
+    eb2.sphere_center = 0.0 0.0 0.0
+    eb2.sphere_has_fluid_inside = 0
+    eb2.small_volfrac = 1.0e-4
+    eb2.maxiter = 200
+
+The `eb2.small_volfrac` controls volume fraction that are deemed too small and eliminated from the EB representation. 
+This operation is done iteratively and the maximum number of iteration is prescribed by `eb2.maxiter`.
+For most applications, a single AMReX object is insufficient to represent the geometry. AMReX enable to combine 
+objects using constructive solid geometry (CSG) in order to create complex geometry. It is up to the user to define
+the combination of basic elements leading to its desired geometry. To switch to a user-defined EB definition, one
+must set:
+
+::
+
+    eb2.geom_type = UserDefined
+
+and then implement the actual geometry definition in a `EBUserDefined.H` file located in the run folder (and added
+to the GNUmakefile using `CEXE_headers += EBUserDefined.H`). An example of such implementation is available in the 
+``Exec/Case/ChallengeProblem`` folder. Example of more generic EB problems are also found in the ``Exec/RegTest/EB_*``
+folders.
+
+In addition to the input keys presented above, a set of `PeleLMeX`-specific keys are available in order to control refinement at the EB:
+
+::
+
+    peleLM.refine_EB_type = Static
+    peleLM.refine_EB_max_level = 1
+    peleLM.refine_EB_buffer = 2.0
+
+By default, the EB is refined to the `amr.max_level`, which can lead to undesirably high number of cells
+close to the EB when the physics of interest might be elsewhere. The above lines enable to limit the
+EB-level to level 1 (must be below `amr.max_level`) and a derefinement strategy is adopted to ensure
+that fine-grid patches do not cross the EB boundary. The last parameter set a safety margin to increase
+how far the derefinement is applied in order to account for grid-patches diagonals and proper nesting contraint.
+Note that the parameter do not ensure coarse-fine/EB crossings are avoided and the code will fail when this happens.
 
 Linear solvers
 --------------

Docs/source/Model.rst

@@ -14,7 +14,7 @@ setup and control of `PeleLMeX` in later sections. `PeleLMeX` is a non-subcyclin
 Overview of `PeleLMeX`
 ----------------------
 
-`PeleLMeX` evolves chemically reacting low Mach number flows with block-structured adaptive mesh refinement (AMR). The code depends upon the `AMReX <https://github.com/AMReX-Codes/amrex>` library to provide the underlying data structures, and tools to manage and operate on them across massively parallel computing architectures. `PeleLMeX` also utilizes the source code and algorithmic infrastructure of `AMReX-Hydro <https://github.com/AMReX-Codes/AMReX-Hydro>`. `PeleLMeX` borrows heavily from `PeleLM <https://github.com/AMReX-Combustion/PeleLM>`.  The core algorithms in `PeleLM` are described in the following papers:
+`PeleLMeX` evolves chemically reacting low Mach number flows with block-structured adaptive mesh refinement (AMR). The code depends upon the `AMReX <https://github.com/AMReX-Codes/amrex>`_ library to provide the underlying data structures, and tools to manage and operate on them across massively parallel computing architectures. `PeleLMeX` also utilizes the source code and algorithmic infrastructure of `AMReX-Hydro <https://github.com/AMReX-Codes/AMReX-Hydro>`_. `PeleLMeX` borrows heavily from `PeleLM <https://github.com/AMReX-Combustion/PeleLM>`_. The core algorithms in `PeleLM` are described in the following papers:
 
 * *A conservative, thermodynamically consistent numerical approach for low Mach number combustion. I. Single-level integration*, A. Nonaka, J. B. Bell, and M. S. Day, *Combust. Theor. Model.*, **22** (1) 156-184 (2018)
 
@@ -27,7 +27,7 @@ Overview of `PeleLMeX`
 * *A Conservative Adaptive Projection Method for the Variable Density Incompressible Navier-Stokes Equations,* A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, and M. L. Welcome, *J. Comp. Phys.*, **142** 1-46 (1998)
 
 The low Mach number flow equations
-----------------------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 `PeleLMeX` solves the reacting Navier-Stokes flow equations in the *low Mach number* regime, where the characteristic fluid velocity is small compared to the sound speed, and the effect of acoustic wave propagation is unimportant to the overall dynamics of the system. Accordingly, acoustic wave propagation can be mathematically removed from the equations of motion, allowing for a numerical time step based on an advective CFL condition, and this leads to an increase in the allowable time step of order :math:`1/M` over an explicit, fully compressible method (:math:`M` is the Mach number).  In this mathematical framework, the total pressure is decomposed into the sum of a spatially constant (ambient) thermodynamic pressure :math:`P_0` and a perturbational pressure, :math:`\pi({\vec x})` that drives the flow.  Under suitable conditions, :math:`\pi/P_0 = \mathcal{O} (M^2)`. 
 
@@ -145,4 +145,27 @@ and simplified velocity constraint,
 
      \nabla \cdot \boldsymbol{u} = \delta S - \delta \theta \frac{\overline S}{\overline \theta} .
 
-    
+Geometry with Embedded Boundaries
+---------------------------------
+
+`PeleLMeX` relies on `AMReX's implementation<https://amrex-codes.github.io/amrex/docs_html/EB_Chapter.html>`_ of 
+the Embedded Boundaries (EB) approach to represent geometrical objects. In this approach, the underlying computational 
+mesh is uniform and block-structured, but the boundary of the irregular-shaped computational domain conceptually cuts 
+through this mesh. Each cell in the mesh becomes labeled as regular, cut or covered, and the finite-volume 
+based discretization methods traditionally used in AMReX applications need to be modified to incorporate these cell shapes.
+AMReX provides the necessary EB data structures, including volume and area fractions, surface normals and centroids, 
+as well as local connectivity information. The fluxes are then modified to account for the apperture opening between adjacent
+cells and the additional EB-fluxes are included when constructing the cell flux divergences.
+
+A common problem arising with EB is the presence of the small cut-cells which can either introduce undesirable constraint on 
+the explicit time step size or lead to numerical instabilities if not accounterd for. `PeleLMeX` relies on a combination of 
+classical flux redistribution (FRD) (Pember et al, 1995) and state redistribution (SRD) (Giuliani et al., 2022) to circumvent the issue.
+In particular, explicit advective fluxes are treated using SRD while explicit diffusion fluxes (appearing in the SDC context) 
+are treated with FRD. Note that implicit diffusion fluxes are not redistributed as AMReX's linear operators are EB-aware.
+
+The use of AMReX's multigrid linear solver introduces contraint on the complexity of the geometry `PeleLMeX` is able to handle. The
+efficiency of the multigrid appraoch relies on generating coarse version of the linear problem. If the geometry includes thin elements
+(such as tube or plate) or narrow channels, coarsening of the geometry is rapidly limited by the occurence of multi-cut cells (not 
+supported by AMReX) and the linear solvers are no longer able to robustly tackle projections and implicit diffusion solves. AMReX
+include an interface to HYPRE which can help circumvent the issue by sending the coarse-level geometry directly to HYPRE algebraic
+multigrid solvers. More details on how to use HYPRE is provided in control Section.

Docs/source/Validation.rst

@@ -1,17 +1,56 @@
-PeleLMeX validations
-====================
+PeleLMeX Verification & Validations
+===================================
 
 This section is work-in-progress.
 
+Laminar Poiseuille flow
+~~~~~~~~~~~~~~~~~~~~~~~
+
+The laminar pipe flow or Poiseuille flow, is a basic test case for wall bounded flows.
+In the present configuration, the geometry consist of a circular channel of radius :math:`R` = 1 cm 
+aligned with the :math:`x`-direction, where no-slip boundary conditions are imposed on
+EB surface. The flow is periodic in the :math:`x`-direction and a background pressure
+gradient :math:`dp /dx` is used to drive the flow.
+
+The exact solution at steady state is:
+
+.. math::
+   u(r) = \frac{G}{4 \mu} (R^2 - r^2)
+
+where :math:`G = -dp/dx`, and :math:`\mu` is the dynamic viscosity.
+The test case can be found in ``Exec/RegTests/EB_PipeFlow``, where
+the input parameters are very similar to the PeleC counterpart of 
+this case.
+
+The steady-state :math:`x`-velocity profiles accross the pipe diameter
+at increasing resolution are plotted along with the theorerical profile:
+
+.. figure:: images/validations/Poiseuille3D/PoiseuilleVelProf.png
+   :align: center
+   :figwidth: 60%
+
+A more quantitative evaluation of PeleLMeX results is obtained by calculating
+the L2 norm of the error against the analytical profile:
+
+
+.. figure:: images/validations/Poiseuille3D/PoiseuilleConvergence.png
+   :align: center
+   :figwidth: 60%
+
+showing second-order convergence for this diffusion dominated flow.
+
 Taylor-Green vortex breakdown
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The Taylor-Green vortex breakdown case is a classical CFD test case
-described in `here <>`
+described in `here <https://www1.grc.nasa.gov/research-and-engineering/hiocfd/>`
+(case C3.3).
 
 Building and running
 ####################
 
+The test case can be found in ``Exec/RegTests/TaylorGreen``.
+
 .. code-block:: bash
 
    $ make -j 16 DIM=3 USE_MPI=TRUE TPL

Docs/source/conf.py

@@ -22,16 +22,15 @@
 author = 'PeleTeam'
 
 # The full version, including alpha/beta/rc tags
-release = '21.11'
+release = '22.11'
 
 
 # -- General configuration ---------------------------------------------------
 
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
-extensions = ['sphinx.ext.autodoc'
-]
+extensions = ['sphinx.ext.autodoc']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']

Docs/source/index.rst

@@ -9,8 +9,7 @@ Welcome to PeleLMeX's documentation!
 `PeleLMeX` is the non-subcycling version of `PeleLM <https://amrex-combustion.github.io/PeleLM/>`_, an adaptive-mesh low Mach number hydrodynamics 
 code for reacting flows. If you need help or have questions, please join the users `forum <https://groups.google.com/forum/#!forum/pelelmusers>`_.
 The documentation pages appearing here are distributed with the code in the ``Docs`` folder as "restructured text" files.  The html is built
-automatically with certain pushes to the `PeleLM` GibHub repository and are maintained online by `ReadTheDocs <https://pelelmex.readthedocs.io/en/latest>`_.
-A local version can also be built as follows ::
+automatically with certain pushes to the `PeleLMeX` GibHub repository. A local version can also be built as follows ::
 
     cd ${PELELMEX_HOME}/Docs
     make html

Exec/RegTests/EB_PipeFlow/input.3d-Poiseuille

@@ -0,0 +1,68 @@
+#----------------------DOMAIN DEFINITION------------------------
+geometry.is_periodic = 1 0 0                # For each dir, 0: non-perio, 1: periodic
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ
+geometry.prob_lo     = 0.0  -0.01 -0.01        # x_lo y_lo (z_lo)
+geometry.prob_hi     = 0.04  0.01  0.01  # x_hi y_hi (z_hi)
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# Interior, Inflow, Outflow, Symmetry,
+# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm
+peleLM.lo_bc = Interior NoSlipWallAdiab SlipWallAdiab
+peleLM.hi_bc = Interior NoSlipWallAdiab SlipWallAdiab
+
+
+#-------------------------AMR CONTROL----------------------------
+amr.n_cell          = 32 16 16         # Level 0 number of cells in each direction   
+amr.v               = 1                # AMR verbose
+amr.max_level       = 0                # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2          # refinement ratio
+amr.regrid_int      = 2                # how often to regrid
+amr.n_error_buf     = 1 1 2 2          # number of buffer cells in error est
+amr.grid_eff        = 0.7              # what constitutes an efficient grid
+amr.blocking_factor = 8               # block factor in grid generation (min box size)
+amr.max_grid_size   = 64               # max box size
+
+#--------------------------- Problem -------------------------------
+prob.T_mean = 300.0
+prob.P_mean = 101325.0
+prob.meanFlowMag = 18.0
+prob.meanFlowDir = 1
+
+#-------------------------PeleLM CONTROL----------------------------
+peleLM.v = 2
+peleLM.incompressible = 1
+peleLM.rho = 1.0
+peleLM.mu = 0.0576
+peleLM.gradP0 = -82944.0 0.0 0.0
+peleLM.do_temporals = 1
+
+#amr.restart = chk01000
+amr.check_int = 500
+amr.plot_int = 100
+amr.max_step = 20000
+amr.dt_shrink = 1.0
+amr.stop_time = 0.015
+amr.cfl = 0.7
+amr.derive_plot_vars = avg_pressure mag_vort
+
+#------------------------- EB SETUP -----------------------------
+eb2.geom_type = cylinder
+eb2.cylinder_radius     = 0.01
+eb2.cylinder_direction  = 0
+eb2.cylinder_center     = 0.0 0.0 0.0
+eb2.cylinder.internal_flow    = true
+eb2.cylinder_has_fluid_inside = 1
+eb2.small_volfrac = 1.0e-4
+
+#--------------------REFINEMENT CONTROL------------------------
+amr.refinement_indicators = vort
+amr.vort.max_level     = 2
+amr.vort.vorticity_greater = 500
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+nodal_proj.verbose = 0
+nodal_proj.mg_max_coarsening_level = 2
+amrex.fpe_trap_invalid = 1
+amrex.fpe_trap_zero = 1
+amrex.fpe_trap_overflow = 1

Utils/Plots/plotPoiseuille.py

@@ -0,0 +1,91 @@
+import os
+import argparse
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.backends.backend_pdf import PdfPages
+
+# A script to plot data from the EB_PipeFlow case
+# when running a simple Poiseuille flow.
+# Flow properties are hardcoded to: G = 82944.0,
+# mu = 0.0576, radius = 0.01, umax = 36.0
+# Data from LMeX run at increasing resolution must be extracted
+# using fextract a-priori and stored in nr* folders.
+
+cmap = [
+    "#EE2E2F",
+    "#008C48",
+    "#185AA9",
+    "#F47D23",
+    "#662C91",
+    "#A21D21",
+    "#B43894",
+    "#010202",
+]
+
+def theory_ooa(order, res, orig):
+    return orig * (res[0] / res) ** order
+
+def eval_u_exact(r, dr):
+    return 82944.0 / (4.0 * 0.0576) * (0.01 ** 2 - r ** 2)
+
+if __name__ == "__main__":
+    radius = 0.01
+    umax = 36
+    r = np.linspace(-radius, radius, 200)
+    u_exact = eval_u_exact(r, 0.0001)
+
+    # Plot the exact data (pointwise, not cell-averaged)
+    plt.figure("u")
+    plt.rc('font', size=12)
+    plt.plot(
+        r / radius, u_exact/umax, lw=2, color=cmap[-1], label="Exact"
+    )
+
+    # Read in LMeX data
+    resolution=[8,16,32,64]
+
+    errors = np.zeros((2, len(resolution)))
+    for k, r in enumerate(resolution):
+        f = open('./nr{}/nr{}prof.dat'.format(r,r), 'r')
+        rlabel="res{}".format(r)
+        rad = []
+        u_x = []
+        for line in f:
+            line = line.strip()
+            columns = line.split()
+            rad.append(float(columns[0]))
+            u_x.append(float(columns[1]))
+        plt.plot(
+            np.array(rad)/radius, np.array(u_x)/umax, linestyle='--', lw=1, color=cmap[k], label=rlabel
+        )
+        errors[0, k] = r
+        errors[1, k] = np.sqrt(
+            np.sum(
+                (
+                    np.array(u_x) / umax
+                    - eval_u_exact(np.array(rad), 0.0001) / umax
+                )
+                ** 2
+            )
+            / r
+        )
+
+    plt.xlabel(r"$r / R$",fontsize=14)
+    plt.ylabel(r"$u / Umax$",fontsize=14)
+    plt.grid(which='both',color='#E6E3E3', linestyle=':', linewidth=1.0)
+    plt.legend(bbox_to_anchor=(0.5, 0.4), loc=1, borderaxespad=0.)
+    plt.savefig("PoiseuilleVelProf.png")
+
+    plt.figure("error")
+    plt.loglog(
+        errors[0, :],
+        errors[1, :],
+        label="PeleLMeX",
+    )
+    p2 = theory_ooa(2, errors[0, :], 0.95*errors[1, 0])
+    plt.loglog(errors[0, :], p2, color='k', linestyle='--', label=f"2nd-order")
+    plt.grid(which='both',color='#E6E3E3', linestyle=':', linewidth=1.0)
+    plt.xlabel("Resolution")
+    plt.ylabel("Error L2-norm")
+    plt.legend(bbox_to_anchor=(0.9, 0.9), loc=1, borderaxespad=0.)
+    plt.savefig("PoiseuilleConvergence.png")