From 2afa0434d43fc2fbc68c01f1f1c63882636f9d19 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 5 Dec 2024 08:10:10 -0500 Subject: [PATCH 1/6] update documentation on integrators --- sphinx_docs/source/ode_integrators.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/sphinx_docs/source/ode_integrators.rst b/sphinx_docs/source/ode_integrators.rst index 6a2286ec1e..5304e3e673 100644 --- a/sphinx_docs/source/ode_integrators.rst +++ b/sphinx_docs/source/ode_integrators.rst @@ -43,6 +43,8 @@ the allowed options are: the `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. @@ -131,6 +133,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 +168,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``) From 1a81667deb6d95ae63f58ee85c25d2baf1b48e7a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 5 Dec 2024 10:51:09 -0500 Subject: [PATCH 2/6] lots more integrator doc fixest --- integration/integrator_setup_strang.H | 2 - sphinx_docs/source/integrators.rst | 99 +++++++++++++------- sphinx_docs/source/nse.rst | 125 +++++++++++++++++--------- sphinx_docs/source/refs.bib | 48 ++++++++++ 4 files changed, 196 insertions(+), 78 deletions(-) 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/sphinx_docs/source/integrators.rst b/sphinx_docs/source/integrators.rst index 77c6b93063..7c1164066d 100644 --- a/sphinx_docs/source/integrators.rst +++ b/sphinx_docs/source/integrators.rst @@ -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. + +#. 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& ydot) + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void actual_rhs(burn_t& state, amrex::Array1D& 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 + 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 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/nse.rst b/sphinx_docs/source/nse.rst index f0dfb51cf7..031847facf 100644 --- a/sphinx_docs/source/nse.rst +++ b/sphinx_docs/source/nse.rst @@ -2,6 +2,12 @@ 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 @@ -23,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 @@ -61,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 +` 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 @@ -86,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` +:math:`\bar{A}`, and :math:`\langle B/A \rangle` resulting from the full 125 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` @@ -129,34 +143,56 @@ 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` : + * Compute the plasma neutrino losses, $\epsilon_{\nu,\mathrm{thermal}}$ - .. math:: + * 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). - (Y_e)_{\rm out} = (Y_e)_{\rm in} + \Delta t \frac{dY_e}{dt} + 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``. - * :math:`\bar{A}_{\rm out}` is simply the value returned from the table - - * 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:: + + 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 weak rate neutrino losses can be disabled by ``integrator.nse_include_enu_weak``. - where :math:`\eta` is an inertia term < 1 to prevent the energy changing too much in one set. + * Predict $\Uc^\prime$ to the midpoint in time, $n+1/2$ and construct + $[\Rb(\Uc^\prime)]^{n+1/2}$. - * the new binding energy for the zone is then: + * 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} - * update the mass fractions, :math:`X_k`, using the values from the table + + 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. * if we are not in NSE: @@ -168,35 +204,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 + +For a zone to be consider in NSE, we require $\rho$ > ``network.rho_nse`` and *either* -* $T$ > ``T_nse`` together with the composition check +* $T$ > ``network.T_nse`` together with the composition check -* $T$ > ``T_always_nse`` +* $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_\mathrm{group})` < ``C_nse`` +* :math:`X(C_\mathrm{group})` < ``network.C_nse`` -* :math:`X(O_\mathrm{group})` < ``O_nse`` +* :math:`X(O_\mathrm{group})` < ``network.O_nse`` -* :math:`X(Si_\mathrm{group})` < ``Si_nse`` +* :math:`X(Si_\mathrm{group})` < ``network.Si_nse`` -* :math:`X(Fe_\mathrm{group}) + X(He_\mathrm{group})` > ``He_Fe_nse`` +* :math:`X(Fe_\mathrm{group}) + X(He_\mathrm{group})` > ``network.He_Fe_nse`` @@ -205,9 +244,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` diff --git a/sphinx_docs/source/refs.bib b/sphinx_docs/source/refs.bib index 4b92f50b21..3e871c8c05 100644 --- a/sphinx_docs/source/refs.bib +++ b/sphinx_docs/source/refs.bib @@ -625,3 +625,51 @@ @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}, +url = {https://www.sciencedirect.com/science/article/pii/S0092640X01908654}, +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 From d4ee01fe31f34e8667a5e594987e623fa5ac2072 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 5 Dec 2024 11:14:14 -0500 Subject: [PATCH 3/6] finish --- sphinx_docs/source/ode_integrators.rst | 89 +++++++++++++++++++------- 1 file changed, 66 insertions(+), 23 deletions(-) diff --git a/sphinx_docs/source/ode_integrators.rst b/sphinx_docs/source/ode_integrators.rst index 5304e3e673..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 @@ -62,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`` @@ -83,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 ========== @@ -238,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. From 4dc4356fe4dd543b2ecc3abf293f73641297ec89 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 5 Dec 2024 11:26:48 -0500 Subject: [PATCH 4/6] some fixes --- sphinx_docs/source/nse.rst | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/sphinx_docs/source/nse.rst b/sphinx_docs/source/nse.rst index 031847facf..bf77f55691 100644 --- a/sphinx_docs/source/nse.rst +++ b/sphinx_docs/source/nse.rst @@ -118,7 +118,7 @@ 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 +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. @@ -229,13 +229,13 @@ The composition check considers the following nuclei groups: and we then say that a composition supports NSE if: -* :math:`X(C_\mathrm{group})` < ``network.C_nse`` +* :math:`X(\mathrm{C}_\mathrm{group})` < ``network.C_nse`` -* :math:`X(O_\mathrm{group})` < ``network.O_nse`` +* :math:`X(\mathrm{O}_\mathrm{group})` < ``network.O_nse`` -* :math:`X(Si_\mathrm{group})` < ``network.Si_nse`` +* :math:`X(\mathrm{Si}_\mathrm{group})` < ``network.Si_nse`` -* :math:`X(Fe_\mathrm{group}) + X(He_\mathrm{group})` > ``network.He_Fe_nse`` +* :math:`X(\mathrm{Fe}_\mathrm{group}) + X(\mathrm{He}_\mathrm{group})` > ``network.He_Fe_nse`` @@ -289,8 +289,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: @@ -319,13 +319,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 From f453433e58332b459b8c40675864e430b6028ee3 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 5 Dec 2024 11:34:01 -0500 Subject: [PATCH 5/6] fix refs --- sphinx_docs/source/refs.bib | 1 - 1 file changed, 1 deletion(-) diff --git a/sphinx_docs/source/refs.bib b/sphinx_docs/source/refs.bib index 3e871c8c05..d0f234fbc5 100644 --- a/sphinx_docs/source/refs.bib +++ b/sphinx_docs/source/refs.bib @@ -669,7 +669,6 @@ @article{langanke:2001 year = {2001}, issn = {0092-640X}, doi = {https://doi.org/10.1006/adnd.2001.0865}, -url = {https://www.sciencedirect.com/science/article/pii/S0092640X01908654}, 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 From 0d7e54afa6e5afddccb66db0f402a09b12edf713 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 7 Dec 2024 07:20:15 -0500 Subject: [PATCH 6/6] address reviewer comments --- sphinx_docs/source/integrators.rst | 2 +- sphinx_docs/source/nse.rst | 12 +++++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/sphinx_docs/source/integrators.rst b/sphinx_docs/source/integrators.rst index 7c1164066d..80fa229d05 100644 --- a/sphinx_docs/source/integrators.rst +++ b/sphinx_docs/source/integrators.rst @@ -138,7 +138,7 @@ the functions ``integrator_setup()`` and ``integrator_cleanup()``. off the energy offset, the energy stored is now just that generated by reactions. -#. Normalize the abundances so they sum to 1. +#. 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 diff --git a/sphinx_docs/source/nse.rst b/sphinx_docs/source/nse.rst index bf77f55691..74d55f0ae8 100644 --- a/sphinx_docs/source/nse.rst +++ b/sphinx_docs/source/nse.rst @@ -16,7 +16,7 @@ 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 @@ -97,7 +97,7 @@ NSE Table Outputs The NSE table provides values for the auxiliary composition, :math:`\bar{A}`, and :math:`\langle B/A \rangle` -resulting from the full 125 nuclei network. It also provides a set of 19 +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 @@ -163,19 +163,19 @@ we provide an outline: .. math:: - 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}}) + 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:: - It only makes sense to include the weak rate neutrino losses, $\epsilon{\nu,\mathrm{react}}$, + 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 weak rate neutrino losses can be disabled by ``integrator.nse_include_enu_weak``. + 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}$. @@ -194,6 +194,8 @@ we provide an outline: * Update the total energy. + * 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: * integrate the ``aprox19`` network as usual