Merge branch 'development'

.github/workflows/sedov-compare.yml

@@ -0,0 +1,39 @@
+name: Sedov
+
+on: [pull_request]
+jobs:
+  Sedov-3d:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+        with:
+          fetch-depth: 0
+
+      - name: Get submodules
+        run: |
+          git submodule update --init
+          cd external/Microphysics
+          git fetch; git checkout development
+          cd ../amrex
+          git fetch; git checkout development
+          cd ../..
+
+      - name: Install dependencies
+        run: |
+          sudo apt-get update -y -qq
+          sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0 gfortran>=9.3.0
+
+      - name: Compile Sedov
+        run: |
+          cd Exec/hydro_tests/Sedov
+          make DEBUG=TRUE USE_MPI=FALSE -j 4
+
+      - name: Run Sedov
+        run: |
+          cd Exec/hydro_tests/Sedov
+          ./Castro3d.gnu.DEBUG.ex inputs.3d.sph max_step=10 amr.max_level=2 amr.plot_files_output=0 amr.checkpoint_files_output=0
+
+      - name: Compare to stored output
+        run: |
+          cd Exec/hydro_tests/Sedov
+          diff grid_diag.out ci-benchmarks/grid_diag.out

CHANGES.md

@@ -1,3 +1,27 @@
+# 22.06
+
+   * castro.stopping_criterion_field and castro.stopping_criterion_value have
+     been added; these allow you to stop the simulation once a certain threshold
+     has been exceeded (for example, if the temperature gets too hot). (#2209)
+
+   * The option castro.show_center_of_mass has been removed. If castro.v = 1
+     and castro.sum_interval > 0, then the center of mass will automatically
+     be included with the other diagnostic sums that are displayed. (#2176)
+
+   * The option castro.state_in_rotating_frame has been removed. The default
+     behavior continues to be that when rotation is being used, fluid variables
+     are measured with respect to the rotating frame. (#2172)
+
+   * The default for use_pslope has been changed to 0 -- disabling this.
+     use_pslope enables reconstruction that knows about HSE for the PLM
+     (castro.ppm_type = 0) implementation.  Since that method is not the
+     default, it is unlikely that this has been used.  This change is being
+     done to allow for a PPM implementation to be added without changing
+     the default behavior of that method. (#2205)
+
+   * The ``castro.use_pslope`` functionality to well-balance HSE has been
+     extended to PPM (#2202)
+
 # 22.05
 
    * A new option castro.hydro_memory_footprint_ratio has been added which

Docs/source/FlowChart.rst

@@ -63,7 +63,7 @@ The time-integration method used is controlled by
     described above that uses the CTU hydro advection and an ODE
     reaction solve.  Note: because this requires a different set of
     state variables, you must compile with ``USE_SIMPLIFIED_SDC = TRUE`` for this
-    method to work (in particular, this defines ``PRIM_SPECIES_HAVE_SOURCES``).
+    method to work.
 
 .. index:: USE_SIMPLIFIED_SDC, USE_TRUE_SDC
 

Docs/source/Hydrodynamics.rst

@@ -568,7 +568,7 @@ There are four major steps in the hydrodynamics update:
 
 #. Doing the conservative update
 
-.. index:: castro.do_hydro, castro.add_ext_src, castro.do_sponge, castro.normalize_species, castro.spherical_star, castro.show_center_of_mass
+.. index:: castro.do_hydro, castro.add_ext_src, castro.do_sponge, castro.normalize_species, castro.spherical_star
 
 Each of these steps has a variety of runtime parameters that
 affect their behavior. Additionally, there are some general
@@ -597,8 +597,6 @@ runtime parameters for hydrodynamics:
    function is then used to set the values outside the domain in
    implementing the boundary conditions.
 
--  ``castro.show_center_of_mass``: (0 or 1; default: 0)
-
 .. index:: castro.small_dens, castro.small_temp, castro.small_pres
 
 Several floors are imposed on the thermodynamic quantities to prevet unphysical

Docs/source/faq.rst

@@ -172,10 +172,10 @@ Managing Runs
    .. note::
 
       The parameter ``amr.message_int`` controls how often the
-      existence of these files is checked; by default it is 10, so the
-      check will be done at the end of every timestep that is a
-      multiple of 10.  Set that to 1 in your inputs file if you’d like
-      it to check every timestep.
+      existence of these files is checked; by default it is 1, so the
+      check will be done at the end of every timestep, but you can
+      set it to some other integer to check only timesteps that are a
+      multiple of that number.
 
 #. *How can I output plotfiles in single precision?*
 

Docs/source/hse.rst

@@ -0,0 +1,169 @@
+***********************
+Hydrostatic Equilibrium
+***********************
+
+There are a few places where we take care to maintain HSE.  These differ
+by integration method:
+
+* CTU: this does characteristic tracing of the interface states and is
+  used by Strang and simplified-SDC.  CTU requires that we predict the
+  interface states at the midpoint in time -- this can make enforcing
+  HSE much more difficult.
+
+* true-SDC: this uses a method of lines time reconstruction of the interface
+  states -- they are only spatially reconstructed, not extrapolated in time.
+
+
+Reconstruction
+==============
+
+Piecewise linear
+----------------
+
+.. index:: castro.use_pslope, castro.pslope_cutoff_density
+
+
+Piecewise linear reconstruction can be used with both CTU and
+true-SDC.
+
+For the 4th-order MC limiter, if we set ``castro.use_pslope = 1``,
+then we subtract off the hydrostatic pressure before limiting.
+The idea here is that the hydrostatic pressure balances gravity
+so it is not available to generate waves or cause over/undershoots.
+We define the excess pressure as :math:`\tilde{p}`:
+
+.. math::
+
+   \tilde{p}_i = p_i - p_{\mathrm{hse},i}
+
+This is designed to work with HSE discretized as:
+
+.. math::
+
+   p_{i+1} = p_i + \frac{1}{4} \Delta x (\rho_i + \rho_{i+1}) (g_i + g_{i+1})
+
+The 4th order MC limiter uses information in zones
+:math:`\{i-2,i-1,i,i+1,i+2\}`.  We pick zone ``i`` as the reference
+and integrate from ``i`` to the 4 other zone centers, and construct
+:math:`\{\tilde{p}_{i-2}, \tilde{p}_{i-1}, \tilde{p}_{i}, \tilde{p}_{i+1}, \tilde{p}_{i+2}\}`.  We then limit on this, giving us the change in the excess
+pressure over the zone, :math:`\Delta \tilde{p}_i`.  Finally, we 
+construct the total slope (including the hydrostatic part) as:
+
+.. math::
+
+   \Delta p_i = \Delta \tilde{p}_i + \rho_i g_i \Delta x
+
+.. note::
+
+   This can be disabled at low densities by setting ``castro.pslope_cutoff_density``.
+
+The remainder of the PLM algorithm is unchanged.  Since this was just
+the reconstruction part, this step is identical between CTU and true-SDC.
+
+
+
+PPM
+---
+
+The PPM version of ``use_pslope`` is essentially the same as the PLM
+version described above.  We first compute the excess pressure,
+:math:`\tilde{p}`, and then fit a cubic to the 4 cells surrounding an
+interface to get the initial interface values.  These become the left
+and right values of the parabola in each zone.  The PPM limiting is
+then done on the parabola, again working with :math:`\tilde{p}`.
+Finally, the parabola values are updated to include the hydrostatic
+pressure.
+
+
+Fully fourth-order method
+-------------------------
+
+The 4th order accurate true SDC solver does not appear to need any
+well-balancing.  It maintains HSE to very high precision, as shown in
+:cite:`castro-sdc`.
+
+
+Boundary conditions
+===================
+
+reflecting
+----------
+
+For piecewise linear reconstruction, in ``slope.H``, when using the
+4th order MC slopes (``castro.plm_limiter = 2``), we impose special
+ghost cell values on the normal velocity at reflecting boundaries to
+ensure that the velocity goes to 0 at the boundary. This follows
+:cite:`saltzman:1994` , page 162 (but note that they have a sign
+error).  Only the normal velocity is treated specially.
+
+
+
+
+HSE
+---
+
+.. index:: castro.hse_zero_vels, castro.hse_reflect_vels, castro.hse_interp_temp, castro.hse_fixed_temp
+
+For hydrostatic boundary conditions, we follow the method from
+:cite:`ppm-hse`.  Essentially, this starts with the last
+zone inside of the boundary and integrates HSE into the ghost cells,
+keeping the density either constant or extrapolating it.  Together
+the discretized HSE equation and the EOS yield the density and pressure
+in the ghost cells.
+
+.. warning::
+
+   The HSE boundary condition only works with constant gravity at the moment.
+
+To enable this, we set the appropriate boundary's ``xl_ext_bc_type``, ``xr_ext_bc_type``,
+``yl_ext_bc_type``, ``yr_ext_bc_type``, ``zl_ext_bc_type``,
+``zr_ext_bc_type`` to ``1``.
+
+We then control the behavior via the following options.
+
+For the velocity, we have:
+
+* ``hse_zero_vels`` : all 3 components of the velocity in the ghost
+  cell are set to ``0``.
+
+* ``hse_reflect_vels`` : the normal velocity is reflected and the transverse
+  velocity components are given a zero gradient.
+
+If neither of these are set, then all components of the velocity are
+simply given a zero gradient.
+
+The temperature in the ghost cells is controlled by:
+
+* ``hse_interp_temp`` : if this is set to ``1``, then we fill the
+  temperatures in the ghost cells via linear extrapolation, using the
+  2 interior zones just inside the domain.  Otherwise, we take the
+  temperature in the ghost cells to be constant.
+
+* ``hse_fixed_temp`` : if this is positive, then we set the
+  temperature in the ghost cells to the value specified.  This
+  requires ``hse_interp_temp = 0``.
+
+
+
+Interface states at reflecting boundary
+=======================================
+
+For all methods, we enforce the reflecting condition on the interface
+states directly by reflecting the state just inside the domain to
+overwrite the state on the reflecting boundary just outside of the
+domain.  This is done for all variables (flipping the sign on the
+normal velocity state).  This is especially important for
+reconstruction that used a one-sided stencil (like the 4th order
+method).
+
+
+Test problems
+=============
+
+``Castro/Exec/gravity_tests/hse_convergence_general`` can be used to
+test the different HSE approaches.  This sets up a 1-d X-ray burst
+atmosphere (based on the ``flame_wave`` setup).  Richardson
+extrapolation can be used to measure the convergence rate (or just
+look at how the peak velocity changes).
+
+

Docs/source/index.rst

@@ -36,6 +36,7 @@ https://github.com/amrex-astro/Castro
    build_system
    debugging
    Hydrodynamics
+   hse
    mhd
    gravity
    diffusion

Docs/source/problem_setups.rst

@@ -77,7 +77,7 @@ problems are:
 
    * ``Sod``: A one-dimensional shock tube setup, including the
      classic Sod problem. This setup was used in the original Castro
-     paper.
+     paper :cite:`castro_I`.
 
    * ``Sod_stellar``: A version of the Sod shock tube for the general
      stellar equation of state. This setup and the included inputs
@@ -147,7 +147,7 @@ problems are:
      Maestro reaction paper :cite:`maestro:III`.
 
    * ``reacting_convergence``: a simple reacting hydrodynamics problem for measuring convergence,
-     used in :cite:`castro-sdc`.
+     used in :cite:`castro-sdc` and :cite:`strang_rnaas`
 
 * ``science``:
 
@@ -161,7 +161,7 @@ problems are:
      was used for the testing in :cite:`eiden:2020`.
 
    * ``flame_wave``: this is a model of a flame propagating across a neutron star as a model for
-     an X-ray burst.  This was presented in :cite:`eiden:2020`.
+     an X-ray burst.  This was presented in :cite:`eiden:2020` and :cite:`harpole2021dynamics`.
 
    * ``nova``: this models convection at the base of an accreted layer
      on a white dwarf as a model of a nova.

Docs/source/refs.bib

@@ -1006,3 +1006,60 @@ @article{amrex-joss
   title = {AMReX: a framework for block-structured adaptive mesh refinement},
   journal = {Journal of Open Source Software}
 }
+
+
+@article{saltzman:1994,
+	title = {An {Unsplit} {3D} {Upwind} {Method} for {Hyperbolic} {Conservation} {Laws}},
+	volume = {115},
+	journal = {Journal of Computational Physics},
+	author = {Saltzman, Jeff},
+	year = {1994},
+	pages = {153--168},
+}
+
+
+@ARTICLE{ppm-hse,
+       author = {{Zingale}, M. and {Dursi}, L.~J. and {ZuHone}, J. and {Calder}, A.~C. and {Fryxell}, B. and {Plewa}, T. and {Truran}, J.~W. and {Caceres}, A. and {Olson}, K. and {Ricker}, P.~M. and {Riley}, K. and {Rosner}, R. and {Siegel}, A. and {Timmes}, F.~X. and {Vladimirova}, N.},
+        title = "{Mapping Initial Hydrostatic Models in Godunov Codes}",
+      journal = {\apjs},
+     keywords = {Hydrodynamics, Methods: Numerical, Stellar Dynamics, Astrophysics},
+         year = 2002,
+        month = dec,
+       volume = {143},
+       number = {2},
+        pages = {539-565},
+          doi = {10.1086/342754},
+archivePrefix = {arXiv},
+       eprint = {astro-ph/0208031},
+ primaryClass = {astro-ph},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2002ApJS..143..539Z},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@article{strang_rnaas,
+	doi = {10.3847/2515-5172/abf3cb},
+	url = {https://doi.org/10.3847/2515-5172/abf3cb},
+	year = 2021,
+	month = {apr},
+	publisher = {American Astronomical Society},
+	volume = {5},
+	number = {4},
+	pages = {71},
+	author = {M. Zingale and M. P. Katz and D. E. Willcox and A. Harpole},
+	title = {Practical Effects of Integrating Temperature with Strang Split Reactions},
+	journal = {Research Notes of the {AAS}},
+	abstract = {Many astrophysical environments involve convective or explosive flows driven by thermonuclear reactions (Type Ia supernovae, classical novae, X-ray bursts, stellar evolution). Simulation codes need to accurately capture the interactions between reactions and hydrodynamics to produce realistic models of these events. For astrophysical reacting flows, operator splitting is commonly used to couple hydrodynamics and reactions. Each process operates independent of one another, but by staggering the updates in a symmetric fashion (via Strang splitting) second order accuracy in time can be achieved. However, approximations are often made to the reacting system, including the choice of whether or not to integrate temperature with the species. Here we demonstrate through a simple convergence test that integrating an energy equation together with reactions achieves the best convergence when modeling reactive flows with Strang splitting. Additionally, second order convergence cannot be achieved without integrating an energy or temperature equation.},
+  subject = {algorithms: compressible hydro}
+}
+
+@article{harpole2021dynamics,
+      title={Dynamics of Laterally Propagating Flames in X-ray Bursts. II. Realistic Burning & Rotation}, 
+      author={A. Harpole and N. M. Ford and K. Eiden and M. Zingale and D. E. Willcox and Y. Cavecchi and M. P. Katz},
+      journal = {\apj},
+      volume = {912},
+      number = {1},
+      pages = {36},
+      year = 2021,
+      subject = {X-ray bursts},
+      url = {https://ui.adsabs.harvard.edu/abs/2021arXiv210200051H/abstract}
+}

Docs/source/rotation.rst

@@ -68,9 +68,6 @@ The main parameters that affect rotation are:
    include the forcing from the time derivative of the rotation
    frequency (default: 1)
 
--  ``castro.state_in_rotating_frame`` : whether state
-   variables are measured in the rotating frame (default: 1)
-
 -  ``castro.rot_source_type`` : method of updating the
    energy during a rotation update (default: 4)