Merge branch 'development' into doxygen_fixes

zingale · Jan 16, 2025 · 75a104f · 75a104f
2 parents b168435 + aeef432
commit 75a104f
Show file tree

Hide file tree

Showing 259 changed files with 19,307 additions and 4,250 deletions.
diff --git a/.github/workflows/docs-test.yml b/.github/workflows/docs-test.yml
@@ -35,10 +35,10 @@ jobs:
 
       - name: Build docs
         run: |
-          cd sphinx_docs/
+          cd Docs/
           make SPHINXOPTS='-v -W --keep-going' html
 
       - name: Link check
         run: |
-          cd sphinx_docs/
+          cd Docs/
           make SPHINXOPTS=-v linkcheck
diff --git a/.gitignore b/.gitignore
@@ -78,7 +78,7 @@ extern.f90
 # sphinx
 build/
 doxy_files
-sphinx_docs/source/runtime_parameters.rst
+Docs/source/runtime_parameters.rst
 
 
 # C++ parameters

diff --git a/.mailmap b/.mailmap
@@ -9,4 +9,9 @@ Abigail Bishop <[email protected]> abigailbishop <ms.abigailbishop@gmai
 Sydney Andrews <[email protected]> sandshrews <[email protected]>
 Xinlong Li <[email protected]> Xinlong <[email protected]>
 Xinlong Li <[email protected]> sailoridy <[email protected]>
-Xinlong Li <[email protected]> sailor <[email protected]>
+Xinlong Li <[email protected]> sailor <[email protected]>
+Xinlong Li <[email protected]>
+Eric T. Johnson <[email protected]>
+Zhi Chen <[email protected]>
+Zhi Chen <[email protected]>
+Doreen Fan <[email protected]>
diff --git a/sphinx_docs/Doxyfile → Docs/Doxyfile b/sphinx_docs/Doxyfile → Docs/Doxyfile
diff --git a/sphinx_docs/Makefile → Docs/Makefile b/sphinx_docs/Makefile → Docs/Makefile
diff --git a/sphinx_docs/converting_docs.txt → Docs/converting_docs.txt b/sphinx_docs/converting_docs.txt → Docs/converting_docs.txt
diff --git a/sphinx_docs/notes.txt → Docs/notes.txt b/sphinx_docs/notes.txt → Docs/notes.txt
diff --git a/sphinx_docs/pandoc.sh → Docs/pandoc.sh b/sphinx_docs/pandoc.sh → Docs/pandoc.sh
diff --git a/sphinx_docs/rp.py → Docs/rp.py b/sphinx_docs/rp.py → Docs/rp.py
@@ -46,6 +46,8 @@ def pretty_category(path):
         return f"CONDUCTIVITY_DIR={subdir}"
     if path.startswith("integration/"):
         return f"INTEGRATOR_DIR={subdir}"
+    if path.startswith("opacity/"):
+        return f"OPACITY_DIR={subdir}"
     return path
 
 
@@ -65,9 +67,9 @@ def make_rest_table(param_files):
 
         # open the file
         try:
-            f = open(pf, "r")
-        except IOError:
-            sys.exit("ERROR: {} does not exist".format(pf))
+            f = open(pf)
+        except OSError:
+            sys.exit(f"ERROR: {pf} does not exist")
 
         descr = r""
 

diff --git a/sphinx_docs/source/CNO_He_burn.png → Docs/source/CNO_He_burn.png b/sphinx_docs/source/CNO_He_burn.png → Docs/source/CNO_He_burn.png
diff --git a/sphinx_docs/source/ECSN.png → Docs/source/ECSN.png b/sphinx_docs/source/ECSN.png → Docs/source/ECSN.png
diff --git a/...x_docs/source/_static/theme_overrides.css → Docs/source/_static/theme_overrides.css b/...x_docs/source/_static/theme_overrides.css → Docs/source/_static/theme_overrides.css
diff --git a/...nx_docs/source/_templates/dev_layout.html → Docs/source/_templates/dev_layout.html b/...nx_docs/source/_templates/dev_layout.html → Docs/source/_templates/dev_layout.html
diff --git a/sphinx_docs/source/_templates/layout.html → Docs/source/_templates/layout.html b/sphinx_docs/source/_templates/layout.html → Docs/source/_templates/layout.html
diff --git a/sphinx_docs/source/autodiff.rst → Docs/source/autodiff.rst b/sphinx_docs/source/autodiff.rst → Docs/source/autodiff.rst
diff --git a/sphinx_docs/source/cno_extras_hide_alpha.png → Docs/source/cno_extras_hide_alpha.png b/sphinx_docs/source/cno_extras_hide_alpha.png → Docs/source/cno_extras_hide_alpha.png
diff --git a/sphinx_docs/source/comprehensive_tests.rst → Docs/source/comprehensive_tests.rst b/sphinx_docs/source/comprehensive_tests.rst → Docs/source/comprehensive_tests.rst
diff --git a/sphinx_docs/source/conf.py → Docs/source/conf.py b/sphinx_docs/source/conf.py → Docs/source/conf.py
@@ -68,6 +68,7 @@ def get_version():
 
 # bibtex
 bibtex_bibfiles = ["refs.bib"]
