update the ODE integrators docs on tolerances (#1703)

this regenerates the plot comparing tolerances and includes it
now as a script under burn_cell

Docs/source/ode_integrators.rst

@@ -146,46 +146,80 @@ Tolerances
 Tolerances dictate how accurate the ODE solver must be while solving
 equations during a simulation.  Typically, the smaller the tolerance
 is, the more accurate the results will be.  However, if the tolerance
-is too small, the code may run for too long or the ODE solver will
-never converge.  In these simulations, ``rtol`` values will set the
-relative tolerances and ``atol`` values will set the absolute tolerances
-for the ODE solver.  Often, one can find and set these values in an
-input file for a simulation.
-
-:numref:`fig:tolerances` shows the results of a simple simulation using the
-burn_cell unit test to determine
-what tolerances are ideal for simulations.
-For this investigation, it was assumed that a run with a tolerance of :math:`10^{-12}`
-corresponded to an exact result,
-so it is used as the basis for the rest of the tests.
-From the figure, one can infer that the :math:`10^{-3}` and :math:`10^{-6}` tolerances
-do not yield the most accurate results
-because their relative error values are fairly large.
-However, the test with a tolerance of :math:`10^{-9}` is accurate
-and not so low that it takes incredible amounts of computer time,
-so :math:`10^{-9}` should be used as the default tolerance in future simulations.
+is too small, the code may run for too long, the ODE solver will
+never converge, or it might require at timestep that underflows.
+
+.. index:: integrator.rtol_spec, integrator.rtol_enuc, integrator.atol_spec, integrator.atol_enuc
+
+There are separate tolerances for the mass fractions and the energy,
+and there are both relative and absolute tolerances which act together.
+The tolerances are:
+
+* ``integrator.rtol_spec`` : the relative tolerance for the species
+  (mass fractions when running with Strang and partial densities when
+  running with SDC).
+
+* ``integrator.rtol_enuc`` : the relative tolerance on the energy
+  (specific internal energy when running with Strang, internal energy
+  density when running with SDC).
+
+* ``integrator.atol_spec`` : the absolute tolerance for the species
+  (this is always interpreted in terms of mass fraction and the appropriate
+  density weighting will be added for SDC).
+
+* ``integrator.atol_enuc`` : the absolute tolerance for energy -- this
+  is generally not interesting, since the energy is so large and therefore
+  best served via a relative tolerance.
+
+The tolerances are combined, e.g. for species, as:
+
+.. math::
+
+   \epsilon_{\mathrm{total}, k} = \epsilon_\mathrm{abs} + \epsilon_\mathrm{rel} |X_k|
+
+so if the mass fraction, $X_k$, is very small, then the absolute tolerance
+will set the error that the integrator tries to achieve.  If the mass fraction
+is large, $\mathcal{O}(1)$, then the relative tolerance dominates.
+
+Some suggestions when setting tolerances:
+
+.. index:: integrator.X_reject_buffer
+
+* If a burn does not converge with one type of Jacobian (analytic or
+  numerical) then it may do better with the other type.  This can be
+  automated via the ``integrator.use_burn_retry`` mechanism described
+  above.
+
+* Sometimes a burn completes better if the absolute tolerances are
+  made even smaller -- this will require the integrator to track trace
+  species better which can help with equilibrium better.
+
+* The VODE integrator has additional logic meant to ensure that
+  species don't change too much per timestep.  This is controlled by
+  ``integrator.X_reject_buffer``.  If a species $k$, has a mass
+  fraction $X_k > \mbox{X_reject_buffer} \cdot \mbox{atol_spec}$ then
+  we reject a VODE timestep if the mass fraction changes by more than
+  a factor of 4 in a single VODE timestep and we try again.  This is
+  all done internally to VODE.  Making ``X_reject_buffer`` larger will
+  allow it to ignore more trace species.
+
+Below is a comparison of how the tolerances affect the nucleosynthesis.
+This is run using ``burn_cell`` and the ``aprox13`` network.  Four separate
+runs were done, using tolerances of $10^{-3}$, $10^{-5}$, $10^{-8}$, and $10^{-12}$
+(all 4 tolerance parameters were set to the same value).  The run with the tightest
+tolerances ($10^{-12}$) was taken as the reference and relative errors were
+computed with respect to it.  The scripts for this are in ``Microphysics/unit_test/burn_cell/compare_tolerances/``.
 
 .. _fig:tolerances:
-.. figure:: tolerances.png
-   :alt: Relative error plot
+.. figure:: tolerance-compare.png
+   :alt: Relative error in mass fractions
    :width: 100%
 
    Relative error of runs with varying tolerances as compared
    to a run with an ODE tolerance of :math:`10^{-12}`.
 
-The integration tolerances for the burn are controlled by
-``rtol_spec`` and  ``rtol_enuc``,
-which are the relative error tolerances for
-:eq:`eq:spec_integrate` and :eq:`eq:enuc_integrate`,
-respectively. There are corresponding
-``atol`` parameters for the absolute error tolerances. Note that
-not all integrators handle error tolerances the same way—see the
-sections below for integrator-specific information.
-
-The absolute error tolerances are set by default
-to :math:`10^{-12}` for the species, and a relative tolerance of :math:`10^{-6}`
-is used for the temperature and energy.
-
+We see that using a tolerance of $10^{-5}$ generally gives reasonable mass
+fractions.  Looser than this can produce large errors.
 
 Controlling Species $\sum_k X_k = 1$
 ====================================
@@ -222,6 +256,7 @@ constraint on the intermediate states during the integration.
   This is enabled by default.
 
 
+
 Retry Mechanism
 ===============
 

Docs/source/runtime_parameters.rst

@@ -767,30 +767,6 @@ namespace: ``network``
 
 
 
-**NETWORK_DIR=he-burn.bak/he-burn-18a:**
-
-+---------------------------------------+---------------------------------------------------------+------------------------------+
-| parameter                             | description                                             | default value                |
-+=======================================+=========================================================+==============================+
-| ``disable_p_C12_to_N13``              |                                                         | 0                            |
-+---------------------------------------+---------------------------------------------------------+------------------------------+
-| ``disable_He4_N13_to_p_O16``          |                                                         | 0                            |
-+---------------------------------------+---------------------------------------------------------+------------------------------+
-
-
-
-**NETWORK_DIR=he-burn.bak/he-burn-22a:**
-
-+---------------------------------------+---------------------------------------------------------+------------------------------+
-| parameter                             | description                                             | default value                |
-+=======================================+=========================================================+==============================+
-| ``disable_p_C12_to_N13``              |                                                         | 0                            |
-+---------------------------------------+---------------------------------------------------------+------------------------------+
-| ``disable_He4_N13_to_p_O16``          |                                                         | 0                            |
-+---------------------------------------+---------------------------------------------------------+------------------------------+
-
-
-
 **NETWORK_DIR=metal_chem:**
 
 +---------------------------------------+---------------------------------------------------------+------------------------------+

unit_test/burn_cell/compare_tolerances/README

@@ -0,0 +1,10 @@
+# Comparing tolerances
+
+This is used the make the plot in the docs that compares tolerances.
+
+Compile `burn_cell` with `aprox13` and run the tests using the script
+`run_tol_compare.sh`.  This will run 4 different instances with different
+tolerances.
+
+A plot comparing the relative error of the simulations (using the run with
+the tightest tolerances as the reference) can be made with `plot-tols.py`.

unit_test/burn_cell/compare_tolerances/plot-tols.py

@@ -0,0 +1,49 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+# read in the output from different tolerances
+tola = np.loadtxt("state_over_time_tol_a.txt")
+tolb = np.loadtxt("state_over_time_tol_b.txt")
+tolc = np.loadtxt("state_over_time_tol_c.txt")
+ref = np.loadtxt("state_over_time_tol_d.txt")
+
+labela = r"tol = $10^{-3}$"
+labelb = r"tol = $10^{-5}$"
+labelc = r"tol = $10^{-8}$"
+
+# compute the relative error with respect to the
+# tighest tolerance
+dta = np.abs(tola - ref) / np.abs(ref)
+dtb = np.abs(tolb - ref) / np.abs(ref)
+dtc = np.abs(tolc - ref) / np.abs(ref)
+
+# plot the output
+fig, ax = plt.subplots()
+
+itemp = 1
+ihe4 = 2
+ic12 = 3
+isi28 = 7
+ini56 = 14
+
+fig, axes = plt.subplots(2, 2)
+
+for n, (ax, idx, name) in enumerate(zip(axes.flat, [ihe4, ic12, isi28, ini56], ["He4", "C12", "Si28", "Ni56"])):
+
+    ax.plot(tola[:, 0], dta[:, idx], label=labela)
+    ax.plot(tolb[:, 0], dtb[:, idx], label=labelb)
+    ax.plot(tolc[:, 0], dtc[:, idx], label=labelc)
+
+    ax.set_xlabel("time (s)")
+    ax.set_ylabel(f"relative error {name}")
+
+    ax.set_xscale("log")
+    ax.set_yscale("log")
+
+    if n == 0:
+        ax.legend(fontsize="small")
+
+fig.tight_layout()
+fig.set_size_inches(8, 7)
+
+fig.savefig("tolerance-compare.png")

unit_test/burn_cell/compare_tolerances/run_tol_compare.sh

@@ -0,0 +1,20 @@
+#!/bin/sh
+
+RUNPARAMS="
+unit_test.tmax=1.e-4
+unit_test.nsteps=1000
+unit_test.density=1.e8
+unit_test.temperature=3.e9"
+
+./main3d.gnu.ex inputs_aprox13 ${RUNPARAMS} integrator.rtol_spec=1.0e-3 integrator.rtol_enuc=1.0e-3 integrator.atol_spec=1.0e-3 integrator.atol_enuc=1.0e-3
+mv state_over_time.txt state_over_time_tol_a.txt
+
+./main3d.gnu.ex inputs_aprox13 ${RUNPARAMS} integrator.rtol_spec=1.0e-5 integrator.rtol_enuc=1.0e-5 integrator.atol_spec=1.0e-5 integrator.atol_enuc=1.0e-5
+mv state_over_time.txt state_over_time_tol_b.txt
+
+./main3d.gnu.ex inputs_aprox13 ${RUNPARAMS} integrator.rtol_spec=1.0e-8 integrator.rtol_enuc=1.0e-8 integrator.atol_spec=1.0e-8 integrator.atol_enuc=1.0e-8
+mv state_over_time.txt state_over_time_tol_c.txt
+
+./main3d.gnu.ex inputs_aprox13 ${RUNPARAMS} integrator.rtol_spec=1.0e-12 integrator.rtol_enuc=1.0e-12 integrator.atol_spec=1.0e-12 integrator.atol_enuc=1.0e-12
+mv state_over_time.txt state_over_time_tol_d.txt
+