diff --git a/.clang-tidy b/.clang-tidy
index 1e228494af..2bd4310f7a 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -19,6 +19,7 @@ Checks: >
        -misc-include-cleaner,
        -misc-non-private-member-variables-in-classes,
        -misc-use-anonymous-namespace,
+       -misc-use-internal-linkage,
      modernize-*,
        -modernize-avoid-c-arrays,
        -modernize-use-trailing-return-type,
@@ -29,12 +30,14 @@ Checks: >
        -readability-avoid-const-params-in-decls,
        -readability-braces-around-statements,
        -readability-else-after-return,
+       -readability-enum-initial-value,
        -readability-function-cognitive-complexity,
        -readability-function-size,
        -readability-identifier-length,
        -readability-implicit-bool-conversion,
        -readability-isolate-declaration,
        -readability-magic-numbers,
+       -readability-math-missing-parentheses,
        -readability-named-parameter,
        -readability-simplify-boolean-expr,
      mpi-*,
diff --git a/.github/workflows/burn_cell_metal_chem.yml b/.github/workflows/burn_cell_metal_chem.yml
index 6c587eebf2..13ec5f98ea 100644
--- a/.github/workflows/burn_cell_metal_chem.yml
+++ b/.github/workflows/burn_cell_metal_chem.yml
@@ -31,7 +31,7 @@ jobs:
       - name: Compile
         run: |
           cd unit_test/burn_cell_metal_chem
-          make -j 2
+          make -j 4
 
       - name: Run and compare outputs for different Z values, also including cosmic ray ionization
         run: |
diff --git a/.github/workflows/burn_cell_primordial_chem.yml b/.github/workflows/burn_cell_primordial_chem.yml
index 04f04643e2..7b2ba7a781 100644
--- a/.github/workflows/burn_cell_primordial_chem.yml
+++ b/.github/workflows/burn_cell_primordial_chem.yml
@@ -31,7 +31,7 @@ jobs:
       - name: Compile and run
         run: |
           cd unit_test/burn_cell_primordial_chem
-          make -j 2
+          make -j 4
           ./main1d.gnu.DEBUG.ex inputs_primordial_chem amrex.fpe_trap_{invalid,zero,overflow}=1 > test.out
 
       - name: Print backtrace
diff --git a/.github/workflows/c-linter.yml b/.github/workflows/c-linter.yml
index 56d6bc7235..30e4baa8c8 100644
--- a/.github/workflows/c-linter.yml
+++ b/.github/workflows/c-linter.yml
@@ -26,7 +26,7 @@ jobs:
         uses: AMReX-Astro/cpp-linter-action@main
         with:
           build_path: 'unit_test/test_react'
-          make_options: '-j 2 USE_OMP=FALSE USE_MPI=FALSE USE_CUDA=FALSE DEBUG=TRUE'
+          make_options: '-j 4 USE_OMP=FALSE USE_MPI=FALSE USE_CUDA=FALSE DEBUG=TRUE'
           ignore_files: 'amrex|util/gcem'
           header_filter: '(/conductivity/|/constants/|/EOS/|/integration/|/interfaces/|/networks/|/neutrinos/|/nse_solver/|/opacity/|/rates/|/screening/|/util/|^\./).*\.H$'
           config_file: ${GITHUB_WORKSPACE}/.clang-tidy
diff --git a/.github/workflows/castro-development.yml b/.github/workflows/castro-development.yml
index 451f4c2037..273e32a7c6 100644
--- a/.github/workflows/castro-development.yml
+++ b/.github/workflows/castro-development.yml
@@ -47,7 +47,7 @@ jobs:
         export MICROPHYSICS_HOME=${PWD}
 
         cd Castro/Exec/science/flame_wave/
-        make -j2 CCACHE=ccache USE_MPI=FALSE
+        make -j 4 CCACHE=ccache USE_MPI=FALSE
 
         ccache -s
         du -hs ~/.cache/ccache
@@ -64,7 +64,7 @@ jobs:
         export MICROPHYSICS_HOME=${PWD}
 
         cd Castro/Exec/science/subchandra/
-        make -j2 CCACHE=ccache USE_MPI=FALSE
+        make -j 4 CCACHE=ccache USE_MPI=FALSE
 
         ccache -s
         du -hs ~/.cache/ccache
diff --git a/.github/workflows/castro.yml b/.github/workflows/castro.yml
index 0ea9864aca..4a5b1d58e8 100644
--- a/.github/workflows/castro.yml
+++ b/.github/workflows/castro.yml
@@ -51,7 +51,7 @@ jobs:
         export MICROPHYSICS_HOME=${PWD}
 
         cd Castro/Exec/science/flame_wave/
-        make -j2 CCACHE=ccache USE_MPI=FALSE
+        make -j 4 CCACHE=ccache USE_MPI=FALSE
 
         ccache -s
         du -hs ~/.cache/ccache
@@ -68,7 +68,7 @@ jobs:
         export MICROPHYSICS_HOME=${PWD}
 
         cd Castro/Exec/science/subchandra/
-        make -j2 CCACHE=ccache USE_MPI=FALSE
+        make -j 4 CCACHE=ccache USE_MPI=FALSE
 
         ccache -s
         du -hs ~/.cache/ccache
diff --git a/.github/workflows/clang-tidy.yml b/.github/workflows/clang-tidy.yml
index 6482717f89..6563da5f78 100644
--- a/.github/workflows/clang-tidy.yml
+++ b/.github/workflows/clang-tidy.yml
@@ -26,10 +26,10 @@ jobs:
 
       - name: Install dependencies
         run: |
-          .github/workflows/dependencies_clang-tidy-apt-llvm.sh 17
+          .github/workflows/dependencies_clang-tidy-apt-llvm.sh 19
 
       - name: Compile burn_cell_sdc
         run: |
           cd unit_test/burn_cell_sdc
-          make USE_MPI=FALSE USE_CLANG_TIDY=TRUE CLANG_TIDY=clang-tidy-17 CLANG_TIDY_WARN_ERROR=TRUE -j 4
+          make USE_MPI=FALSE USE_CLANG_TIDY=TRUE CLANG_TIDY=clang-tidy-19 CLANG_TIDY_WARN_ERROR=TRUE -j 4
 
diff --git a/.github/workflows/cmake_build_cell_primordial_chem.yml b/.github/workflows/cmake_build_cell_primordial_chem.yml
index 844ab4c8d7..6063770722 100644
--- a/.github/workflows/cmake_build_cell_primordial_chem.yml
+++ b/.github/workflows/cmake_build_cell_primordial_chem.yml
@@ -18,5 +18,5 @@ jobs:
         run: |
           mkdir build && cd build
           cmake .. -DBUILD_UNIT_TEST=true -DBUILD_AMReX=true
-          make -j2
+          make -j 4
           ctest --output-on-failure
diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml
index 7e418249a8..fab16b93fe 100644
--- a/.github/workflows/cuda.yml
+++ b/.github/workflows/cuda.yml
@@ -32,17 +32,17 @@ jobs:
         run: |
           export PATH=/usr/local/nvidia/bin:/usr/local/cuda/bin:${PATH}
           cd unit_test/test_react
-          make NETWORK_DIR=aprox13 USE_CUDA=TRUE COMP=gnu USE_MPI=FALSE -j 2
+          make NETWORK_DIR=aprox13 USE_CUDA=TRUE COMP=gnu USE_MPI=FALSE -j 4
 
       - name: compile test_react (ignition_reaclib/URCA-simple)
         run: |
           export PATH=/usr/local/nvidia/bin:/usr/local/cuda/bin:${PATH}
           cd unit_test/test_react
           make realclean
-          make NETWORK_DIR=ignition_reaclib/URCA-simple USE_CUDA=TRUE COMP=gnu USE_MPI=FALSE -j 2
+          make NETWORK_DIR=ignition_reaclib/URCA-simple USE_CUDA=TRUE COMP=gnu USE_MPI=FALSE -j 4
 
       - name: compile test_nse_net (ase)
         run: |
           export PATH=/usr/local/nvidia/bin:/usr/local/cuda/bin:${PATH}
           cd unit_test/test_nse_net