+bibtex_default_style = "unsrt"
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']

diff --git a/sphinx_docs/source/data_structures.rst → Docs/source/data_structures.rst b/sphinx_docs/source/data_structures.rst → Docs/source/data_structures.rst
diff --git a/sphinx_docs/source/design.rst → Docs/source/design.rst b/sphinx_docs/source/design.rst → Docs/source/design.rst
@@ -12,6 +12,8 @@ and the generic solvers:
 
 * ``constants``: fundamental constants
 
+* ``Docs``: the sphinx source for this documentation
+
 * ``EOS/``: the various equations of state
 
 * ``integration/``: the ODE integration routines used for the
@@ -24,7 +26,7 @@ and the generic solvers:
   righthand side of the network, as the actual integrators are decoupled from
   the network.
 
-* ``neutrinos/``: neutino loss source terms for the network energy equation.
+* ``neutrinos/``: neutrino loss source terms for the network energy equation.
 
 * ``opacity/``: opacity routines for radiation transport
 
@@ -34,8 +36,6 @@ and the generic solvers:
 * ``screening/``: common electron screening factors used by some of the
   reaction networks.
 
-* ``sphinx_docs``: the sphinx source for this documentation
-
 * ``unit_test/``: self-contained unit tests for Microphysics. These don’t
   need any application code to build, but will require AMReX.
 

diff --git a/sphinx_docs/source/doxygen_network.png → Docs/source/doxygen_network.png b/sphinx_docs/source/doxygen_network.png → Docs/source/doxygen_network.png
diff --git a/sphinx_docs/source/eos.rst → Docs/source/eos.rst b/sphinx_docs/source/eos.rst → Docs/source/eos.rst
diff --git a/sphinx_docs/source/getting_started.rst → Docs/source/getting_started.rst b/sphinx_docs/source/getting_started.rst → Docs/source/getting_started.rst
@@ -27,7 +27,7 @@ to set the ``AMREX_HOME`` environment variable to point to the
 
 .. index:: burn_cell
 
-A good unit test to start with is ``burn_cell`` -- this is simply a
+A good unit test to start with is ``burn_cell`` --- this is simply a
 one-zone burn.  In ``Microphysics/`` do:
 
 .. prompt:: bash

diff --git a/sphinx_docs/source/index.rst → Docs/source/index.rst b/sphinx_docs/source/index.rst → Docs/source/index.rst
@@ -10,10 +10,11 @@ AMReX-Astro Microphysics is a collection of microphysics routines (equations of
 reaction networks, ...) and utilities (ODE integrators, NSE solvers)
 for astrophysical simulation codes.
 
-The original design was to support the `AMReX
-<https://github.com/amrex-codes/amrex>`_ codes `CASTRO
-<https://github.com/amrex-astro/Castro>`_ and MAESTRO (now `MAESTROeX
-<https://github.com/amrex-astro/MAESTROeX>`_). These all have a
+The original design was to support codes based on the `AMReX
+<https://github.com/amrex-codes/amrex>`_ adaptive mesh refinement library :cite:`amrex_joss`,
+including `CASTRO
+<https://github.com/amrex-astro/Castro>`_ :cite:`castro_I`, `MAESTROeX
+<https://github.com/amrex-astro/MAESTROeX>`_ :cite:`maestroex`, and Quokka :cite:`quokka`. These all have a
 consistent interface and the separate Microphysics repository allows
 them to share the same equation of state, reaction networks, and more.
 Later, Microphysics was adopted by the `Quokka <https://github.com/quokka-astro/quokka>`_
