From a397a43531e8d48a93bec811c051a5cefecaf853 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 26 Jun 2024 12:04:32 -0400 Subject: [PATCH 1/7] add a runtime parameter to disable the T/e correction term in Jacobians this is the dT/dX |_e term that appears in all species --- integration/_parameters | 5 +++++ integration/integrator_rhs_sdc.H | 32 ++++++++++++++++++-------------- 2 files changed, 23 insertions(+), 14 deletions(-) diff --git a/integration/_parameters b/integration/_parameters index f3900aaca8..60c00dd76a 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -103,3 +103,8 @@ nse_include_enu_weak bool 1 # for the linear algebra, do we allow pivoting? linalg_do_pivoting bool 1 + +# when converting the Jacobian from (rho, Y, T) to (rho, X, e), do +# we add a correction to the species terms reflecting that +# dT/dX |_e is not zero? +correct_jacobian_for_const_e bool 1 diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index b6272fa5ea..f4a96e7332 100644 --- a/integration/integrator_rhs_sdc.H +++ b/integration/integrator_rhs_sdc.H @@ -154,26 +154,30 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd) // now correct the species derivatives // this constructs dy/dX_k |_e = dy/dX_k |_T - e_{X_k} |_T dy/dT / c_v - eos_re_extra_t eos_state; - eos_state.rho = state.rho; - eos_state.T = state.T; - eos_state.e = state.e; - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = state.xn[n]; - } + if (integrator_rp::correct_jacobian_for_const_e) { + + eos_re_extra_t eos_state; + eos_state.rho = state.rho; + eos_state.T = state.T; + eos_state.e = state.e; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = state.xn[n]; + } #ifdef AUX_THERMO - // make the aux data consistent with the state X's - set_aux_comp_from_X(eos_state); + // make the aux data consistent with the state X's + set_aux_comp_from_X(eos_state); #endif - eos(eos_input_re, eos_state); + eos(eos_input_re, eos_state); - eos_xderivs_t eos_xderivs = composition_derivatives(eos_state); + eos_xderivs_t eos_xderivs = composition_derivatives(eos_state); - for (int m = 1; m <= neqs; m++) { - for (int n = 1; n <= NumSpec; n++) { - pd(m, n) -= eos_xderivs.dedX[n-1] * pd(m, net_ienuc); + for (int m = 1; m <= neqs; m++) { + for (int n = 1; n <= NumSpec; n++) { + pd(m, n) -= eos_xderivs.dedX[n-1] * pd(m, net_ienuc); + } } + } // apply scale_system scaling (if needed) From 626d852f2f5868963c3ea74ae4993b5319334a4a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 26 Jun 2024 15:37:41 -0400 Subject: [PATCH 2/7] add the corrections to Strang --- integration/integrator_rhs_sdc.H | 16 ++++++----- integration/integrator_rhs_strang.H | 41 +++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 7 deletions(-) diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index f4a96e7332..6a7ffd846d 100644 --- a/integration/integrator_rhs_sdc.H +++ b/integration/integrator_rhs_sdc.H @@ -122,15 +122,17 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd) // The Jacobian from the nets is in terms of dYdot/dY, but we want // it was dXdot/dX, so convert here. - for (int n = 1; n <= NumSpec; n++) { - for (int m = 1; m <= neqs; m++) { - pd(n,m) = pd(n,m) * aion[n-1]; + if (!use_number_density) { + for (int n = 1; n <= NumSpec; n++) { + for (int m = 1; m <= neqs; m++) { + pd(n,m) = pd(n,m) * aion[n-1]; + } } - } - for (int m = 1; m <= neqs; m++) { - for (int n = 1; n <= NumSpec; n++) { - pd(m,n) = pd(m,n) * aion_inv[n-1]; + for (int m = 1; m <= neqs; m++) { + for (int n = 1; n <= NumSpec; n++) { + pd(m,n) = pd(m,n) * aion_inv[n-1]; + } } } diff --git a/integration/integrator_rhs_strang.H b/integration/integrator_rhs_strang.H index 7901eca3ba..5b5db6f70a 100644 --- a/integration/integrator_rhs_strang.H +++ b/integration/integrator_rhs_strang.H @@ -153,6 +153,47 @@ void jac ([[maybe_unused]] const amrex::Real time, BurnT& state, T& int_state, M } } + // The system we integrate has the form (rho X_k, rho e) + + // pd is now of the form: + // + // SFS / d(rho X1dot)/dX1 d(rho X1dit)/dX2 ... 1/cv d(rho X1dot)/dT \ // + // | d(rho X2dot)/dX1 d(rho X2dot)/dX2 ... 1/cv d(rho X2dot)/dT | + // SFS-1+nspec | ... | + // SEINT \ d(rho Edot)/dX1 d(rho Edot)/dX2 ... 1/cv d(rho Edot)/dT / + // + // SFS SEINT + + // now correct the species derivatives + // this constructs dy/dX_k |_e = dy/dX_k |_T - e_{X_k} |_T dy/dT / c_v + + if (integrator_rp::correct_jacobian_for_const_e) { + + eos_re_extra_t eos_state; + eos_state.rho = state.rho; + eos_state.T = state.T; + eos_state.e = state.e; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = state.xn[n]; + } +#ifdef AUX_THERMO + // make the aux data consistent with the state X's + set_aux_comp_from_X(eos_state); +#endif + + eos(eos_input_re, eos_state); + + eos_xderivs_t eos_xderivs = composition_derivatives(eos_state); + + for (int m = 1; m <= neqs; m++) { + for (int n = 1; n <= NumSpec; n++) { + pd(m, n) -= eos_xderivs.dedX[n-1] * pd(m, net_ienuc); + } + } + + } + + // scale the energy derivatives if (scale_system) { From 5dc740920a976569fb64820db31466ce465a20a0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 26 Jun 2024 15:45:56 -0400 Subject: [PATCH 3/7] fix compilation --- integration/integrator_rhs_sdc.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index 6a7ffd846d..1625166155 100644 --- a/integration/integrator_rhs_sdc.H +++ b/integration/integrator_rhs_sdc.H @@ -122,7 +122,7 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd) // The Jacobian from the nets is in terms of dYdot/dY, but we want // it was dXdot/dX, so convert here. - if (!use_number_density) { + if (!integrator_rp::use_number_density) { for (int n = 1; n <= NumSpec; n++) { for (int m = 1; m <= neqs; m++) { pd(n,m) = pd(n,m) * aion[n-1]; From da319042d5821b8989004f409b4c7c35a8f2c72f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 26 Jun 2024 16:59:00 -0400 Subject: [PATCH 4/7] fix the fix --- integration/integrator_rhs_sdc.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index 1625166155..01d43ba02f 100644 --- a/integration/integrator_rhs_sdc.H +++ b/integration/integrator_rhs_sdc.H @@ -122,7 +122,7 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd) // The Jacobian from the nets is in terms of dYdot/dY, but we want // it was dXdot/dX, so convert here. - if (!integrator_rp::use_number_density) { + if (!integrator_rp::use_number_densities) { for (int n = 1; n <= NumSpec; n++) { for (int m = 1; m <= neqs; m++) { pd(n,m) = pd(n,m) * aion[n-1]; From 09aa7b8c8021c093b242b88911efe3f07ef0e22f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 28 Jun 2024 11:53:10 -0400 Subject: [PATCH 5/7] update jacobian docs --- sphinx_docs/source/integrators.rst | 111 +++++++++++++++++++++++++---- sphinx_docs/source/sdc.rst | 5 ++ 2 files changed, 103 insertions(+), 13 deletions(-) diff --git a/sphinx_docs/source/integrators.rst b/sphinx_docs/source/integrators.rst index 1e7b5fe547..0f1bcb9355 100644 --- a/sphinx_docs/source/integrators.rst +++ b/sphinx_docs/source/integrators.rst @@ -197,12 +197,16 @@ fields to fill are: * ``state.ydot(net_ienuc)`` : the change in the internal energy from the net, :math:`de/dt` -The righthand side routine is assumed to return the change in *molar fractions*, -:math:`dY_k/dt`. These will be converted to the change in mass fractions, :math:`dX_k/dt` -by the wrappers that call the righthand side routine for the integrator. -If the network builds the RHS in terms of mass fractions directly, :math:`dX_k/dt`, then -these will need to be converted to molar fraction rates for storage, e.g., -:math:`dY_k/dt = A_k^{-1} dX_k/dt`. +.. important:: + + The righthand side routine is assumed to return the change in + *molar fractions*, :math:`dY_k/dt`. These will be converted to the + change in mass fractions, :math:`dX_k/dt` by the wrappers that call + the righthand side routine for the integrator. If the network + builds the RHS in terms of mass fractions directly, + :math:`dX_k/dt`, then these will need to be converted to molar + fraction rates for storage, e.g., :math:`dY_k/dt = A_k^{-1} + dX_k/dt`. Righthand side wrapper ^^^^^^^^^^^^^^^^^^^^^^ @@ -259,15 +263,42 @@ The Jacobian matrix elements are stored in ``jac`` as: * ``jac(net_ienuc, net_ienuc)`` : :math:`d(\dot{e})/de` -The form looks like: +.. important:: + + The Jacobian returned by the network is assumed to be in terms of molar fractions. + However, we do convert the temperature derivative to an energy derivative already + in the network by multipling by $(1/c_v)$. + +The form of the Jacobian return by the integrator looks like: .. math:: \left ( \begin{matrix} - \ddots & \vdots & & \vdots \\ - \cdots & \partial \dot{Y}_m/\partial Y_n & \cdots & \partial \dot{Y}_m/\partial e \\ - & \vdots & \ddots & \vdots \\ - \cdots & \partial \dot{e}/\partial Y_n & \cdots & \partial \dot{e}/\partial e \\ + \dfrac{\partial \dot{Y}_1}{\partial Y_1} & + \dfrac{\partial \dot{Y}_1}{\partial Y_2} & + \cdots & + \dfrac{\partial \dot{Y}_1}{\partial Y_\mathrm{NumSpec}} & + \dfrac{1}{c_v} \dfrac{\partial \dot{Y}_1}{\partial T} \\ + % + \dfrac{\partial \dot{Y}_2}{\partial Y_1} & + \dfrac{\partial \dot{Y}_2}{\partial Y_2} & + \cdots & + \dfrac{\partial \dot{Y}_2}{\partial Y_\mathrm{NumSpec}} & + \dfrac{1}{c_v} \dfrac{\partial \dot{Y}_2}{\partial T} \\ + % + \vdots & \vdots & \ddots & \vdots & \vdots \\ + % + \dfrac{\partial \dot{Y}_\mathrm{NumSpec}}{\partial Y_1} & + \dfrac{\partial \dot{Y}_\mathrm{NumSpec}}{\partial Y_2} & + \cdots & + \dfrac{\partial \dot{Y}_\mathrm{NumSpec}}{\partial Y_\mathrm{NumSpec}} & + \dfrac{1}{c_v} \dfrac{\partial \dot{Y}_\mathrm{NumSpec}}{\partial T} \\ + % + \dfrac{\partial \dot{e}}{\partial Y_1} & + \dfrac{\partial \dot{e}}{\partial Y_2} & + \cdots & + \dfrac{\partial \dot{e}}{\partial Y_\mathrm{NumSpec}} & + \dfrac{1}{c_v} \dfrac{\partial \dot{e}}{\partial T} \end{matrix} \right ) @@ -291,16 +322,70 @@ flow is (for VODE): wrapper above which did the ``clean_state`` and ``update_thermodynamics`` calls. -.. index:: integrator.react_boost +.. index:: integrator.react_boost, integrator.correct_jacobian_for_const_e #. call ``integrator_to_burn()`` to update the ``burn_t`` #. call ``actual_jac()`` to have the network fill the Jacobian array -#. convert the derivative to be mass-fraction-based +#. convert the derivative to be mass-fraction-based. + + Since $Y_k = X_k/A_k$, we have $\partial/\partial X_k = A_k^{-1} \partial/\partial Y_k$. + + We transform by: + + * multiplying all rows of the form $\partial Y_m / \partial \star$ by $A_m$ (where $\star$ is either a molar fraction or $T$/$e$). + + * multiplying all columns of the form $\partial \star / \partial Y_n$ by $1/A_n$. + +#. add correction + terms proportional to :math:`\partial e/\partial X_k |_{\rho, T, X_{j,j\ne k}}` + if ``integrator.correct_jacobian_for_const_e`` is ``1``. + + The system we integrate is $(X_k, e)$, but the derivatives we took in the analytic Jacobian + were in terms of $T$ and not $e$. So we need to correct for the fact that for some quantity + $q$, + + .. math:: + + \left . \frac{\partial q}{\partial X_k} \right |_e \ne \left . \frac{\partial q}{\partial X_k} \right |_T + + If we write $q = q(\rho, T(\rho, X_k, e), X_k)$ then we find that: + + .. math:: + + \left . \frac{\partial q}{\partial X_k} \right |_{\rho, e, X_{j,j\ne k}} = + \left . \frac{\partial q}{\partial X_k} \right |_{\rho, T, X_{j,j\ne k}} - \frac{e_{X_k}}{c_v} \left . \frac{\partial T}{\partial X_k} \right |_{\rho, e, X_{j,j\ne k}} + + where :math:`e_{X_k} = \partial e / \partial X_k |_{\rho, T, X_{j,j\ne k}}`. + + This correction term is described in :cite:`castro_simple_sdc`. #. apply any boosting to the rates if ``integrator.react_boost`` > 0 +The final form of the Jacobian is: + +.. math:: + \left ( + \begin{matrix} + \dfrac{\partial \dot{X}_1}{\partial X_1} - \dfrac{e_{X_1}}{c_v} \dfrac{\partial \dot{X}_1}{\partial T} & + \dfrac{\partial \dot{X}_1}{\partial X_2} - \dfrac{e_{X_2}}{c_v} \dfrac{\partial \dot{X}_1}{\partial T} & + \cdots & + \dfrac{1}{c_v} \dfrac{\partial \dot{X}_1}{\partial T} \\ + % + \dfrac{\partial \dot{X}_2}{\partial X_1} - \dfrac{e_{X_1}}{c_v} \dfrac{\partial \dot{X}_2}{\partial T} & + \dfrac{\partial \dot{X}_2}{\partial X_2} - \dfrac{e_{X_2}}{c_v} \dfrac{\partial \dot{X}_2}{\partial T} & + \cdots & + \dfrac{1}{c_v} \dfrac{\partial \dot{X}_2}{\partial T} \\ + % + \vdots & \vdots & \ddots & \vdots \\ + % + \dfrac{\partial \dot{e}}{\partial X_1} - \dfrac{e_{X_1}}{c_v} \dfrac{\partial \dot{e}}{\partial T} & + \dfrac{\partial \dot{e}}{\partial X_2} - \dfrac{e_{X_2}}{c_v} \dfrac{\partial \dot{e}}{\partial T} & + \cdots & + \dfrac{1}{c_v} \dfrac{\partial \dot{e}}{\partial T} + \end{matrix} + \right ) diff --git a/sphinx_docs/source/sdc.rst b/sphinx_docs/source/sdc.rst index 3cbf8e7c32..b13af3de79 100644 --- a/sphinx_docs/source/sdc.rst +++ b/sphinx_docs/source/sdc.rst @@ -239,3 +239,8 @@ for a reaction network. system was more complicated, and we included density in ${\bf w}$. This is not needed, and we use the Jacobian defined in :cite:`castro_simple_sdc` instead. + +.. note:: + + The final form of the Jacobian is identical to the one used in + Strang integration. From 5683362180d23e5c21d846e852443b4a7e52ad1c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 28 Jun 2024 12:48:20 -0400 Subject: [PATCH 6/7] more fixes --- .github/workflows/burn_cell.yml | 2 +- integration/_parameters | 2 +- sphinx_docs/source/integrators.rst | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/burn_cell.yml b/.github/workflows/burn_cell.yml index 0cd9ddbd17..9004b4dc71 100644 --- a/.github/workflows/burn_cell.yml +++ b/.github/workflows/burn_cell.yml @@ -67,7 +67,7 @@ jobs: - name: Run burn_cell (VODE, ignition_chamulak) run: | cd unit_test/burn_cell - ./main3d.gnu.ex inputs_ignition_chamulak > test.out + ./main3d.gnu.ex inputs_ignition_chamulak integrator.correct_jacobian_for_const_e=0 > test.out - name: Compare to stored output (VODE, ignition_chamulak) run: | diff --git a/integration/_parameters b/integration/_parameters index 60c00dd76a..51426eca66 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -106,5 +106,5 @@ linalg_do_pivoting bool 1 # when converting the Jacobian from (rho, Y, T) to (rho, X, e), do # we add a correction to the species terms reflecting that -# dT/dX |_e is not zero? +# dT/dX at constant e is not zero? correct_jacobian_for_const_e bool 1 diff --git a/sphinx_docs/source/integrators.rst b/sphinx_docs/source/integrators.rst index 0f1bcb9355..ddd146a7e5 100644 --- a/sphinx_docs/source/integrators.rst +++ b/sphinx_docs/source/integrators.rst @@ -267,7 +267,7 @@ The Jacobian matrix elements are stored in ``jac`` as: The Jacobian returned by the network is assumed to be in terms of molar fractions. However, we do convert the temperature derivative to an energy derivative already - in the network by multipling by $(1/c_v)$. + in the network by multiplying by $(1/c_v)$. The form of the Jacobian return by the integrator looks like: From 4ee930200041e6daf1d67932863ac2a6ae1480bd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 31 Jul 2024 14:22:26 -0400 Subject: [PATCH 7/7] sync --- integration/integrator_rhs_sdc.H | 2 +- integration/integrator_rhs_strang.H | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index 271fff9426..98ba3847d5 100644 --- a/integration/integrator_rhs_sdc.H +++ b/integration/integrator_rhs_sdc.H @@ -128,7 +128,7 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd) } } } - + // apply fudge factor: if (integrator_rp::react_boost > 0.0_rt) { diff --git a/integration/integrator_rhs_strang.H b/integration/integrator_rhs_strang.H index 8e554e54d4..1b5a3be83f 100644 --- a/integration/integrator_rhs_strang.H +++ b/integration/integrator_rhs_strang.H @@ -176,9 +176,9 @@ void jac ([[maybe_unused]] const amrex::Real time, BurnT& state, T& int_state, M eos_xderivs_t eos_xderivs = composition_derivatives(eos_state); - for (int m = 1; m <= neqs; m++) { - for (int n = 1; n <= NumSpec; n++) { - pd(m, n) -= eos_xderivs.dedX[n-1] * pd(m, net_ienuc); + for (int jcol = 1; jcol <= NumSpec;++jcol) { + for (int irow = 1; irow <= neqs; ++irow) { + pd(irow, jcol) -= eos_xderivs.dedX[jcol-1] * pd(irow, net_ienuc); } }