-          make USE_CUDA=TRUE COMP=gnu USE_MPI=FALSE -j 2
+          make USE_CUDA=TRUE COMP=gnu USE_MPI=FALSE -j 4
diff --git a/.github/workflows/macos_build_cell_metal_chem.yml b/.github/workflows/macos_build_cell_metal_chem.yml
index 3758f43cae..7880dabafc 100644
--- a/.github/workflows/macos_build_cell_metal_chem.yml
+++ b/.github/workflows/macos_build_cell_metal_chem.yml
@@ -17,5 +17,5 @@ jobs:
         run: |
           mkdir build && cd build
           cmake .. -DBUILD_UNIT_TEST_MC=true -DBUILD_AMReX=true
-          make -j2
+          make -j 4
           ctest --output-on-failure
diff --git a/.github/workflows/macos_build_cell_primordial_chem.yml b/.github/workflows/macos_build_cell_primordial_chem.yml
index 3e7698432d..0390c426c1 100644
--- a/.github/workflows/macos_build_cell_primordial_chem.yml
+++ b/.github/workflows/macos_build_cell_primordial_chem.yml
@@ -17,5 +17,5 @@ jobs:
         run: |
           mkdir build && cd build
           cmake .. -DBUILD_UNIT_TEST_PC=true -DBUILD_AMReX=true
-          make -j2
+          make -j 4
           ctest --output-on-failure
diff --git a/.github/workflows/test_neutrinos.yml b/.github/workflows/test_neutrinos.yml
index 46cb99cd8c..4ab1a80cc8 100644
--- a/.github/workflows/test_neutrinos.yml
+++ b/.github/workflows/test_neutrinos.yml
@@ -29,7 +29,7 @@ jobs:
       - name: Build the fextrema tool
         run: |
           cd external/amrex/Tools/Plotfile
-          make programs=fextrema -j 2
+          make programs=fextrema -j 4
 
       - name: Compile
         run: |
diff --git a/.github/workflows/test_partition_functions.yml b/.github/workflows/test_partition_functions.yml
new file mode 100644
index 0000000000..2752920610
--- /dev/null
+++ b/.github/workflows/test_partition_functions.yml
@@ -0,0 +1,44 @@
+name: test_partition_functions
+
+on: [pull_request]
+jobs:
+  test_partition_functions:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+
+      - name: Get AMReX
+        run: |
+          mkdir external
+          cd external
+          git clone https://github.com/AMReX-Codes/amrex.git
+          cd amrex
+          git checkout development
+          echo 'AMREX_HOME=$(GITHUB_WORKSPACE)/external/amrex' >> $GITHUB_ENV
+          echo $AMREX_HOME
+          if [[ -n "${AMREX_HOME}" ]]; then exit 1; fi
+          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
+
+      - name: Compile
+        run: |
+          cd unit_test/test_part_func
+          make clean
+          make -j 4
+
+      - name: Run test_part_func
+        run: |
+          cd unit_test/test_part_func
+          ./main3d.gnu.ex > test.out
+
+      - name: Compare to stored output
+        run: |
+          cd unit_test/test_part_func
+          diff -I "^Initializing AMReX" -I "^AMReX" -I "^reading in reaclib rates" test.out ci-benchmarks/part_func.out
+
diff --git a/.github/workflows/test_rhs.yml b/.github/workflows/test_rhs.yml
index 4654200513..21b4a9dcd8 100644
--- a/.github/workflows/test_rhs.yml
+++ b/.github/workflows/test_rhs.yml
@@ -29,7 +29,7 @@ jobs:
       - name: Build the fextrema tool
         run: |
           cd external/amrex/Tools/Plotfile
-          make programs=fextrema -j 2
+          make programs=fextrema -j 4
 
       - name: Compile, test_rhs (VODE, ignition_simple)
         run: |
diff --git a/CHANGES.md b/CHANGES.md
index ab4a276748..224a101517 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,10 +1,14 @@
-# 24.11
+# 24.12
 