@@ -60,6 +61,7 @@ of state routines.
    networks
    templated_networks
    screening
+   neutrinos
 
 .. toctree::
    :maxdepth: 1
@@ -71,6 +73,13 @@ of state routines.
    nse
    sdc
 
+.. toctree::
+   :maxdepth: 1
+   :caption: Util / external libraries
+   :hidden:
+
+   util
+
 .. toctree::
    :maxdepth: 1
    :caption: Unit tests

diff --git a/sphinx_docs/source/integrators.rst → Docs/source/integrators.rst b/sphinx_docs/source/integrators.rst → Docs/source/integrators.rst
@@ -14,13 +14,13 @@ The equations we integrate to do a nuclear burn are:
    :label: eq:spec_integrate
 
 .. math::
-   \frac{de}{dt} = f(\rho,X_k,T)
+   \frac{de}{dt} = \epsilon(\rho,X_k,T)
    :label: eq:enuc_integrate
 
 Here, :math:`X_k` is the mass fraction of species :math:`k`, :math:`e`
 is the specific nuclear energy created through reactions. Also needed
 are density :math:`\rho`, temperature :math:`T`, and the specific
-heat. The function :math:`f` provides the energy release from
+heat. The function :math:`\epsilon` provides the energy release from
 reactions and can often be expressed in terms of the instantaneous
 reaction terms, :math:`\dot{X}_k`. As noted in the previous section,
 this is implemented in a network-specific manner.
@@ -46,7 +46,7 @@ energy. This allows us to easily call the EOS during the burn to obtain the temp
       :label: eq:spec_n_integrate
 
    .. math::
-      \frac{de}{dt} = f(\rho,n_k,T)
+      \frac{de}{dt} = \epsilon(\rho,n_k,T)
       :label: eq:enuc_n_integrate
 
    The effect of this flag in the integrators is that we don't worry
@@ -355,7 +355,7 @@ with the initial temperature, density, and composition:
 As the system is integrated, :math:`e` is updated to account for the
 nuclear energy release (and thermal neutrino losses),
 
-.. math:: e(t) = e_0 + \int_{t_0}^t f(\dot{Y}_k) dt
+.. math:: e(t) = e_0 + \int_{t_0}^t \epsilon(\dot{Y}_k) dt
 
 .. note::
 

diff --git a/sphinx_docs/source/mathsymbols.tex → Docs/source/mathsymbols.tex b/sphinx_docs/source/mathsymbols.tex → Docs/source/mathsymbols.tex
diff --git a/sphinx_docs/source/networks-overview.rst → Docs/source/networks-overview.rst b/sphinx_docs/source/networks-overview.rst → Docs/source/networks-overview.rst
diff --git a/sphinx_docs/source/networks.rst → Docs/source/networks.rst b/sphinx_docs/source/networks.rst → Docs/source/networks.rst
@@ -26,10 +26,10 @@ is stored as ``mion(:)`` in the network.
 
 The energy release per gram is converted from the rates as:
 
-.. math:: \edot = -N_A c^2 \sum_k \frac{dY_k}{dt} M_k - \edotnu
+.. math:: \epsilon = -N_A c^2 \sum_k \frac{dY_k}{dt} M_k - \epsilon_\nu
 
 where :math:`N_A` is Avogadro’s number (to convert this to “per gram”)
-and :math:`\edotnu` is the neutrino loss term.
+and :math:`\edotnu` is the neutrino loss term (see :ref:`neutrino_loss`).
 
 
 ``general_null``
@@ -137,7 +137,7 @@ network is interpolated based on the density. It is stored as the
 binding energy (ergs/g) *per nucleon*, with a sign convention that
 binding energies are negative. The energy generation rate is then:
 
-.. math:: \edot = q \frac{dX(\isotm{C}{12})}{dt} = q A_{\isotm{C}{12}} \frac{dY(\isotm{C}{12})}{dt}
+.. math:: \epsilon = q \frac{dX(\isotm{C}{12})}{dt} = q A_{\isotm{C}{12}} \frac{dY(\isotm{C}{12})}{dt}
 
 (this is positive since both :math:`q` and :math:`dY/dt` are negative)
 
