Merge branch 'development' into joss_paper

.github/workflows/cuda.yml

@@ -26,7 +26,7 @@ jobs:
           cd ../..
 
       - name: Dependencies
-        run: .github/workflows/dependencies/dependencies_nvcc11.sh
+        run: .github/workflows/dependencies/dependencies_nvcc.sh 11.7
 
       - name: compile test_react (aprox13)
         run: |

.github/workflows/dependencies/dependencies_hip.sh

@@ -40,6 +40,7 @@ echo 'export PATH=/opt/rocm/llvm/bin:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/r
 
 # we should not need to export HIP_PATH=/opt/rocm/hip with those installs
 
+sudo apt-get clean
 sudo apt-get update
 
 # Ref.: https://rocmdocs.amd.com/en/latest/Installation_Guide/Installation-Guide.html#installing-development-packages-for-cross-compilation
@@ -56,7 +57,9 @@ sudo apt-get install -y --no-install-recommends \
     roctracer-dev   \
     rocprofiler-dev \
     rocrand-dev     \
-    rocprim-dev
+    rocfft-dev      \
+    rocprim-dev     \
+    rocsparse-dev
 
 # hiprand-dev is a new package that does not exist in old versions
 sudo apt-get install -y --no-install-recommends hiprand-dev || true

.github/workflows/dependencies/dependencies_nvcc.sh

@@ -0,0 +1,41 @@
+#!/usr/bin/env bash
+#
+# Copyright 2020-2022 Axel Huebl
+#
+# License: BSD-3-Clause-LBNL
+
+set -eu -o pipefail
+
+# `man apt.conf`:
+#   Number of retries to perform. If this is non-zero APT will retry
+#   failed files the given number of times.
+echo 'Acquire::Retries "3";' | sudo tee /etc/apt/apt.conf.d/80-retries
+
+sudo apt-get -qqq update
+sudo apt-get install -y \
+    build-essential     \
+    ca-certificates     \
+    cmake               \
+    g++                 \
+    gfortran            \
+    gnupg               \
+    libopenmpi-dev      \
+    openmpi-bin         \
+    pkg-config          \
+    wget
+
+VERSION_DOTTED=${1-12.0} && VERSION_DASHED=$(sed 's/\./-/' <<< $VERSION_DOTTED)
+curl -O https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/cuda-keyring_1.0-1_all.deb
+sudo dpkg -i cuda-keyring_1.0-1_all.deb
+sudo apt-get update
+sudo apt-get install -y \
+    cuda-command-line-tools-$VERSION_DASHED \
+    cuda-compiler-$VERSION_DASHED           \
+    cuda-cupti-dev-$VERSION_DASHED          \
+    cuda-minimal-build-$VERSION_DASHED      \
+    cuda-nvml-dev-$VERSION_DASHED           \
+    cuda-nvtx-$VERSION_DASHED               \
+    libcufft-dev-$VERSION_DASHED            \
+    libcurand-dev-$VERSION_DASHED           \
+    libcusparse-dev-$VERSION_DASHED
+sudo ln -s cuda-$VERSION_DOTTED /usr/local/cuda

.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

.gitignore

@@ -78,7 +78,7 @@ extern.f90
 # sphinx
 build/
 doxy_files
-sphinx_docs/source/runtime_parameters.rst
+Docs/source/runtime_parameters.rst
 
 
 # C++ parameters

CHANGES.md