-  * a new metal chemistry network was added (#1648, #1650, #1651)
+  * documentation improvements (#1661, #1667, #1670)
 
-  * add dust temperature to primordial chem (#1649)
+  * optimize tabular NSE EOS calls (#1668)
 
-  * documentation updates (#1652)
+  * CI fixes (#1666, #1671, #1675) and new partition function CI
+    (#1673)
+
+  * `burn_cell` can now initialize all mass fractions to be equal
+    (#1665)
 
 # 24.10
 
diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H
index 055ddc1394..20c6747224 100644
--- a/integration/VODE/vode_dvhin.H
+++ b/integration/VODE/vode_dvhin.H
@@ -122,17 +122,13 @@ void dvhin (BurnT& state, DvodeT& vstate, amrex::Real& H0, int& NITER, int& IER)
 
         }
 
-        // Iteration done. Apply bounds, bias factor, and sign. Then exit.
+        // Iteration done. Apply bounds and bias factor. Then exit.
         H0 = hnew * 0.5_rt;
-        if (H0 < HLB) {
-            H0 = HLB;
-        }
-        if (H0 > HUB) {
-            H0 = HUB;
-        }
+        H0 = std::clamp(H0, HLB, HUB);
 
     }
 
+    // apply sign
     H0 = std::copysign(H0, vstate.tout - vstate.t);
     NITER = iter;
     IER = 0;
diff --git a/integration/integrator_setup_strang.H b/integration/integrator_setup_strang.H
index 4afc1552bc..66df898327 100644
--- a/integration/integrator_setup_strang.H
+++ b/integration/integrator_setup_strang.H
@@ -90,8 +90,6 @@ IntegratorT integrator_setup (BurnT& state, amrex::Real dt, bool is_retry)
 
     burn_to_integrator(state, int_state);
 
-    // Save the initial composition, temperature, and energy for our later diagnostics.
-
     return int_state;
 }
 
diff --git a/networks/aprox13/actual_network.H b/networks/aprox13/actual_network.H
index 937585dd88..f678123274 100644
--- a/networks/aprox13/actual_network.H
+++ b/networks/aprox13/actual_network.H
@@ -134,7 +134,7 @@ namespace NSE_INDEX
 #endif
 
 namespace Rates {
-    enum NetworkRates {
+    enum NetworkRates: std::uint8_t {
         He4_He4_He4_to_C12 = 1,
         C12_He4_to_O16,
         C12_C12_to_Ne20_He4,
diff --git a/networks/aprox21/actual_network.H b/networks/aprox21/actual_network.H
index dcade97f08..9fca91944a 100644
--- a/networks/aprox21/actual_network.H
+++ b/networks/aprox21/actual_network.H
@@ -152,7 +152,7 @@ namespace NSE_INDEX
 #endif
 
 namespace Rates {
-    enum NetworkRates {
+    enum NetworkRates : std::uint8_t {
         H1_H1_to_He3 = 1,
         H1_H1_H1_to_He3,
         P_to_N,
diff --git a/networks/iso7/actual_network.H b/networks/iso7/actual_network.H
index ea60bd0351..82f81313cc 100644
--- a/networks/iso7/actual_network.H
+++ b/networks/iso7/actual_network.H
@@ -114,7 +114,7 @@ namespace NSE_INDEX
 #endif
 
 namespace Rates {
-    enum NetworkRates {
+    enum NetworkRates : std::uint8_t {
         C12_He4_to_O16 = 1,
         He4_He4_He4_to_C12,
         C12_C12_to_Ne20_He4,
diff --git a/networks/rprox/actual_network.H b/networks/rprox/actual_network.H
index c16d12dc61..4afd4fc8bb 100644
--- a/networks/rprox/actual_network.H
+++ b/networks/rprox/actual_network.H
@@ -112,7 +112,7 @@ namespace network
 
 namespace Rates
 {
-    enum NetworkRates
+    enum NetworkRates : std::uint8_t
     {
         He4_He4_He4_to_C12 = 1,
         C12_H1_to_N13,
diff --git a/networks/subch_base b/networks/subch_base
deleted file mode 120000
index 9bd398baed..0000000000
--- a/networks/subch_base
+++ /dev/null
@@ -1 +0,0 @@
-he-burn/he-burn-18a/
\ No newline at end of file
diff --git a/networks/subch_simple b/networks/subch_simple
deleted file mode 120000
index 43aa9e14b0..0000000000
--- a/networks/subch_simple
+++ /dev/null
@@ -1 +0,0 @@
-he-burn/he-burn-22a/
\ No newline at end of file
diff --git a/nse_tabular/nse_eos.H b/nse_tabular/nse_eos.H
index 425e0d70df..1006889325 100644
--- a/nse_tabular/nse_eos.H
+++ b/nse_tabular/nse_eos.H
@@ -65,7 +65,7 @@ nse_T_abar_from_e(const amrex::Real rho, const amrex::Real e_in, const amrex::Re
         amrex::Real abar_old = nse_state.abar;
 
         // call the EOS with the initial guess for T
-        eos_extra_t eos_state;
+        eos_re_extra_t eos_state;
         eos_state.rho = rho;
         eos_state.T = T;
         eos_state.aux[iye] = Ye;
@@ -156,7 +156,7 @@ nse_rho_abar_from_e(const amrex::Real T, const amrex::Real e_in, const amrex::Re
         amrex::Real abar_old = nse_state.abar;
 
         // call the EOS with the initial guess for rho
-        eos_extra_t eos_state;
+        eos_re_extra_t eos_state;
         eos_state.rho = rho;
         eos_state.T = T;
         eos_state.aux[iye] = Ye;
@@ -247,7 +247,7 @@ nse_T_abar_from_p(const amrex::Real rho, const amrex::Real p_in, const amrex::Re
         amrex::Real abar_old = nse_state.abar;
 
         // call the EOS with the initial guess for T
-        eos_extra_t eos_state;
+        eos_rep_extra_t eos_state;
         eos_state.rho = rho;
         eos_state.T = T;
         eos_state.aux[iye] = Ye;
@@ -338,7 +338,7 @@ nse_rho_abar_from_p(const amrex::Real T, const amrex::Real p_in, const amrex::Re
         amrex::Real abar_old = nse_state.abar;
 
         // call the EOS with the initial guess for rho
-        eos_extra_t eos_state;
+        eos_rep_extra_t eos_state;
         eos_state.rho = rho;
         eos_state.T = T;
         eos_state.aux[iye] = Ye;
diff --git a/sphinx_docs/source/data_structures.rst b/sphinx_docs/source/data_structures.rst
index a76cf1f0cc..8a4c4ccc3c 100644
--- a/sphinx_docs/source/data_structures.rst
+++ b/sphinx_docs/source/data_structures.rst
@@ -13,6 +13,8 @@ EOS
 ``eos_t``
 ---------
 
+.. index:: eos_t
+
 The main data structure for interacting with the EOS is ``eos_t``.
 This is a collection of data specifying the microphysical state of the
 fluid that we are evaluating. This has many components. For a
@@ -49,6 +51,8 @@ Networks
 ``burn_t``
 ----------
 
+.. index:: burn_t
+
 The main data structure for interacting with the reaction networks is
 ``burn_t``. This holds the composition (mass fractions), thermodynamic
 state, and a lot of internal information used by the reaction network
@@ -89,7 +93,9 @@ the user will only need to fill/use the following information:
 ``rate_t``, ``rate_fr_t``
 -------------------------
 
-The ``rate_t`` and ``rate_fr_t`` structures are used internally in a network to pass the
+.. index:: rate_t
+
+The ``rate_t`` structure is used internally in a network to pass the
 raw reaction rate information (usually just the temperature-dependent
 terms) between various subroutines. It does not come out of the
 network-specific righthand side or Jacobian routines.
@@ -97,6 +103,8 @@ network-specific righthand side or Jacobian routines.
 ``burn_type.H``
 ---------------
 
+.. index:: burn_type.H
+
 In addition to defining the ``burn_t`` type, the header ``burn_type.H``
 also defines integer indices into the solution vector that can be used
 to access the different components of the state:
diff --git a/sphinx_docs/source/eos.rst b/sphinx_docs/source/eos.rst
index f49d4987fa..dee9ee809f 100644
--- a/sphinx_docs/source/eos.rst
+++ b/sphinx_docs/source/eos.rst
@@ -9,6 +9,8 @@ an EOS module in case you want to build your own.
 Available Equations of State
 ============================
 
+.. index:: eos_t
+
 The following equations of state are available in Microphysics.
 Except where noted, each of these EOSs will provide the full
 thermodynamic data (including all derivatives) in the ``eos_t``
@@ -22,7 +24,7 @@ equation of state:
 
 .. math:: p = (\gamma - 1) \rho e.
 
-:math:`\gamma` is specified by the runtime parameter ``eos_gamma``. For
+:math:`\gamma` is specified by the runtime parameter ``eos.eos_gamma``. For
 an ideal gas, this represents the ratio of specific heats. The gas is
 assumed to be ideal, with the pressure given by
 
@@ -30,12 +32,12 @@ assumed to be ideal, with the pressure given by
 
 where :math:`k` is Boltzmann’s constant and :math:`\mu` is the mean molecular
 weight, calculated from the composition, :math:`X_k`. This EOS assumes
-the gas is either completely neutral (``assume_neutral = T``),
+the gas is either completely neutral (``eos.assume_neutral = 1``),
 giving:
 
 .. math:: \mu^{-1} = \sum_k \frac{X_k}{A_k}
 
-or completely ionized (``assume_neutral = F``), giving:
+or completely ionized (``eos.assume_neutral = 0``), giving:
 
 .. math:: \mu^{-1} = \sum_k \left ( 1 + Z_k \right ) \frac{X_k}{A_k}
 
@@ -43,11 +45,6 @@ The entropy comes from the Sackur-Tetrode equation. Because of the
 complex way that composition enters into the entropy, the entropy
 formulation here is only correct for a :math:`\gamma = 5/3` gas.
 
-Note that the implementation provided in Microphysics is the same as
-the version shipped with MAESTRO, but more general than the
-``gamma_law`` EOS provided with CASTRO. CASTRO’s default EOS only
-fills the thermodynamic information in ``eos_t`` that is required
-by the hydrodynamics module in CASTRO.
 
 polytrope
 ---------
@@ -63,21 +60,21 @@ only independent variable; there is no temperature dependence. The
 user either selects from a set of predefined options reflecting
 physical polytropes (e.g. a non-relativistic, fully degenerate
 electron gas) or inputs their own values for :math:`K` and :math:`\gamma`
-via ``polytrope_K`` and ``polytrope_gamma``.
+via ``eos.polytrope_K`` and ``eos.polytrope_gamma``.
 
-The runtime parameter ``polytrope_type`` selects the pre-defined
+The runtime parameter ``eos.polytrope_type`` selects the pre-defined
 polytropic relations. The options are:
 
--  ``polytrope_type = 1``: sets :math:`\gamma = 5/3` and
+-  ``eos.polytrope_type = 1``: sets :math:`\gamma = 5/3` and
 
    .. math:: K = \left ( \frac{3}{\pi} \right)^{2/3} \frac{h^2}{20 m_e m_p^{5/3}} \frac{1}{\mu_e^{5/3}}
 
-   where :math:`mu_e` is the mean molecular weight per electron, specified via ``polytrope_mu_e``
+   where :math:`mu_e` is the mean molecular weight per electron, specified via ``eos.polytrope_mu_e``
 
    This is the form appropriate for a non-relativistic
    fully-degenerate electron gas.
 
--  ``polytrope_type = 2``: sets :math:`\gamma = 4/3` and
+-  ``eos.polytrope_type = 2``: sets :math:`\gamma = 4/3` and
 
    .. math:: K = \left ( \frac{3}{\pi} \right)^{1/3} \frac{hc}{8 m_p^{4/3}} \frac{1}{\mu_e^{4/3}}
 
@@ -161,12 +158,12 @@ and :math:`p = \rho e (\gamma_\mathrm{effective} - 1)`.
 This equation of state takes several runtime parameters that can set
 the :math:`\gamma_i` for a specific species. The parameters are:
 
--  ``eos_gamma_default``: the default :math:`\gamma` to apply for all
+-  ``eos.eos_gamma_default``: the default :math:`\gamma` to apply for all
    species
 
--  ``species_X_name`` and ``species_X_gamma``: set the
+-  ``eos.species_X_name`` and ``eos.species_X_gamma``: set the
    :math:`\gamma_i` for the species whose name is given as
-   ``species_X_name`` to the value provided by ``species_X_gamma``.
+   ``eos.species_X_name`` to the value provided by ``eos.species_X_gamma``.
    Here, ``X`` can be one of the letters: ``a``, ``b``, or
    ``c``, allowing us to specify custom :math:`\gamma_i` for up to three
    different species.
@@ -209,6 +206,8 @@ appropriate interpolation table from that site to use this.
 Interface and Modes
 ===================
 
+.. index:: eos_t, eos_re_t, eos_rep_t, eos_rh_t, chem_eos_t
+
 The EOS is called as:
 
 .. code:: c++
@@ -244,6 +243,11 @@ The *eos_type* passed in is one of
 
 * ``eos_rep_t`` : expands on ``eos_re_t`` to include pressure information
 
+* ``eos_rh_t`` : expands on ``eos_rep_t`` to include enthalpy information
+
+* ``chem_eos_t`` : adds some quantities needed for the primordial chemistry EOS
+  and explicitly does not include the mass fractions.
+
 In general, you should use the type that has the smallest set of
 information needed, since we optimize out needless quantities at
 compile type (via C++ templating) for ``eos_re_t`` and ``eos_rep_t``.
@@ -260,6 +264,7 @@ compile type (via C++ templating) for ``eos_re_t`` and ``eos_rep_t``.
 Auxiliary Composition
 ---------------------
 
+.. index:: USE_AUX_THERMO
 
 With ``USE_AUX_THERMO=TRUE``, we interpret the composition from the auxiliary variables.
 The auxiliary variables are
@@ -298,6 +303,21 @@ The equation of state also needs :math:`\bar{Z}` which is easily computed as
    \bar{Z} = \bar{A} Y_e
 
 
+Composition Derivatives
+-----------------------
+
+.. index:: eos_extra_t, eos_re_extra_t, eos_rep_extra_t
+
+The derivatives $\partial p/\partial A$, $\partial p/\partial Z$,
+and $\partial e/\partial A$, $\partial e/\partial Z$ are available via
+the ``eos_extra_t``, ``eos_re_extra_t``, ``eos_rep_extra_t``, which
+extends the non-"extra" variants with these additional fields.
+
+The composition derivatives can be used via the ``composition_derivatives()`` function
+in ``eos_composition.H``
+to compute :math:`\partial p/\partial X_k |_{\rho, T, X_j}`, :math:`\partial e/\partial X_k |_{\rho, T, X_j}`, and :math:`\partial h/\partial X_k |_{\rho, T, X_j}`.
+
+
 Initialization and Cutoff Values
 ================================
 
@@ -328,6 +348,3 @@ appropriate time for, say, loading an interpolation table into memory.
 
 The main evaluation routine is called ``actual_eos``. It should
 accept an eos_input and an eos_t state; see Section :ref:`data_structures`.
-
-
-
diff --git a/sphinx_docs/source/getting_started.rst b/sphinx_docs/source/getting_started.rst
index a338f28ec6..b02c121249 100644
--- a/sphinx_docs/source/getting_started.rst
+++ b/sphinx_docs/source/getting_started.rst
@@ -5,6 +5,8 @@ Getting Started
 Standalone
 ==========
 
+.. index:: AMREX_HOME
+
 Microphysics can be used in a "standalone" fashion to run the unit
 tests and explore the behavior of the reaction networks.  The main
 requirement is a copy of AMReX:
@@ -23,6 +25,8 @@ to set the ``AMREX_HOME`` environment variable to point to the
 
 (where you change ``/path/to/amrex`` to your actual path).
 
+.. index:: burn_cell
+
 A good unit test to start with is ``burn_cell`` -- this is simply a
 one-zone burn.  In ``Microphysics/`` do:
 
@@ -45,6 +49,8 @@ Here ``inputs_aprox21`` is the inputs file that sets options.
 Running with AMReX Application Code
 ===================================
 
+.. index:: MICROPHYSICS_HOME
+
 Getting started with Microphysics using either CASTRO or MAESTROeX is
 straightforward. Because the modules here are already in a format that
 the AMReX codes understand, you only need to provide to the code
diff --git a/sphinx_docs/source/index.rst b/sphinx_docs/source/index.rst
index 7b48af3de1..c0e9fc7623 100644
--- a/sphinx_docs/source/index.rst
+++ b/sphinx_docs/source/index.rst
@@ -6,18 +6,36 @@
 AMReX-Astro Microphysics
 ************************
 
-A collection of microphysics routines (equations of state,
+AMReX-Astro Microphysics is a collection of microphysics routines (equations of state,
 reaction networks, ...) and utilities (ODE integrators, NSE solvers)
 for astrophysical simulation codes.
 
-.. toctree::
-   :maxdepth: 1
+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
+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>`_
+simulation code.
+
+While there are a number of unit tests that exercise the functionality,
+Microphysics is primarily intended to be used along with another simulation
+code.   At the moment, the interfaces and
+build stubs are compatible with the AMReX codes and use the AMReX build
+system.
+
+A number of the routines contained here we authored by other people.
+We bundle them here with permission, usually changing the interfaces
+to be compatible with our standardized interface. We in particular
+thank Frank Timmes for numerous reaction networks and his equation
+of state routines.
 
-   preface
 
 .. toctree::
    :maxdepth: 1
    :caption: Microphysics overview
+   :hidden:
 
    getting_started
    design
@@ -28,6 +46,7 @@ for astrophysical simulation codes.
 .. toctree::
    :maxdepth: 1
    :caption: EOS and transport
+   :hidden:
 
    eos
    transport
@@ -35,6 +54,7 @@ for astrophysical simulation codes.
 .. toctree::
    :maxdepth: 1
    :caption: Reaction networks
+   :hidden:
 
    networks-overview
    networks
@@ -44,6 +64,7 @@ for astrophysical simulation codes.
 .. toctree::
    :maxdepth: 1
    :caption: ODE integrators
+   :hidden:
 
    integrators
    ode_integrators
@@ -53,6 +74,7 @@ for astrophysical simulation codes.
 .. toctree::
    :maxdepth: 1
    :caption: Unit tests
+   :hidden:
 
    unit_tests
    comprehensive_tests
@@ -61,13 +83,13 @@ for astrophysical simulation codes.
 .. toctree::
    :maxdepth: 1
    :caption: References
+   :hidden:
 
    zreferences
 
+.. toctree::
+   :maxdepth: 1
+   :caption: Index
+   :hidden:
 
-Indices and tables
-==================
-
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`
+   genindex
diff --git a/sphinx_docs/source/integrators.rst b/sphinx_docs/source/integrators.rst
index cf6c87565d..80fa229d05 100644
--- a/sphinx_docs/source/integrators.rst
+++ b/sphinx_docs/source/integrators.rst
@@ -17,11 +17,13 @@ The equations we integrate to do a nuclear burn are:
    \frac{de}{dt} = f(\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 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.
+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
+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.
 
 In this system, :math:`e` is equal to the total specific internal
 energy. This allows us to easily call the EOS during the burn to obtain the temperature.
@@ -103,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``
@@ -153,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:
 
@@ -223,7 +246,7 @@ flow is (for VODE):
    and zero out the temperature and energy derivatives if we are not integrating
    those quantities.
 
-#. apply any boosting if ``react_boost`` > 0
+#. apply any boosting if ``integrator.react_boost`` > 0
 
 
 Jacobian implementation
@@ -243,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:
 
@@ -314,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:
@@ -329,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 f(\dot{Y}_k) dt
 
-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.
+.. 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.
+
+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``.
diff --git a/sphinx_docs/source/networks-overview.rst b/sphinx_docs/source/networks-overview.rst
index f4beb95e9a..02074c1a73 100644
--- a/sphinx_docs/source/networks-overview.rst
+++ b/sphinx_docs/source/networks-overview.rst
@@ -114,9 +114,11 @@ There are two primary files within each network directory.
    state and (respectively) the time-derivatives and Jacobian
    elements to fill in.
 
-   Note: some networks do not provide an analytic Jacobian and instead
-   rely on the numerical difference-approximation to the Jacobian. In
-   this case, the interface ``actual_jac`` is still needed to compile.
+   .. note::
+
+      Some networks do not provide an analytic Jacobian and instead
+      rely on the numerical difference-approximation to the Jacobian. In
+      this case, the interface ``actual_jac`` is still needed to compile.
 
 Notice that these modules have initialization routines:
 
diff --git a/sphinx_docs/source/networks.rst b/sphinx_docs/source/networks.rst
index e0fd944f00..660054d3a4 100644
--- a/sphinx_docs/source/networks.rst
+++ b/sphinx_docs/source/networks.rst
@@ -56,6 +56,8 @@ for plotfile variables, and the mass number, :math:`A`, and proton number, :math
 
 The name of the inputs file by one of two make variables:
 
+.. index:: NETWORK_INPUTS, GENERAL_NET_INPUTS
+
 * ``NETWORK_INPUTS`` : this is simply the name of the "`.net`" file, without
   any path.  The build system will look for it in the current directory
   and then in ``$(MICROPHYSICS_HOME)/networks/general_null/``.
diff --git a/sphinx_docs/source/nse.rst b/sphinx_docs/source/nse.rst
index 67b646b1f4..74d55f0ae8 100644
--- a/sphinx_docs/source/nse.rst
+++ b/sphinx_docs/source/nse.rst
@@ -2,13 +2,21 @@
 NSE
 ***
 
+.. important::
+
+   NSE is only supported with the simplified-SDC method for
+   coupling hydrodynamics and reactions.  We do not support
+   operator-splitting (Strang) coupling with NSE.
+
 The reaction networks in Microphysics have the ability to use NSE
 instead of integrating the entire network when the conditions are
 appropriate.  There are 2 different implementations of NSE in
 Microphysics, that have slightly different use cases.
 
+.. index:: USE_NSE_TABLE, USE_NSE_NET
+
 * :ref:`tabulated_nse` : this uses a table of NSE abundances given
-  :math:`(\rho, T, Y_e)` generate from a large network (125 isotopes).
+  :math:`(\rho, T, Y_e)` generate from a large network (96 isotopes).
   The table also returns :math:`dY_e/dt` resulting from
   electron-captures, to allow for the NSE state to evolve.  This is
   meant to be used in the cores of massive stars and works only with the
@@ -21,8 +29,11 @@ Microphysics, that have slightly different use cases.
   :math:`\langle B/A\rangle`.  All of the EOS calls will work with
   these quantities.
 
+  This algorithm was described in :cite:`sdc-nse`.
+
   This is enabled via ``USE_NSE_TABLE``
 
+
 * :ref:`self_consistent_nse` : this adds an NSE solver to the network that
   can be called to find the equilibrium abundances of each of the
   species defined in the network.  It works with any of the
@@ -59,18 +70,19 @@ standard ``aprox19`` network with a table for nuclear statistic
 equilibrium resulting from a much larger network at high density and
 temperatures.    This option is enabled by building with:
 
-.. prompt:: bash
+.. code:: bash
 
    NETWORK_DIR=aprox19 USE_NSE_TABLE=TRUE
 
 Composition and EOS
 -------------------
 
-The NSE table was generated using a 125 nuclei reaction network
-(described in :cite:`ma:2013`), and includes electron-capture rates,
-so the compositional quantities it carries, :math:`\bar{A}` and
-:math:`Y_e` and not representable from the 19 isotopes we carry in the
-network.  In particular, it can attain a lower :math:`Y_e` than
+The NSE table was generated using `pynucastro
+<https://pynucastro.github.io/pynucastro/>` using 96 nuclei and
+electron/positron capture/decay rates from :cite:`langanke:2001`.  The
+table takes $Y_e$ as the primary composition variable and provides a
+set of mass fractions that is mapped into those used by ``aprox19``.
+Using the value allows us to attain a lower :math:`Y_e` than
 ``aprox19`` can represent.
 
 For this reason, when we are using the NSE network, we always take the
@@ -84,35 +96,39 @@ NSE Table Outputs
 -----------------
 
 The NSE table provides values for the auxiliary composition,
-:math:`Y_e`, :math:`\bar{A}`, and :math:`\langle B/A \rangle`
-resulting from the full 125 nuclei network.   It also provides a set of 19
+:math:`\bar{A}`, and :math:`\langle B/A \rangle`
+resulting from the full 96 nuclei network.   It also provides a set of 19
 :math:`X_k` that map into the isotopes carried by ``aprox19``.
-
-
 These three quantities are stored as ``aux`` data in the network and
 are indexed as ``iye``, ``iabar``, and ``ibea``.  Additionally, when
 coupling to hydrodynamics, we need to advect these auxiliary
 quantities.
 
-For Strang split coupling of hydro and reactions, :math:`DX_k/Dt = 0`,
-and our evolution equations are:
+The evolution equations for the auxiliary variables are:
 
 .. math::
 
    \begin{align*}
-   \frac{DY_e}{Dt} &= \sum_k \frac{Z_k}{A_k} \frac{DX_k}{Dt} = 0 \\
-   \frac{D}{Dt} \frac{1}{\bar{A}} &= - \frac{1}{\bar{A}^2} \frac{D\bar{A}}{Dt} = \sum_k \frac{1}{A_k} \frac{DX_k}{Dt} = 0 \rightarrow \frac{D\bar{A}}{Dt} = 0 \\
-   \frac{D}{Dt} \left (\frac{B}{A} \right ) &= \sum_k \frac{B_k}{A_k} \frac{DX_k}{Dt} = 0
+   \frac{DY_e}{Dt} &= \sum_k \frac{Z_k}{A_k} \dot{\omega}_k \\
+   \frac{D\bar{A}}{Dt} &= -\bar{A}^2 \sum_k \frac{1}{A_k} \dot{\omega}_k \\
+   \frac{D}{Dt} \left (\frac{B}{A} \right ) &= \sum_k \frac{B_k}{A_k} \dot{\omega}_k
    \end{align*}
 
 Therefore each of these auxiliary equations obeys an advection equation
 in the hydro part of the advancement.
 
+The table also provides $dY_e/dt$, $(d\langle
+B/A\rangle/dt)_\mathrm{weak}$, and $\epsilon_{\nu,\mathrm{react}}$, the
+weak rate neutrino losses.  These quantities are used to update the
+thermodynamic state as we integrate.
 
 NSE Flow
 --------
 
-The basic flow of a simulation using ``aprox19`` + the NSE table is as follows:
+.. index:: integrator.nse_deriv_dt_factor, integrator.nse_include_enu_weak
+
+The time integration algorithm is described in detail in :cite:`sdc-nse`.  Here
+we provide an outline:
 
 * initialize the problem, including :math:`X_k`
 
@@ -127,34 +143,58 @@ The basic flow of a simulation using ``aprox19`` + the NSE table is as follows:
 
   * if we are in an NSE region:
 
-    * use :math:`\rho`, :math:`T`, and :math:`Y_e` to call the table.
-      This returns: :math:`dY_e/dt`, :math:`(B/A)_{\rm out}`, and :math:`\bar{A}_{\rm out}`.
+    * Compute the initial temperature given $\rho$, $e$, and $Y_e$,
+      using an EOS inversion algorithm that understands NSE (in
+      particular that the composition depends on $T$ in NSE)
 
-    * update :math:`Y_e` [#fY]_ :
+    * Compute the plasma neutrino losses, $\epsilon_{\nu,\mathrm{thermal}}$
 
-      .. math::
-
-         (Y_e)_{\rm out} = (Y_e)_{\rm in} + \Delta t \frac{dY_e}{dt}
+    * Use :math:`\rho`, :math:`T`, and :math:`Y_e` to evaluate the NSE
+      state and construct $[\Rb(\Uc^\prime)]^n$, the source term from reactions to the
+      reduced conserved state $\Uc^\prime$ (this is the state used by the SDC algorithm
+      and includes the internal energy density, mass fractions, and auxiliary variables).
 
-    * :math:`\bar{A}_{\rm out}` is simply the value returned from the table
+      This is done via finite differencing in time (through a step
+      $\tau \ll \Delta t$), and the reactive sources are constructed
+      to exclude the advective contributions.  The size of $\tau$ is
+      controlled via ``integrator.nse_deriv_dt_factor``.
 
-    * the energy generation rate, :math:`e_{\rm nuc}` is:
+      In particular, the energy source is constructed as:
 
       .. math::
 
-         e_{\rm nuc} = \eta \left [ \left ( \frac{B}{A} \right )_{\rm out} -
-                                    \left ( \frac{B}{A} \right )_{\rm in} \right ] * \frac{1.602 \times 10^{-6}  {\rm erg}}{{\rm MeV}} N_A \frac{1}{\Delta t}
+         R(\rho e) = N_A \frac{\Delta (\rho \langle B/A\rangle)}{\tau} + N_A \Delta m_{np} c^2 \rho \frac{dY_e}{dt} - \rho (\epsilon_{\nu,\mathrm{thermal}} + \epsilon_{\nu,\mathrm{react}})
+
+      where $\Delta m_{np}$ is the difference between the neutron and H atom mass.
 
+      .. important::
 
-      where :math:`\eta` is an inertia term < 1 to prevent the energy changing too much in one set.
+         It only makes sense to include the weak rate neutrino losses, $\epsilon_{\nu,\mathrm{react}}$,
+         if the initial model that you are using in your simulation also included those losses.
+         Otherwise, the energy loss from our NSE table will likely be too great and that simulation
+         will not be in equilibrium.  This is an issue, for example, when using a MESA model
+         constructed with ``aprox21``, which does not have all of the weak rates we model here.
 
-    * the new binding energy for the zone is then:
+         The weak rate neutrino losses can be disabled by ``integrator.nse_include_enu_weak=0``.
+
+    * Predict $\Uc^\prime$ to the midpoint in time, $n+1/2$ and construct
+      $[\Rb(\Uc^\prime)]^{n+1/2}$.
+
+    * Do the final update to time $n$ as:
 
       .. math::
 
-         \left ( \frac{B}{A} \right )_{\rm out}  = \left ( \frac{B}{A} \right )_{\rm in} + \eta \left [ \left ( \frac{B}{A} \right )_{\rm out} - \left ( \frac{B}{A} \right )_{\rm in} \right ]
+         \Uc^{\prime,n+1/2} = \Uc^{\prime,n} + \frac{\Delta t}{2} [\Advs{\Uc^\prime}]^{n+1/2} + \frac{\Delta t}{2} [\Rb(\Uc^\prime)]^{n+1/2}
+
+
+      where $[\Advs{\Uc^\prime}]^{n+1/2}$ are the advective updates carried by the SDC
+      algorithm.
+
+    * Compute the energy generation rate from the change in internal energy from $\Uc^{\prime,n}$ to $\Uc^{\prime,n+1}$, excluding advection.
+
+    * Update the total energy.
 
-    * update the mass fractions, :math:`X_k`, using the values from the table
+    * Set the mass fractions carried on the grid from the NSE table (with the new temperature and $Y_e$).
 
   * if we are not in NSE:
 
@@ -166,35 +206,38 @@ The basic flow of a simulation using ``aprox19`` + the NSE table is as follows:
 NSE check
 ---------
 
-For a zone to be consider in NSE, we require $\rho$ > ``rho_nse`` and *either*
+.. index:: network.rho_nse, network.T_nse, network.T_always_nse
+.. index:: network.He_Fe_nse, network.C_nse, network.O_nse, network.Si_nse
 
-* $T$ > ``T_nse`` together with the composition check
+For a zone to be consider in NSE, we require $\rho$ > ``network.rho_nse`` and *either*
 
-* $T$ > ``T_always_nse``
+* $T$ > ``network.T_nse`` together with the composition check
+
+* $T$ > ``network.T_always_nse``
 
 where we assume that ``T_always_nse`` > ``T_nse``.
 
 The composition check considers the following nuclei groups:
 
-* ``He_group``: atomic numbers 1 to 2 (H to He)
+* He-group: atomic numbers 1 to 2 (H to He)
 
-* ``C_group``: atomic numbers 6 to 7 (C to N)
+* C-group: atomic numbers 6 to 7 (C to N)
 
-* ``O_group``: atomic number 8 (O)
+* O-group: atomic number 8 (O)
 
-* ``Si_group``: atomic number 14 (Si)
+* Si-group: atomic number 14 (Si)
 
-* ``Fe_group``: atomic numbers 24 to 30 (Cr to Zn)
+* Fe-group: atomic numbers 24 to 30 (Cr to Zn)
 
 and we then say that a composition supports NSE if:
 
-* :math:`X(C_group)` < ``C_nse``
+* :math:`X(\mathrm{C}_\mathrm{group})` < ``network.C_nse``
 
-* :math:`X(O_group)` < ``O_nse``
+* :math:`X(\mathrm{O}_\mathrm{group})` < ``network.O_nse``
 
-* :math:`X(Si_group)` < ``Si_nse``
+* :math:`X(\mathrm{Si}_\mathrm{group})` < ``network.Si_nse``
 
-* :math:`X(Fe_group) + X(He_group)` > ``He_Fe_nse``
+* :math:`X(\mathrm{Fe}_\mathrm{group}) + X(\mathrm{He}_\mathrm{group})` > ``network.He_Fe_nse``
 
 
 
@@ -203,9 +246,9 @@ NSE table ranges
 
 The NSE table was created for:
 
-* :math:`9 < \log_{10}(T) < 10.4`
+* :math:`9.4 < \log_{10}(T) < 10.4`
 * :math:`7 < \log_{10}(\rho) < 10`
-* :math:`0.4 < Y_e < 0.5`
+* :math:`0.43 < Y_e < 0.5`
 
 
 
@@ -248,8 +291,8 @@ different syntax.
 The overall framework is constructed following :cite:`Kushnir_2020` with slight
 variations. The overview of the steps we take are the following:
 
-* Minimum Temperature Check: require ``T > T_min_nse``, where ``T_min_nse`` is
-  a runtime parameter with a default value ``T_min_nse = 4.0e9``.
+* Minimum Temperature Check: require ``T > nse.T_min_nse``, where ``nse.T_min_nse`` is
+  a runtime parameter with a default value ``nse.T_min_nse = 4.0e9``.
 
 * Mass Abundance Check: compare the current mass abundances of the nuclei to
   the NSE mass fractions. A detailed criteria are the following:
@@ -278,13 +321,13 @@ variations. The overview of the steps we take are the following:
 
   .. math::
 
-    \epsilon_{abs} = Y^i - Y^i_{NSE} < \mbox{nse_abs_tol}
+    \epsilon_{abs} = Y^i - Y^i_{NSE} < \mbox{nse.nse_abs_tol}
 
   .. math::
 
-    \epsilon_{rel} = \frac{\epsilon_{abs}}{Y^i} < \mbox{nse_rel_tol}
+    \epsilon_{rel} = \frac{\epsilon_{abs}}{Y^i} < \mbox{nse.nse_rel_tol}
 
-  where ``nse_rel_tol = 0.2`` and ``nse_abs_tol = 0.005`` by default.
+  where ``nse.nse_rel_tol = 0.2`` and ``nse.nse_abs_tol = 0.005`` by default.
 
 
 * **Removed** :cite:`Kushnir_2020` also requires a fast reaction cycle that
@@ -441,14 +484,3 @@ to the self-consistent nse check:
   the subsequent NSE checks. This is mainly to avoid unnecessary computations
   of computing the NSE mass fractions when the current temperature is too low.
   This is set to 4.0e9 by default.
-
-
-.. rubric:: Footnotes
-
-.. [#fY] The table actually provides the weak rate, which is the sum
-   of all electron capture and positron decay rates times the
-   appropriate abundances minus a similar rate for the beta decay and
-   positron capture, [wrate] = [rectot] + [rpdtot] - [redtot] - [rpctot]
-
-   So if electron capture dominates, then [wrate] is positive and this should
-   be subtracted from :math:`Y_e`.
diff --git a/sphinx_docs/source/ode_integrators.rst b/sphinx_docs/source/ode_integrators.rst
index 6a2286ec1e..fa7eb2b14c 100644
--- a/sphinx_docs/source/ode_integrators.rst
+++ b/sphinx_docs/source/ode_integrators.rst
@@ -1,9 +1,11 @@
 .. _ch:networks:integrators:
 
-*********************
-Available Integrators
-*********************
+***************
+ODE Integrators
+***************
 
+Available integrators
+=====================
 
 We use a high-order implicit ODE solver for integrating the reaction
 system.  A few alternatives, including first order implicit and explicit integrators are also
@@ -43,6 +45,8 @@ the allowed options are:
   the `Gershgorin circle theorem <https://en.wikipedia.org/wiki/Gershgorin_circle_theorem>`_
   is used instead.
 
+.. index:: integrator.use_jacobian_caching
+
 * ``VODE``: the VODE :cite:`vode` integration package.  We ported this
   integrator to C++ and removed the non-stiff integration code paths.
 
@@ -60,12 +64,6 @@ robust.
 
 .. index:: integrator.scale_system
 
-.. important::
-
-   The integrator will not abort if it encounters trouble.  Instead it will
-   set ``burn_t burn_state.success = false`` on exit.  It is up to the
-   application code to handle the failure.
-
 .. note::
 
    The runtime parameter ``integrator.scale_system``
@@ -81,6 +79,67 @@ robust.
 
    This option currently does not work with the ForwardEuler or QSS integrators.
 
+Timestep selection
+==================
+
+All of the integrators will select the timestep internally to meet the
+tolerances.  There are 2 controls that affect timestepping:
+
+* ``integrator.ode_max_dt`` : sets the maximum allowed timestep
+
+* ``integrator.ode_max_steps`` : sets the maximum number of steps
+  the integrator is allowed to take.  If it exceeds this, then
+  it will return an error.
+
+
+Linear algebra
+==============
+
+All implicit integrators use the LINPACK LU decomposition routines.
+
+For the templated networks (``aprox13``, ``aprox19``, ...) the implementation
+is done using ``consexpr`` loops over the equations and no pivoting is allowed.
+
+.. index:: integrator.linalg_do_pivoting
+
+For the other networks (usually pynucastro networks), the implementation is
+provided in ``Microphysics/util/linpack.H`` and is templated on the number
+of equations.  Pivoting can be disabled by setting ``integrator.linalg_do_pivoting=0``.
+
+Integration errors
+==================
+
+.. important::
+
+   The integrator will not abort if it encounters trouble.  Instead it will
+   set ``burn_t burn_state.success = false`` on exit.  It is up to the
+   application code to handle the failure.
+
+The ``burn_t`` ``error_code`` field will provide an error code that can be
+used to interpret the failure.  The current codes are:
+
++-------+----------------------------------------------------------+
+| code  | meaning                                                  |
++=======+==========================================================+
+| 1     | success                                                  |
++-------+----------------------------------------------------------+
+| -1    | invalid inputs                                           |
++-------+----------------------------------------------------------+
+| -2    | underflow in computing  $\Delta t$                       |
++-------+----------------------------------------------------------+
+| -3    | spectral radius estimation did not converge              |
++-------+----------------------------------------------------------+
+| -4    | too many steps needed                                    |
++-------+----------------------------------------------------------+
+| -5    | unable to meet the accuracy demanded by the tolerances   |
++-------+----------------------------------------------------------+
+| -6    | non-convergence in the corrector iteration               |
++-------+----------------------------------------------------------+
+| -7    | LU decomposition failed                                  |
++-------+----------------------------------------------------------+
+| -100  | entered NSE                                              |
++-------+----------------------------------------------------------+
+
 Tolerances
 ==========
 
@@ -131,6 +190,8 @@ is used for the temperature and energy.
 Controlling Species $\sum_k X_k = 1$
 ====================================
 
+.. index:: integrator.renormalize_abundances, integrator.SMALL_X_SAFE, integrator.do_species_clip
+
 The ODE integrators don't know about the constraint that
 
 $$\sum_k X_k = 1$$
@@ -164,6 +225,8 @@ constraint on the intermediate states during the integration.
 Retry Mechanism
 ===============
 
+.. index:: integrator.ode_max_steps
+
 Integration can fail for a number of reasons.  Some of the errors you may see are:
 
 1. Not enough steps allowed (``integrator.ode_max_steps``)
@@ -232,17 +295,3 @@ The runtime parameters that come into play when doing the retry are:
    ``integrator.ode_max_steps`` to a small value (like ``10000``) and
    start with the analytic Jacobian (``integrator.jacobian = 1``) and
    then use the retry mechanism to swap the Jacobian on any zones that fail.
-
-
-Overriding Parameter Defaults on a Network-by-Network Basis
-===========================================================
-
-Any network can override or add to any of the existing runtime
-parameters by creating a ``_parameters`` file in the network directory
-(e.g., ``networks/triple_alpha_plus_cago/_parameters``). As noted in
-:doc:`rp_intro`, the fourth column in the ``_parameter``
-file definition is the *priority*. When a duplicate parameter is
-encountered by the scripts writing the runtime parameter header files, the value
-of the parameter with the highest priority is used. So picking a large
-integer value for the priority in a network’s ``_parameter`` file will
-ensure that it takes precedence.
diff --git a/sphinx_docs/source/one_zone_tests.rst b/sphinx_docs/source/one_zone_tests.rst
index ec9706a05e..9352a2e626 100644
--- a/sphinx_docs/source/one_zone_tests.rst
+++ b/sphinx_docs/source/one_zone_tests.rst
@@ -9,6 +9,8 @@ on a single zone to inspect the output directly.
 ``burn_cell``
 =============
 
+.. index:: ``burn_cell``
+
 ``burn_cell`` is a simple one-zone burn that will evolve a state with
 a network for a specified amount of time.  This can be used to
 understand the timescales involved in a reaction sequence or to
diff --git a/sphinx_docs/source/preface.rst b/sphinx_docs/source/preface.rst
deleted file mode 100644
index efb680fb48..0000000000
--- a/sphinx_docs/source/preface.rst
+++ /dev/null
@@ -1,30 +0,0 @@
-*******
-Preface
-*******
-
-Welcome to the AMReX-Astro Microphysics!
-
-In this User’s Guide we describe the microphysics modules designed to
-enable simulations of stellar explosions.
-
-The original design was to support the AMReX codes CASTRO and
-MAESTRO. These all have a consistent interface and are designed to
-provide the users of those codes an easy experience in moving from the
-barebones microphysics modules provided in those codes. For the
-purposes of this user’s guide, the microphysical components we
-currently deal with are the equation of state (EOS) and the nuclear
-burning network.
-
-Microphysics is not a stand-alone code. It is intended to be used in
-conjunction with a simulation code. At the moment, the interfaces and
-build stubs are compatible with the AMReX codes. In many cases we
-will provide test modules that demonstrate a minimal working example
-for how to run the modules (leveraging the AMReX build system). The
-goal is to make the support more general, and extend to other codes
-in the future.
-
-A number of the routines contained here we authored by other people.
-We bundle them here with permission, usually changing the interfaces
-to be compatible with our standardized interface. We in particular
-thank Frank Timmes for numerous reaction networks and his equation
-of state routines.
diff --git a/sphinx_docs/source/refs.bib b/sphinx_docs/source/refs.bib
index 4b92f50b21..d0f234fbc5 100644
--- a/sphinx_docs/source/refs.bib
+++ b/sphinx_docs/source/refs.bib
@@ -625,3 +625,50 @@ @misc{autodiff
     howpublished = {\texttt{https://autodiff.github.io}},
     year = {2018}
 }
+
+
+@article{sdc-nse,
+doi = {10.3847/1538-4357/ad8a66},
+url = {https://dx.doi.org/10.3847/1538-4357/ad8a66},
+year = {2024},
+month = {nov},
+publisher = {The American Astronomical Society},
+volume = {977},
+number = {1},
+pages = {30},
+author = {Michael Zingale and Zhi Chen and Eric T. Johnson and Max P. Katz and Alexander Smith Clark},
+title = {Strong Coupling of Hydrodynamics and Reactions in Nuclear Statistical Equilibrium for Modeling Convection in Massive Stars},
+journal = {The Astrophysical Journal},
+abstract = {We build on the simplified spectral deferred corrections
+                  (SDC) coupling of hydrodynamics and reactions to
+                  handle the case of nuclear statistical equilibrium
+                  (NSE) and electron/positron captures/decays in the
+                  cores of massive stars. Our approach blends a
+                  traditional reaction network on the grid with a
+                  tabulated NSE state from a very large, nuclei
+                  network. We demonstrate how to achieve second-order
+                  accuracy in the simplified-SDC framework when
+                  coupling NSE to hydrodynamics, with the ability to
+                  evolve the star on the hydrodynamics time step. We
+                  discuss the application of this method to convection
+                  in massive stars leading up to core collapse. We
+                  also show how to initialize the initial convective
+                  state from a 1D model in a self-consistent
+                  fashion. All of these developments are done in the
+                  publicly available Castro simulation code and the
+                  entire simulation methodology is fully
+                  GPU-accelerated.}
+}
+
+@article{langanke:2001,
+title = {RATE TABLES FOR THE WEAK PROCESSES OF pf-SHELL NUCLEI IN STELLAR ENVIRONMENTS},
+journal = {Atomic Data and Nuclear Data Tables},
+volume = {79},
+number = {1},
+pages = {1-46},
+year = {2001},
+issn = {0092-640X},
+doi = {https://doi.org/10.1006/adnd.2001.0865},
+author = {K. LANGANKE and G. MARTÍNEZ-PINEDO},
+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.}
+}
\ No newline at end of file
diff --git a/sphinx_docs/source/screening.rst b/sphinx_docs/source/screening.rst
index 1d99a925c7..d1739aeb86 100644
--- a/sphinx_docs/source/screening.rst
+++ b/sphinx_docs/source/screening.rst
@@ -2,6 +2,8 @@
 Screening of Reaction Rates
 ***************************
 
+.. index:: SCREEN_METHOD
+
 Screening of reaction rates can be computed using several different methods,
 controlled by the make parameter ``SCREEN_METHOD``.  For example,
 
@@ -51,4 +53,12 @@ The options are:
 
 Runtime Options
 ----------------
-* ``screening.enable_chabrier1998_quantum_corr = 1`` in the input file enables an additional quantum correction term added to the screening factor when ``SCREEN_METHOD=chabrier1998``. This is disabled by default since ``chabrier1998`` is often used along with ``USE_NSE_NET=TRUE``, and the NSE solver doesn't include quantum corrections.
+
+.. index:: screening.enable_chabrier1998_quantum_corr
+
+* ``screening.enable_chabrier1998_quantum_corr = 1`` in the input file
+  enables an additional quantum correction term added to the screening
+  factor when ``SCREEN_METHOD=chabrier1998``. This is disabled by
+  default since ``chabrier1998`` is often used along with
+  ``USE_NSE_NET=TRUE``, and the NSE solver doesn't include quantum
+  corrections.
diff --git a/sphinx_docs/source/sdc.rst b/sphinx_docs/source/sdc.rst
index 3cbf8e7c32..38821a2a10 100644
--- a/sphinx_docs/source/sdc.rst
+++ b/sphinx_docs/source/sdc.rst
@@ -33,6 +33,7 @@ hydrodynamical sources),
 and :math:`\Rb(\Uc)`
 is the reaction source term.
 
+.. index:: USE_TRUE_SDC, USE_SIMPLIFIED_SDC
 
 .. note::
 
diff --git a/sphinx_docs/source/unit_tests.rst b/sphinx_docs/source/unit_tests.rst
index a95a18de72..5a36b2c4af 100644
--- a/sphinx_docs/source/unit_tests.rst
+++ b/sphinx_docs/source/unit_tests.rst
@@ -43,6 +43,8 @@ Tests are divided into three categories:
 Comprehensive tests
 ===================
 
+.. index:: test_aprox_rates, test_conductivity, test_eos, test_jac, test_neutrino_cooling, test_react, test_rhs, test_screening_templated, test_sdc
+
 Each of these tests sets up a cube of data, $(\rho, T, X_k)$, with the
 range of $T$ and $\rho$, and the species to focus on for $X_k$ controlled
 by options in the input file.
@@ -113,6 +115,8 @@ by options in the input file.
 One-zone tests
 ==============
 
+.. index:: burn_cell, burn_cell_primordial_chem, burn_cell_sdc, eos_cell, jac_cell, nse_table_cell, test_ase, test_part_func
+
 * ``burn_cell`` :
 
   given a $\rho$, $T$, and $X_k$, integrate a reaction network through a specified time
@@ -157,6 +161,8 @@ One-zone tests
 Infrastructure tests
 ====================
 
+.. index:: test_linear_algebra, test_nse_interp, test_parameters, test_sdc_vode_rhs
+
 * ``test_linear_algebra`` :
 
   create a diagonally dominant matrix, multiply it by a test vector, $x$,
diff --git a/unit_test/burn_cell/burn_cell.H b/unit_test/burn_cell/burn_cell.H
index ad7771e7bc..dd35a9fa95 100644
--- a/unit_test/burn_cell/burn_cell.H
+++ b/unit_test/burn_cell/burn_cell.H
@@ -23,6 +23,7 @@ void burn_cell_c()
 
         if (unit_test_rp::init_species_all_equal) {
             massfractions[n-1] = 1.0_rt / static_cast<Real>(NumSpec);
+            massfractions[n-1] = 1.0_rt / static_cast<amrex::Real>(NumSpec);
         } else {
             massfractions[n-1] = get_xn(n);
 
diff --git a/unit_test/test_part_func/ci-benchmarks/part_func.out b/unit_test/test_part_func/ci-benchmarks/part_func.out
new file mode 100644
index 0000000000..6bb827da07
--- /dev/null
+++ b/unit_test/test_part_func/ci-benchmarks/part_func.out
@@ -0,0 +1,9 @@
+Initializing AMReX (24.10-20-gb9d549bcf4a6)...
+AMReX (24.10-20-gb9d549bcf4a6) initialized
+starting the single zone burn...
+reading in network electron-capture / beta-decay tables...
+temperature = 5000000000
+Ni56: 1.010868 2.337172328e-11
+Fe52: 1.743116 3.226066771e-10
+spins: 2 1 1 8 1
+AMReX (24.10-20-gb9d549bcf4a6) finalized