Merge branch 'development' into update_unit_test_docs

Docs/source/index.rst

@@ -60,6 +60,7 @@ of state routines.
    networks
    templated_networks
    screening
+   neutrinos
 
 .. toctree::
    :maxdepth: 1
@@ -71,6 +72,13 @@ of state routines.
    nse
    sdc
 
+.. toctree::
+   :maxdepth: 1
+   :caption: Util / external libraries
+   :hidden:
+
+   util
+
 .. toctree::
    :maxdepth: 1
    :caption: Unit tests

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::
 

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

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}$

Docs/source/refs.bib

@@ -688,3 +688,16 @@ @article{langanke:2001
 abstract = {The weak interaction rates in stellar environments are computed for pf-shell nuclei in the mass range A=45–65 using large-scale shell-model calculations. The calculated capture and decay rates take into consideration the latest experimental energy levels and log ft-values. The rates are tabulated at the same grid points of density and temperature as those used by Fuller, Fowler, and Newman for densities ρY e =10–1011 g/cm3 and temperatures T=107–1011 K, and hence are relevant for both types of supernovae (Type Ia and Type II). Effective 〈ft〉 values for capture rates and average neutrino (antineutrino) energies are also given to facilitate the use of interpolated rates in stellar evolution codes.}
 }
 
+@ARTICLE{itoh:1996,
+       author = {{Itoh}, Naoki and {Hayashi}, Hiroshi and {Nishikawa}, Akinori and {Kohyama}, Yasuharu},
+        title = "{Neutrino Energy Loss in Stellar Interiors. VII. Pair, Photo-, Plasma, Bremsstrahlung, and Recombination Neutrino Processes}",
+      journal = {\apjs},
+     keywords = {DENSE MATTER, ELEMENTARY PARTICLES, RADIATION MECHANISMS: NONTHERMAL, STARS: INTERIORS, METHODS: NUMERICAL},
+         year = 1996,
+        month = feb,
+       volume = {102},
+        pages = {411},
+          doi = {10.1086/192264},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/1996ApJS..102..411I},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}

Docs/source/sdc.rst

@@ -147,7 +147,7 @@ The reactions don’t modify the total density, :math:`\rho`, or momentum,
       \Adv{\rho X_k}^{n+1/2} \\ \Adv{\rho e}^{n+1/2} \\
    \end{array} \right ) +
    \left (
-      \begin{array}{c} \rho \omegadot_k \\ \rho \Sdot \end{array}
+      \begin{array}{c} \rho \omegadot_k \\ \rho \epsilon \end{array}
    \right )
 
 Here the advective terms are piecewise-constant (in time)

Docs/source/util.rst

@@ -0,0 +1,46 @@
+******************************
+Helper Functions and Libraries
+******************************
+
+The ``util/`` directory contains a number of external libraries and simple
+utilities that are used by the different components of Microphysics.
+
+* ``approx_math/`` : these are a set of headers that implement
+  approximations to ``atan()``, ``exp()``, ``log()``, and ``pow()``.
+  These can be much faster than the C++ library versions, especially
+  on GPUs.
+
+* ``autodiff/`` : this is a clone of the C++ autodiff library from
+  https://github.com/autodiff/autodiff
+
+  The header ``microphysics_autodiff.H`` provides a set of interfaces
+  for working with the AMReX datatypes and interfacing with the
+  autodiff library.
+
+* ``build_scripts/`` : a set of python scripts used during the build
+  process to parse the runtime parameters.
+
+* ``cj_detonation/`` : a simple routine to compute the Chapman-Jouguet
+  detonation speed for one of our networks.
+
+* ``esum.H`` : an implementation of the exact sum algorithm based on the
+  msum algorithm by Raymond Hettinger.  It is generated automatically
+  by the ``esum_cxx.py`` script and creates implementations for exact
+  numbers of terms (``esum3()``, ``esum4()``, ...)
+
+* ``gcem/`` : a templated math library that provides implementations of
+  the standard library math functions that can be used in ``constexpr``
+  expressions.  This is from https://github.com/kthohr/gcem
+
+  Some of the constants are redefined in 64-bit floating point in
+  ``microphysics_math.H`` to avoid ``long double`` issues on some
+  architectures.
+
+* ``hybrj/`` : a C++ port of the MINPACK hybrid Powell minimization function
+  to zero a set of functions.
+
+* ``linpack.H`` : a C++ port of the LINPACK ``dgesl`` and ``dgefa`` LU
+  decomposition Gaussian elimination routines.
+
+* ``microphysics_sort.H`` : a set of sorting routines for
+  ``amrex::Array1D`` data.

EOS/breakout/actual_eos.H

@@ -21,6 +21,10 @@ void actual_eos_init ()
         gamma_const = 5.0_rt / 3.0_rt;
     }
 
+    // this EOS assumes that has a 1/mu entry -- make sure that is valid
+    if (aux_names_cxx[AuxZero::iinvmu] != "invmu") {
+        amrex::Error("invalid aux state for breakout EOS");
+    }
 }
 
 template <typename I>
@@ -52,7 +56,7 @@ void actual_eos (I input, T& state)
 
     // Calculate mu. This is the only difference between
     // this EOS and gamma_law.
-    state.mu = 1.0_rt / state.aux[1];
+    state.mu = 1.0_rt / state.aux[AuxZero::iinvmu];
 
     switch (input) {