@@ -1,3 +1,11 @@
+# 25.01
+
+  * update HIP/CUDA dependences to include sparse libraries (#1686)
+
+  * rename `Opacity_dir` -> `OPACITY_DIR` (#1679)
+
+  * update the integration and NSE docs (#1682)
+
 # 24.12
 
   * documentation improvements (#1661, #1667, #1670)

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
@@ -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.
 

sphinx_docs/source/index.rst → 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

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
@@ -105,30 +105,50 @@ passed into the integration routines.  For this reason, we often need
 to pass both the specific integrator's type (e.g. ``dvode_t``) and
 ``burn_t`` objects into the lower-level network routines.
 
-The overall flow of the integrator is (using VODE as the example):
+Below we outline the overall flow of the integrator (using VODE as the
+example).  Most of the setup and cleanup after calling the particular
+integration routine is the same for all integrators, and is handled by
+the functions ``integrator_setup()`` and ``integrator_cleanup()``.
 
-#. Call the EOS directly on the input ``burn_t`` state using :math:`\rho` and :math:`T` as inputs.
+.. index:: integrator.scale_system, burn_to_integrator, integrator_to_burn
+.. index:: integrator.call_eos_in_rhs, integrator.subtract_internal_energy, integrator.burner_verbose
+
+#. Call the EOS directly on the input ``burn_t`` state using
+   :math:`\rho` and :math:`T` as inputs.
+
+#. Scale the absolute energy tolerance if we are using
+   ``integrator.scale_system``
 
 #. Fill the integrator type by calling ``burn_to_integrator()`` to create a
    ``dvode_t``.
 
-#. call the ODE integrator, ``dvode()``, passing in the ``dvode_t`` _and_ the
+#. Save the initial thermodynamic state for diagnostics and optionally
+   subtracting off the initial energy later.
+
+#. Call the ODE integrator, ``dvode()``, passing in the ``dvode_t`` *and* the
    ``burn_t`` --- as noted above, the auxiliary information that is
    not part of the integration state will be obtained from the
    ``burn_t``.
 
-#. subtract off the energy offset---we now store just the energy released
-   in the ``dvode_t`` integration state.
+#. Convert back to a ``burn_t`` by calling ``integrator_to_burn``
 
-#. convert back to a ``burn_t`` by calling ``integrator_to_burn``
+#. Recompute the temperature if we are using ``integrator.call_eos_in_rhs``.
 
-#. normalize the abundances so they sum to 1.
+#. If we set ``integrator.subtract_internal_energy``, then subtract
+   off the energy offset, the energy stored is now just that generated
+   by reactions.
+
+#. Normalize the abundances so they sum to 1 (except if ``integrator.use_number_density`` is set).
+
+#. Output statistics on the integration if we set ``integrator.burner_verbose``.
+   This is not recommended for big simulations, as it will output information
+   for every zone's burn.
 
 .. index:: integrator.subtract_internal_energy
 
-.. note::
+.. important::
 
-   Upon exit, ``burn_t burn_state.e`` is the energy *released* during
+   By default, upon exit, ``burn_t burn_state.e`` is the energy *released* during
    the burn, and not the actual internal energy of the state.
 
    Optionally, by setting ``integrator.subtract_internal_energy=0``
@@ -155,7 +175,8 @@ The righthand side of the network is implemented by
 
 .. code-block:: c++
 
-   void actual_rhs(burn_t& state, Array1D<Real, 1, neqs>& ydot)
+   AMREX_GPU_HOST_DEVICE AMREX_INLINE
+   void actual_rhs(burn_t& state, amrex::Array1D<amrex::Real, 1, neqs>& ydot)
 
 All of the necessary integration data comes in through state, as:
 
@@ -245,7 +266,11 @@ The analytic Jacobian is specific to each network and is provided by
 
 .. code-block:: c++
 
-   void actual_jac(burn_t& state, MathArray2D<1, neqs, 1, neqs>& jac)
+   template<class MatrixType>
+   AMREX_GPU_HOST_DEVICE AMREX_INLINE
+   void actual_jac(const burn_t& state, MatrixType& jac)
+
+where the ``MatrixType`` is most commonly ``MathArray2D<1, neqs, 1, neqs>``
 
 The Jacobian matrix elements are stored in ``jac`` as:
 
@@ -316,13 +341,9 @@ Thermodynamics and :math:`e` Evolution
 ======================================
 
 The thermodynamic equation in our system is the evolution of the internal energy,
-:math:`e`.
-
-.. note::
-
-   When the system is integrated in an operator-split approach, the
-   energy equation accounts for only the nuclear energy release and
-   not pdV work.
+:math:`e`.  During the course of the integration, we ensure that the temperature stay
+below the value ``integrator.MAX_TEMP`` (defaulting to ``1.e11``) by clamping the
+temperature if necessary.
 
 At initialization, :math:`e` is set to the value from the EOS consistent
 with the initial temperature, density, and composition:
@@ -331,28 +352,40 @@ with the initial temperature, density, and composition:
 
    e_0 = e(\rho_0, T_0, {X_k}_0)
 
-In the integration routines, this is termed the *energy offset*.
-
 As the system is integrated, :math:`e` is updated to account for the
-nuclear energy release,
+nuclear energy release (and thermal neutrino losses),
+
+.. math:: e(t) = e_0 + \int_{t_0}^t \epsilon(\dot{Y}_k) dt
+
+.. note::
 
-.. math:: e(t) = e_0 + \int_{t_0}^t f(\dot{Y}_k) dt
+   When the system is integrated in an operator-split approach, the
+   energy equation accounts for only the nuclear energy release and
+   not pdV work.
 
-As noted above, upon exit, we subtract off this initial offset, so ``state.e`` in
-the returned ``burn_t`` type from the ``actual_integrator``
-call represents the energy *release* during the burn.
+If ``integrator.subtract_internal_energy`` is set, then, on exit, we
+subtract off this initial $e_0$, so ``state.e`` in the returned
+``burn_t`` type from the ``actual_integrator`` call represents the
+energy *release* during the burn.
 
-Integration of Equation :eq:`eq:enuc_integrate`
-requires an evaluation of the temperature at each integration step
-(since the RHS for the species is given in terms of :math:`T`, not :math:`e`).
-This involves an EOS call and is the default behavior of the integration.
-Note also that for the Jacobian, we need the specific heat, :math:`c_v`, since we
-usually calculate derivatives with respect to temperature (as this is the form
-the rates are commonly provided in).
+Integration of Equation :eq:`eq:enuc_integrate` requires an evaluation
+of the temperature at each integration step (since the RHS for the
+species is given in terms of :math:`T`, not :math:`e`).  This involves
+an EOS call and is the default behavior of the integration.
+
+Note also that for the Jacobian, we need the specific heat,
+:math:`c_v`, since we usually calculate derivatives with respect to
+temperature (as this is the form the rates are commonly provided in).
 
 .. index:: integrator.call_eos_in_rhs
 
 .. note::
 
-   If desired, the EOS call can be skipped and the temperature and $c_v$ kept
-   frozen over the entire time interval of the integration by setting ``integrator.call_eos_in_rhs=0``.
+   If desired, the EOS call can be skipped and the temperature and
+   $c_v$ kept frozen over the entire time interval of the integration
+   by setting ``integrator.call_eos_in_rhs=0``.
+
+.. index:: integrator.integrate_energy
+
+We also provide the option to completely remove the energy equation from
+the system by setting ``integrator.integrate_energy=0``.

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