@@ -244,7 +244,7 @@ Finally, for the
 energy generation, we take our reaction to release a specific energy,
 :math:`[\mathrm{erg~g^{-1}}]`, of :math:`\qburn`, and our energy source is
 
-.. math:: \edot = -\qburn \frac{dX_f}{dt}
+.. math:: \epsilon = -\qburn \frac{dX_f}{dt}
 
 There are a number of parameters we use to control the constants in
 this network. This is one of the few networks that was designed

diff --git a/Docs/source/neutrinos.rst b/Docs/source/neutrinos.rst
@@ -0,0 +1,47 @@
+.. _neutrino_loss:
+
+***************
+Neutrino Losses
+***************
+
+We model thermal neutrino losses (plasma, photo-, pair-,
+recombination, and Bremsstrahlung) using the method described in
+:cite:`itoh:1996`.  This neutrino loss term is included in all of the
+reaction networks by default, and modifies the energy equation to have
+the form (for Strang splitting):
+
+.. math::
+
+   \frac{de}{dt} = \epsilon - \epsilon_\nu
+
+where $\epsilon_\nu$ are the thermal neutrino losses.
+
+.. note::
+
+   Thermal neutrino losses can be disabled at compile time by setting the make
+   variable ``USE_NEUTRINOS = FALSE``.
+
+The interface for the neutrino loss function is:
+
+.. code:: c++
+
+   template <int do_derivatives>
+   AMREX_GPU_HOST_DEVICE AMREX_INLINE
+   void sneut5(const amrex::Real temp, const amrex::Real den,
+               const amrex::Real abar, const amrex::Real zbar,
+               amrex::Real& snu, amrex::Real& dsnudt, amrex::Real& dsnudd,
+               amrex::Real& dsnuda, amrex::Real& dsnudz)
+
+Here, the template parameter, ``do_derivatives``, can be used to disable the code
+the computes the derivatives of the neutrino loss, for example, if a numerical Jacobian
+is used.  The output is
+
+* ``snu`` : $\epsilon_\nu$, the neutrino loss in erg/g/s
+
+* ``dsnudt`` : $d\epsilon_\nu/dT$
+
+* ``dsnudd`` : $d\epsilon_\nu/d\rho$
+
+* ``dsnuda`` : $d\epsilon_\nu/d\bar{A}$
+
+* ``dsnudz`` : $d\epsilon_\nu/d\bar{Z}$
diff --git a/sphinx_docs/source/nova.png → Docs/source/nova.png b/sphinx_docs/source/nova.png → Docs/source/nova.png
diff --git a/sphinx_docs/source/nova2.png → Docs/source/nova2.png b/sphinx_docs/source/nova2.png → Docs/source/nova2.png
diff --git a/sphinx_docs/source/nse.rst → Docs/source/nse.rst b/sphinx_docs/source/nse.rst → Docs/source/nse.rst
diff --git a/sphinx_docs/source/ode_integrators.rst → Docs/source/ode_integrators.rst b/sphinx_docs/source/ode_integrators.rst → Docs/source/ode_integrators.rst
@@ -146,46 +146,80 @@ Tolerances
 Tolerances dictate how accurate the ODE solver must be while solving
 equations during a simulation.  Typically, the smaller the tolerance
 is, the more accurate the results will be.  However, if the tolerance
