move redist into amrex (AMReX-Codes#3378)

This PR moves all redistribution functionality from AMReX-Hydro to amrex/Src/EB. It also adds re-redistribution to the routines in EB (copied from Tests/EB_CNS). (It also moves the slope routines since those are needed by state redistribution.) This change will require changes in any code that currently uses AMReX-Hydro for redistribution. Additional changes will be required in application codes that want to enable re-redistribution. Once this PR is merged in amrex, all routines in AMReX-Hydro/Slopes and AMReX-Hydro/Redistribution will be deleted. Co-authored-by: Weiqun Zhang <[email protected]>
IlianCS · Jul 3, 2023 · 9fc6354 · 9fc6354
1 parent 22a1a3a
commit 9fc6354
Show file tree

Hide file tree

Showing 47 changed files with 5,484 additions and 752 deletions.
diff --git a/Docs/sphinx_documentation/source/EB.rst b/Docs/sphinx_documentation/source/EB.rst
@@ -4,7 +4,6 @@
 .. role:: fortran(code)
    :language: fortran
 
-
 .. _sec:EB:ebinit:
 
 Initializing the Geometric Database
@@ -387,6 +386,174 @@ testing cell types and getting neighbor information. For example
         end do
     end do
 
+Small Cell Problem and Redistribution
+=====================================
+
+First, we review finite volume discretizations with embedded boundaries as used by
+AMReX-based applications. Then we illustrate the small cell problem.
+
+Finite Volume Discretizations
+-----------------------------
+
+Consider a system of PDEs to advance a conserved quantity :math:`U` with fluxes
+:math:`F`:
+
+.. math:: \frac{\partial U}{\partial t} + \nabla \cdot F = 0.
+  :label: eqn::hypsys
+
+A conservative, finite volume discretization starts with the divergence theorm
+
+.. math:: \int_V \nabla \cdot F dV = \int_{\partial V} F \cdot n dA.
+
+In an embedded boundary cell, the "conservative divergence" is discretized (as
+:math:`D^c(F)`) as follows
+
+.. math::
+  :label: eqn::ebdiv
+
+   D^c(F) = \frac{1}{\kappa h} \left( \sum^D_{d = 1}
+     (F_{d, \mathrm{hi}} \, \alpha_{d, \mathrm{hi}} - F_{d, \mathrm{lo}}\, \alpha_{d, \mathrm{lo}})
+     + F^{EB} \alpha^{EB} \right).
+
+Geometry is discretely represented by volumes (:math:`V = \kappa h^d`) and
+apertures (:math:`A= \alpha h^{d-1}`), where :math:`h` is the (uniform) mesh
+spacing at that AMR level, :math:`\kappa` is the volume fraction and
+:math:`\alpha` are the area fractions. Without multivalued cells the volume
+fractions, area fractions and cell and face centroids (see
+:numref:`fig::volume`) are the only geometric information needed to compute
+second-order fluxes centered at the face centroids, and to infer the
+connectivity of the cells. Cells are connected if adjacent on the Cartesian
+mesh, and only via coordinate-aligned faces on the mesh. If an aperture,
+:math:`\alpha = 0`, between two cells, they are not directly connected to each
+other.
+
+.. raw:: latex
+
+   \begin{center}
+
+.. |a| image:: ./EB/areas_and_volumes.png
+       :width: 100%
+
+.. |b| image:: ./EB/eb_fluxes.png
+       :width: 100%
+
+.. _fig::volume:
+
+.. table:: Illustration of embedded boundary cutting a two-dimensional cell.
+   :align: center
+
+   +-----------------------------------------------------+------------------------------------------------------+
+   |                        |a|                          |                        |b|                           |
+   +-----------------------------------------------------+------------------------------------------------------+
+   | | A typical two-dimensional uniform cell that is    | | Fluxes in a cut cell.                              |
+   | | cut by the embedded boundary. The grey area       | |                                                    |
+   | | represents the region excluded from the           | |                                                    |
+   | | calculation. The portion of the cell faces        | |                                                    |
+   | | faces (labelled with A) through which fluxes      | |                                                    |
+   | | flow are the "uncovered" regions of the full      | |                                                    |
+   | | cell faces. The volume (labelled V) is the        | |                                                    |
+   | | uncovered region of the interior.                 | |                                                    |
+   +-----------------------------------------------------+------------------------------------------------------+
+
+.. raw:: latex
+
+   \end{center}
+
+
+Small Cells And Stability
+-------------------------
+
+In the context of time-explicit advance methods for, say hyperbolic
+conservation laws, a naive discretization in time of :eq:`eqn::hypsys` using
+:eq:`eqn::ebdiv`,
+
+.. math:: U^{n+1} = U^{n} - \delta t D^c(F)
+
+would have a time step constraint :math:`\delta t \sim h \kappa^{1/D}/V_m`,
+which goes to zero as the size of the smallest volume fraction :math:`\kappa` in
+the calculation. Since EB volume fractions can be arbitrarily small, this presents an
+unacceptable constraint. This is the so-called "small cell problem," and AMReX-based
+applications address it with redistribution methods.
+
+Flux Redistribution
+-----------------------------
+
+Consider a conservative update in the form:
+
+.. math:: (\rho \phi)_t + \nabla \cdot ( \rho \phi u) = RHS
+
+For each valid cell in the domain, compute the conservative divergence, :math:`(\nabla \cdot F)^c` ,
+of the convective fluxes, :math:`F`
+
+.. math:: (\nabla \cdot {F})^c_i = \dfrac{1}{\mathcal{V}_i} \sum_{f=1}^{N_f} ({F}_f\cdot{n}_f) A_f
+
+Here :math:`N_f` is the number of faces of cell :math:`i`, :math:`\vec{n}_f` and :math:`A_f`
+are the unit normal and area of the :math:`f` -th face respectively,
+and :math:`\mathcal{V}_i` is the volume of cell :math:`i` given by
+
+.. math:: \mathcal{V}_i = (\Delta x \Delta y \Delta z)\cdot \mathcal{K}_i
+
+where :math:`\mathcal{K}_i` is the volume fraction of cell :math:`i` .
+
+Now, a conservative update can be written as
+
+.. math:: \frac{ \rho^{n+1} \phi^{n+1} - \rho^{n} \phi^{n} }{\Delta t} = - \nabla \cdot{F}^c
+
+For each cell cut by the EB geometry, compute the non-conservative update, :math:`\nabla \cdot {F}^{nc}` ,
+
+.. math:: \nabla\cdot{F}^{nc}_i = \dfrac{\sum\limits_{j\in N(i) } \mathcal{K}_j\nabla \cdot {F}^c_j} {\sum\limits_{j\in N(i) } {\mathcal{K}}_j}
+
+where :math:`N(i)` is the index set of cell :math:`i` and its neighbors.
+
+For each cell cut by the EB geometry, compute the convective update :math:`\nabla \cdot{F}^{EB}` follows:
+
+.. math:: \nabla \cdot{F}^{EB}_i = \mathcal{K}_i\nabla \cdot{F}^{c}_i +(1-\mathcal{K}_i) \nabla \cdot \mathcal{F}^{nc}_i
+
+For each cell cut by the EB geometry, redistribute its mass loss, :math:`\delta M_i` , to its neighbors:
+
+.. math::  \nabla \cdot {F}^{EB}_j :=   \nabla \cdot {F}^{EB}_j + w_{ij}\delta M_i\, \qquad \forall j\in N(i)\setminus i
+
+where the mass loss in cell :math:`i` , :math:`\delta M_i` , is given by
+
+.. math:: \delta M_i =  \mathcal{K}_i(1- \mathcal{K}_i)[ \nabla \cdot {F}^c_i-  \nabla \cdot {F}^{nc}_i]
+
+and the weights, :math:`w_{ij}` , are
+
+.. math:: w_{ij} = \dfrac{1}{\sum\limits_{j\in N(i)\setminus i}  \mathcal{K}_j}
+
+Note that :math:`\nabla \cdot{F}_i^{EB}` gives an update for :math:`\rho \phi` ; i.e.,
+
+.. math:: \frac{(\rho \phi_i)^{n+1} - (\rho \phi_i)^{n} }{\Delta t} = - \nabla \cdot{F}^{EB}_i
+
+Typically, the redistribution neighborhood for each cell is one that can be
+reached via a monotonic path in each coordinate direction of unit length (see,
+e.g., :numref:`fig::redistribution`)
+
+.. raw:: latex
+
+   \begin{center}
+
+.. _fig::redistribution:
+
+.. figure:: ./EB/redist.png
+   :width: 50.0%
+
+   : Redistribution illustration. Excess update distributed to neighbor cells.
+
+.. raw:: latex
+
+   \end{center}
+
+
+State Redistribution
+-----------------------------
+
+For state redistribution we implement the weighted state
+redistribution algorithm as described in Guiliani et al (2021),
+which is available on `arxiv  <https://arxiv.org/abs/2112.12360>`_ .
+This is an extension of the original state redistribution algorithm
+of Berger and Guiliani (2020).
+
 
 Linear Solvers
 ==============

diff --git a/Docs/sphinx_documentation/source/EB/areas_and_volumes.fig b/Docs/sphinx_documentation/source/EB/areas_and_volumes.fig
@@ -0,0 +1,69 @@
+#FIG 3.2  Produced by xfig version 3.2.5c
+Landscape
+Center
+Metric
+A4      
+100.00
+Single
+-2
+1200 2
+2 2 0 7 0 7 14 -1 -1 0.000 0 0 -1 0 0 5
+	 1800 1800 5400 1800 5400 5400 1800 5400 1800 1800
+2 1 0 7 0 7 14 -1 15 0.000 0 0 -1 0 0 4
+	 1800 2925 4275 5400 1800 5400 1800 2925
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 4230 5580 4230 6570
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 5400 5625 5400 6525
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 5625 1800 7650 1800
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 5625 5400 7650 5400
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 1575 1800 225 1800
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 1575 2925 225 2925
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 1800 1575 1800 450
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 5400 1575 5400 225
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 1852 2845 2430 2340
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 0 0 2
+	 4372 5275 4950 4770
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 3060 3510 2160 2610
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 3690 4140 4590 5040
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 3420 6120 4230 6120
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 6300 6120 5400 6120
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 6300 3420 6300 1890
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 6300 4140 6300 5400
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 810 990 810 1890
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 810 3690 810 2970
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 3060 1170 1800 1170
+2 1 0 3 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 4050 1170 5400 1170
+4 0 0 14 -1 18 18 0.0000 4 270 375 4500 6210 Ay\001
+4 0 0 14 -1 18 18 0.0000 4 210 375 6120 3780 Ax\001
+4 0 0 14 -1 18 18 0.0000 4 270 375 3330 1260 Ay\001
+4 0 0 14 -1 18 18 0.0000 4 210 375 720 2340 Ax\001
+4 0 0 14 -1 18 18 0.0000 4 210 555 3150 3780 Aeb\001
+4 0 0 14 -1 18 18 0.0000 4 210 195 3780 2790 V\001
diff --git a/Docs/sphinx_documentation/source/EB/areas_and_volumes.png b/Docs/sphinx_documentation/source/EB/areas_and_volumes.png
diff --git a/Docs/sphinx_documentation/source/EB/eb_fluxes.fig b/Docs/sphinx_documentation/source/EB/eb_fluxes.fig
@@ -0,0 +1,38 @@
+#FIG 3.2  Produced by xfig version 3.2.5c
+Landscape
+Center
+Metric
+A4      
+100.00
+Single
+-2
+1200 2
+2 2 0 7 0 7 14 -1 -1 0.000 0 0 -1 0 0 5
+	 1800 1800 5400 1800 5400 5400 1800 5400 1800 1800
+2 1 0 7 0 7 14 -1 15 0.000 0 0 -1 0 0 4
+	 1800 2925 4275 5400 1800 5400 1800 2925
+2 1 0 4 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 2925 4275 4050 3150
+2 1 0 4 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 720 2475 2745 2475
+2 1 0 4 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 4725 6435 4725 4410
+2 1 0 4 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 4725 3600 6525 3600
+2 1 0 4 0 7 14 -1 -1 0.000 0 0 -1 1 0 2
+	1 1 5.00 135.00 165.00
+	 3600 2475 3600 450
+4 0 0 14 -1 18 28 0.0000 4 345 285 4050 900 F\001
+4 0 0 14 -1 18 28 0.0000 4 345 285 4050 3150 F\001
+4 0 0 14 -1 18 18 0.0000 4 210 405 4320 3240 EB\001
+4 0 0 14 -1 18 18 0.0000 4 270 930 4320 990 Y, high\001
+4 0 0 14 -1 18 28 0.0000 4 345 285 -270 2475 F\001
+4 0 0 14 -1 18 18 0.0000 4 255 735 45 2520 X,low\001
+4 0 0 14 -1 18 28 0.0000 4 345 285 4455 6885 F\001
+4 0 0 14 -1 18 18 0.0000 4 255 825 4725 6975 Y, low\001
+4 0 0 14 -1 18 18 0.0000 4 270 840 6930 3915 X,high\001
+4 0 0 14 -1 18 28 0.0000 4 345 285 6660 3825 F\001
diff --git a/Docs/sphinx_documentation/source/EB/eb_fluxes.png b/Docs/sphinx_documentation/source/EB/eb_fluxes.png
diff --git a/Docs/sphinx_documentation/source/EB/redist.fig b/Docs/sphinx_documentation/source/EB/redist.fig
@@ -0,0 +1,41 @@
+#FIG 3.2  Produced by xfig version 3.2.5c
+Landscape
+Center
+Metric
+A4      
+100.00
+Single
+-2
+1200 2
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 1800 1755 1800 9900
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 3600 1755 3600 9900
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 5400 1755 5400 9900
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 7200 1755 7200 9900
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 0 1755 0 9900
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 8100 6300 0 6300
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 8325 4500 0 4500
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 8550 2700 0 2700
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 7875 9900 0 9900
+2 1 0 7 0 7 56 -1 12 0.000 0 0 -1 0 0 11
+	 0 2925 0 9900 7650 9900 7200 9000 6300 8100 5400 7650
+	 3600 7200 2475 6300 1800 5175 1350 4500 0 2925
+2 1 0 5 0 7 50 -1 -1 0.000 0 0 -1 0 0 2
+	 8100 8100 0 8100
+2 1 0 3 0 7 56 -1 -1 0.000 0 0 -1 1 0 2
+	0 0 3.00 180.00 360.00
+	 3150 6750 2700 5400
+2 1 0 3 0 7 56 -1 -1 0.000 0 0 -1 1 0 2
+	0 0 3.00 180.00 360.00
+	 3150 6750 4500 5175
+2 1 0 3 0 7 56 -1 -1 0.000 0 0 -1 1 0 2
+	0 0 3.00 180.00 360.00
+	 3150 6750 4725 6750
diff --git a/Docs/sphinx_documentation/source/EB/redist.png b/Docs/sphinx_documentation/source/EB/redist.png
diff --git a/Docs/sphinx_documentation/source/LinearSolvers.rst b/Docs/sphinx_documentation/source/LinearSolvers.rst
@@ -509,7 +509,7 @@ passed to the constructor of a linear operator to disable the
 coarsening completely.  In that case the bottom solver is solving the
 residual correction form of the original problem. To build Hypre, follow the next steps:
 
-.. highlight:: c++
+.. highlight:: console
 
 ::
 
@@ -577,7 +577,7 @@ The user is referred to the
 AMReX can also use `PETSc <https://www.mcs.anl.gov/petsc/>`_ as a bottom solver for cell-centered
 problems. To build PETSc, follow the next steps:
 
-.. highlight:: c++
+.. highlight:: console
 
 ::
 

diff --git a/Src/Base/AMReX_FArrayBox.H b/Src/Base/AMReX_FArrayBox.H
@@ -3,13 +3,9 @@
 #define BL_FARRAYBOX_H
 #include <AMReX_Config.H>
 
-#include <AMReX_Box.H>
 #include <AMReX_BaseFab.H>
 #include <AMReX_FabConv.H>
-#include <AMReX_REAL.H>
-#include <AMReX_SPACE.H>
 #include <AMReX_FabFactory.H>
-#include <cmath>
 
 namespace amrex {