Change default to not use slopes in SRD for robustness (#698)

Docs/sphinx/geometry/EB.rst

@@ -12,21 +12,21 @@
 Geometry treatment in PeleC
 ===========================
 
-The treatment of geometric features that do not align along cartesian coordinate directions effectively reduces to 
+The treatment of geometric features that do not align along cartesian coordinate directions effectively reduces to
 determining the correct flux terms at cut-cell interfaces and subsequent update of divergence term in each cell.
-This involves the initialization and query of the necessary AMReX-provided data structures containing the 
-geometry information, and computation of PeleC-specific advection and diffusion operators. The various steps in the 
+This involves the initialization and query of the necessary AMReX-provided data structures containing the
+geometry information, and computation of PeleC-specific advection and diffusion operators. The various steps in the
 process are:
 
-1. Creation of a functional specification of the irregular geometry to embed in the uniform grid. This is done via exact 
+1. Creation of a functional specification of the irregular geometry to embed in the uniform grid. This is done via exact
    function representations of the geometry or implicit functions.
-2. Construction of map of the (continuous) implicit representation of geometry onto the discrete mesh on all AMR levels.  
+2. Construction of map of the (continuous) implicit representation of geometry onto the discrete mesh on all AMR levels.
    This will be a large, complex, distributed data structure.
 3. Communication of the subsets of this large data set to the local cores tasked with building the PeleC operators.
 4. Actual construction of the diffusion and advection components of the PeleC time advance.
 
-AMReX data structures and functions provide for the first 3 steps.  
-Step 4 is implemented using a "method-of-lines" update within PeleC (see section :ref:`MOL<MOL>`). 
+AMReX data structures and functions provide for the first 3 steps.
+Step 4 is implemented using a "method-of-lines" update within PeleC (see section :ref:`MOL<MOL>`).
 
 Embedded Boundary Representation
 --------------------------------
@@ -40,7 +40,7 @@ Embedded Boundary Representation
    Embedded boundary representation of geometry
 
 
-Geometry is treated in PeleC using an embedded boundary (EB) formulation, based on datastructures and algorithmic components provided by AMReX.   In the EB formalism, geometry is represented by volume fractions (:math:`v_l`) 
+Geometry is treated in PeleC using an embedded boundary (EB) formulation, based on datastructures and algorithmic components provided by AMReX.   In the EB formalism, geometry is represented by volume fractions (:math:`v_l`)
 and apertures (:math:`A_l^k`) for each cell :math:`l` that have faces :math:`1,...,k,6`. See :ref:`EB_F` for an illustration where the grey area represents the region excluded from the solution domain and the arrows represent fluxes. The fluid volume in a given cell is given by  (:math:`V_l = v_l\,\,dx\,dy\,dz`=v_l\,dx^3`); it should be noted that the grid spacing along each direction is the same in PeleC.
 
 .. _EB_F:
@@ -60,7 +60,7 @@ and apertures (:math:`A_l^k`) for each cell :math:`l` that have faces :math:`1,.
    Embedded boundary representation of geometry
 
 
-The geometry components in AMReX are used in PeleC to implement a time-explicit integrator based on the method-of-lines.  For the advection and diffusion components of the PeleC time integrator, the time rate of change of the conserved fields, S, in cell :math:`l` can be written as 
+The geometry components in AMReX are used in PeleC to implement a time-explicit integrator based on the method-of-lines.  For the advection and diffusion components of the PeleC time integrator, the time rate of change of the conserved fields, S, in cell :math:`l` can be written as
 
 .. math::
   \frac{dS_l}{dt} = \nabla \cdot F
@@ -71,20 +71,21 @@ Redistribution
 --------------
 .. _Redist:
 
-A straightforward implemention of the finite-volume advance of intensive conserved fields is numerically unstable (this is the well-known "small cell issue") due to presence of 
+A straightforward implemention of the finite-volume advance of intensive conserved fields is numerically unstable (this is the well-known "small cell issue") due to presence of
 the fluid cell volume in the denominator of the conservative divergence (:math:`(DC)_l`):
 
 .. math::
   (DC)_l = \frac{1}{V_l} \sum_{k_l} \left( F_k \cdot n_k A_k \right),
 
 where :math:`k_l` is the number of regular and cut faces surrounding cell :math:`l` and :math:`F_k` is the intensive flux at the centroid of face :math:`k`.
 
-There are a number of ways to deal with this "small cell issue" and the reader is referred to the relevant discussion in `Berger, Marsha, and Andrew Giuliani. "A state redistribution algorithm for finite volume schemes on cut cell meshes." Journal of Computational Physics 428 (2021): 109820 <https://doi.org/10.1016/j.jcp.2020.109820>`_. PeleC supports the different types of redistributions described in the paper with the keyword ``pelec.redistribution_type``, which can have the following values:
+There are a number of ways to deal with this "small cell issue" and the reader is referred to the relevant discussion in `Berger, Marsha, and Andrew Giuliani. "A state redistribution algorithm for finite volume schemes on cut cell meshes." Journal of Computational Physics 428 (2021): 109820 <https://doi.org/10.1016/j.jcp.2020.109820>`_ and `Giuliani, Andrew et al. "A weighted state redistribution algorithm for embedded boundary grids" Journal of Computational Physics 464 (2022): 111305 <https://doi.org/10.1016/j.jcp.2022.111305>`_. PeleC supports the different types of redistributions described in the papers and implemented in AMReX, using the keyword ``pelec.redistribution_type``, which can have the following values:
 
 * ``"NoRedist"``: no redistribution
 * ``"FluxRedist"``: flux redistribution
-* ``"StateRedist"``: state redistribution
-* ``"NewStateRedist"``: improved state redistribution (default)
+* ``"StateRedist"``: weighted state redistribution
+
+The "StateRedist" approach is the default and recommended approach. Part of this algorithm involves using limited, least squares slopes to reconstruct the state within merged cells. The slopes for each state variable are limited individually, which may introduce inconsistency between the internal energy and density that leads to negative temperatures. Therefore, these slopes are set to zero by default in PeleC for robustness. Advanced users wishing to use this feature to maintain 2nd order accuracy at the EB surface may do so by setting ``pelec.eb_srd_max_order = 2`` in their input file.
 
 Mass fractions at the EB
 ------------------------
@@ -108,22 +109,22 @@ Re-redistribution
 
    \(a) an example situation with an EB spanning a coarse-fine boundary, (b) same situation as seen by the coarse level and (c) same situation as seen by the fine level. The cells with the dotted lines are ghost cells.
 
-The redistribution of mass with the use of hybrid divergence method leads 
-to an accounting problem at coarse-fine interfaces that have an EB passing through them, as shown in :ref:`re-redistribution figure<eb_re_redist>` (a). 
-The correct strategy will be to redistribute mass from the coarse mesh on the left side to the fine mesh on the right and vice-versa, when divergence is evaluated on the fly. 
-This strategy is difficult to implement directly into the current algorithmic framework, because flux/residual calculation and time advance are done separately at 
-each level with a ghost-cell treatment at coarse-fine boundaries. Therefore the mass distributed to and from ghost-cells need to be accounted and adjusted after each level has 
+The redistribution of mass with the use of hybrid divergence method leads
+to an accounting problem at coarse-fine interfaces that have an EB passing through them, as shown in :ref:`re-redistribution figure<eb_re_redist>` (a).
+The correct strategy will be to redistribute mass from the coarse mesh on the left side to the fine mesh on the right and vice-versa, when divergence is evaluated on the fly.
+This strategy is difficult to implement directly into the current algorithmic framework, because flux/residual calculation and time advance are done separately at
+each level with a ghost-cell treatment at coarse-fine boundaries. Therefore the mass distributed to and from ghost-cells need to be accounted and adjusted after each level has
 advanced a single time step, which we refer to as re-redistribution. Specifically, four different mass terms need to be accounted for:
 
 *   In :ref:`re-redistribution figure<eb_re_redist>` (b) the left-coarse-real-cell distributes mass to the right-coarse-ghost-cell. This needs to be captured and given to the right-fine-real-cells.
-*   In :ref:`re-redistribution figure<eb_re_redist>` (b) the right-coarse-ghost-cell distributes mass to the left-coarse-real-cell. This needs to be captured and removed from the left-coarse-real-cell update 
+*   In :ref:`re-redistribution figure<eb_re_redist>` (b) the right-coarse-ghost-cell distributes mass to the left-coarse-real-cell. This needs to be captured and removed from the left-coarse-real-cell update
     because the correct distributed mass has to come from the right-fine-real-cells.
 *   In :ref:`re-redistribution figure<eb_re_redist>` (c) the right-fine-real-cells distribute mass to the left-fine-ghost-cells. This needs to be captured and given to the left-coarse-real-cell.
-*   In :ref:`re-redistribution figure<eb_re_redist>` (c) the left-fine-ghost-cells distribute mass to the right-fine-real-cells. This needs to be captured and removed from the right-fine-real-cells update 
+*   In :ref:`re-redistribution figure<eb_re_redist>` (c) the left-fine-ghost-cells distribute mass to the right-fine-real-cells. This needs to be captured and removed from the right-fine-real-cells update
     because the correct distributed mass has to come from the left-coarse-real-cell.
 
-The re-redistribution is implemented as a book-keeping step where the mass distributed are stored during MOL divergence calculation and given to the coarse and fine flux registers to reflux at 
-the end of each time step. The re-redistribution is performed every time the reflux function is called in post_timestep. More details regarding re-redistribution are 
+The re-redistribution is implemented as a book-keeping step where the mass distributed are stored during MOL divergence calculation and given to the coarse and fine flux registers to reflux at
+the end of each time step. The re-redistribution is performed every time the reflux function is called in post_timestep. More details regarding re-redistribution are
 presented in `Pember et al. <https://www.sciencedirect.com/science/article/pii/S0021999185711655>`_. A forthcoming paper will describe the methodology for this procedure when using state redistribution.
 
 
@@ -152,10 +153,10 @@ An array of structures is created on level creation by copying data from the AMR
    :width: 500
 
 
-   Storage for sparse EB structures 
+   Storage for sparse EB structures
+
 
-
-On creation of a new AMRLevel, data is cached from the *dense* AMReX structures in the *sparse* PeleC structures. 
+On creation of a new AMRLevel, data is cached from the *dense* AMReX structures in the *sparse* PeleC structures.
 
 
 Applying boundary and face stencils

Source/Params/_cpp_parameters

@@ -217,7 +217,7 @@ eb_clean_massfrac           bool          true
 eb_clean_massfrac_threshold  Real         0.0
 
 # Max order used for SRD slopes
-eb_srd_max_order             int          2
+eb_srd_max_order             int          0
 
 # Weight types for ML redistribution (0 = 1.0, 1 = energy, 2 = density, 3 = vfrac)
 eb_weights_type              int          2

Source/Params/param_includes/pelec_defaults.H

@@ -64,7 +64,7 @@ bool PeleC::eb_noslip = true;
 std::string PeleC::redistribution_type = "StateRedist";
 bool PeleC::eb_clean_massfrac = true;
 amrex::Real PeleC::eb_clean_massfrac_threshold = 0.0;
-int PeleC::eb_srd_max_order = 2;
+int PeleC::eb_srd_max_order = 0;
 int PeleC::eb_weights_type = 2;
 bool PeleC::do_mms = false;
 std::string PeleC::masa_solution_name = "ad_cns_3d_les";