-is too small, the code may run for too long or the ODE solver will
-never converge.  In these simulations, ``rtol`` values will set the
-relative tolerances and ``atol`` values will set the absolute tolerances
-for the ODE solver.  Often, one can find and set these values in an
-input file for a simulation.
-
-:numref:`fig:tolerances` shows the results of a simple simulation using the
-burn_cell unit test to determine
-what tolerances are ideal for simulations.
-For this investigation, it was assumed that a run with a tolerance of :math:`10^{-12}`
-corresponded to an exact result,
-so it is used as the basis for the rest of the tests.
-From the figure, one can infer that the :math:`10^{-3}` and :math:`10^{-6}` tolerances
-do not yield the most accurate results
-because their relative error values are fairly large.
-However, the test with a tolerance of :math:`10^{-9}` is accurate
-and not so low that it takes incredible amounts of computer time,
-so :math:`10^{-9}` should be used as the default tolerance in future simulations.
+is too small, the code may run for too long, the ODE solver will
+never converge, or it might require at timestep that underflows.
+
+.. index:: integrator.rtol_spec, integrator.rtol_enuc, integrator.atol_spec, integrator.atol_enuc
+
+There are separate tolerances for the mass fractions and the energy,
+and there are both relative and absolute tolerances which act together.
+The tolerances are:
+
+* ``integrator.rtol_spec`` : the relative tolerance for the species
+  (mass fractions when running with Strang and partial densities when
+  running with SDC).
+
+* ``integrator.rtol_enuc`` : the relative tolerance on the energy
+  (specific internal energy when running with Strang, internal energy
+  density when running with SDC).
+
+* ``integrator.atol_spec`` : the absolute tolerance for the species
+  (this is always interpreted in terms of mass fraction and the appropriate
+  density weighting will be added for SDC).
+
+* ``integrator.atol_enuc`` : the absolute tolerance for energy -- this
+  is generally not interesting, since the energy is so large and therefore
+  best served via a relative tolerance.
+
+The tolerances are combined, e.g. for species, as:
+
+.. math::
+
+   \epsilon_{\mathrm{total}, k} = \epsilon_\mathrm{abs} + \epsilon_\mathrm{rel} |X_k|
+
+so if the mass fraction, $X_k$, is very small, then the absolute tolerance
+will set the error that the integrator tries to achieve.  If the mass fraction
+is large, $\mathcal{O}(1)$, then the relative tolerance dominates.
+
+Some suggestions when setting tolerances:
+
+.. index:: integrator.X_reject_buffer
+
+* If a burn does not converge with one type of Jacobian (analytic or
+  numerical) then it may do better with the other type.  This can be
+  automated via the ``integrator.use_burn_retry`` mechanism described
+  above.
+
+* Sometimes a burn completes better if the absolute tolerances are
+  made even smaller -- this will require the integrator to track trace
+  species better which can help with equilibrium better.
+
+* The VODE integrator has additional logic meant to ensure that
+  species don't change too much per timestep.  This is controlled by
+  ``integrator.X_reject_buffer``.  If a species $k$, has a mass
+  fraction $X_k > \mbox{X_reject_buffer} \cdot \mbox{atol_spec}$ then
+  we reject a VODE timestep if the mass fraction changes by more than
+  a factor of 4 in a single VODE timestep and we try again.  This is
+  all done internally to VODE.  Making ``X_reject_buffer`` larger will
+  allow it to ignore more trace species.
+
+Below is a comparison of how the tolerances affect the nucleosynthesis.
+This is run using ``burn_cell`` and the ``aprox13`` network.  Four separate
+runs were done, using tolerances of $10^{-3}$, $10^{-5}$, $10^{-8}$, and $10^{-12}$
+(all 4 tolerance parameters were set to the same value).  The run with the tightest
+tolerances ($10^{-12}$) was taken as the reference and relative errors were
+computed with respect to it.  The scripts for this are in ``Microphysics/unit_test/burn_cell/compare_tolerances/``.
 
 .. _fig:tolerances:
-.. figure:: tolerances.png
-   :alt: Relative error plot
+.. figure:: tolerance-compare.png
+   :alt: Relative error in mass fractions
    :width: 100%
 
    Relative error of runs with varying tolerances as compared
    to a run with an ODE tolerance of :math:`10^{-12}`.
 
-The integration tolerances for the burn are controlled by
-``rtol_spec`` and  ``rtol_enuc``,
-which are the relative error tolerances for
-:eq:`eq:spec_integrate` and :eq:`eq:enuc_integrate`,
-respectively. There are corresponding
-``atol`` parameters for the absolute error tolerances. Note that
-not all integrators handle error tolerances the same way—see the
-sections below for integrator-specific information.
-
-The absolute error tolerances are set by default
-to :math:`10^{-12}` for the species, and a relative tolerance of :math:`10^{-6}`
-is used for the temperature and energy.
-
+We see that using a tolerance of $10^{-5}$ generally gives reasonable mass
+fractions.  Looser than this can produce large errors.
 
 Controlling Species $\sum_k X_k = 1$
 ====================================
@@ -222,6 +256,7 @@ constraint on the intermediate states during the integration.
   This is enabled by default.
 
 
+
 Retry Mechanism
 ===============