work on addressing some issues with the Jacobian scaling (#1479)

this applies if we run with scale_system=1 also as part of this change, allow the Circle theorem to use the numerical Jacobian.
AMReX-Astro · Feb 28, 2024 · adb2600 · adb2600
1 parent 48bf41c
commit adb2600
Show file tree

Hide file tree

Showing 5 changed files with 47 additions and 26 deletions.
diff --git a/integration/integrator_rhs_strang.H b/integration/integrator_rhs_strang.H
@@ -154,10 +154,17 @@ void jac ([[maybe_unused]] const amrex::Real time, BurnT& state, T& int_state, M
     }
 
     // scale the energy derivatives
+
     if (scale_system) {
+        // first the row de/dX
         for (int j = 1; j <= INT_NEQS; ++j) {
             pd(net_ienuc,j) /= state.e_scale;
         }
+
+        // now the column dX/de
+        for (int i = 1; i <= INT_NEQS; ++i) {
+            pd(i,net_ienuc) *= state.e_scale;
+        }
     }
 
     // apply fudge factor:

diff --git a/integration/utils/circle_theorem.H b/integration/utils/circle_theorem.H
@@ -11,6 +11,7 @@
 #include <integrator_rhs_sdc.H>
 #endif
 #include <limits>
+#include <numerical_jacobian.H>
 
 template<typename BurnT, typename T>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
@@ -19,7 +20,14 @@ void circle_theorem_sprad(const amrex::Real time, BurnT& state, T& int_state, am
 
     ArrayUtil::MathArray2D<1, INT_NEQS, 1, INT_NEQS> jac_array;
 
-    jac(time, state, int_state, jac_array);
+    if (integrator_rp::jacobian == 1) {
+        jac(time, state, int_state, jac_array);
+    } else {
+        integrator_to_burn(int_state, state);
+        jac_info_t jac_info;
+        jac_info.h = 0.0_rt;
+        numerical_jac(state, jac_info, jac_array);
+    }
 
     // the Gershgorin circle theorem says that the spectral radius is <
     // max_i ( -a_{ii} + sum_{j,j!=i} |a_{ij}|)

diff --git a/integration/utils/jacobian_utilities.H b/integration/utils/jacobian_utilities.H
@@ -13,9 +13,7 @@ template <typename BurnT>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 Real temperature_to_energy_jacobian (const BurnT& state, const Real& jac_T)
 {
-    Real jac_e = 0.0_rt;
-
-    jac_e = jac_T / state.cv;
+    Real jac_e = jac_T / state.cv;
 
     return jac_e;
 }

diff --git a/integration/utils/numerical_jacobian.H b/integration/utils/numerical_jacobian.H
@@ -12,6 +12,7 @@
 #endif
 #include <integrator_data.H>
 
+
 using namespace integrator_rp;
 
 ///
@@ -204,9 +205,15 @@ void numerical_jac(BurnT& state, const jac_info_t& jac_info, JacNetArray2D& jac)
 
     // scale the energy derivatives
     if (scale_system) {
+        // first the de/dX row
         for (int n = 1; n <= INT_NEQS; ++n) {
             jac(net_ienuc, n) /= state.e_scale;
         }
+
+        // now the dX/de column
+        for (int m = 1; m <= INT_NEQS; ++m) {
+            jac(m, net_ienuc) *= state.e_scale;
+        }
     }
 
     // apply boosting factor:

diff --git a/unit_test/burn_cell/ci-benchmarks/aprox13_RKC_unit_test.out b/unit_test/burn_cell/ci-benchmarks/aprox13_RKC_unit_test.out
@@ -1,4 +1,5 @@
-AMReX (23.06-6-g1e1c433b401c) initialized
+Initializing AMReX (23.08-206-g798ca3356d16)...
+AMReX (23.08-206-g798ca3356d16) initialized
 starting the single zone burn...
 Maximum Time (s): 0.01
 State Density (g/cm^3): 1000000
@@ -32,8 +33,8 @@ RHS at t = 0
   Ni56 3.938787868e-40
 ------------------------------------
 successful? 1
- - Hnuc = 8.558390239e+18
- - added e = 8.558390239e+16
+ - Hnuc = 8.55839024e+18
+ - added e = 8.55839024e+16
  - final T = 1516425860
 ------------------------------------
 e initial = 1.284393683e+17
@@ -43,30 +44,30 @@ new mass fractions:
 He4 0.8760723967
 C12 0.1064099566
 O16 0.0001403204362
-Ne20 8.129701233e-05
-Mg24 0.0002972369574
+Ne20 8.129701234e-05
+Mg24 0.0002972369573
 Si28 0.01113151728
-S32 0.005297014797
-Ar36 0.0005641483161
-Ca40 6.109847027e-06
-Ti44 2.004306621e-09
-Cr48 3.272104798e-14
-Fe52 7.12614193e-20
-Ni56 1.224820021e-26
+S32 0.005297014805
+Ar36 0.0005641483217
+Ca40 6.109847214e-06
+Ti44 2.004306753e-09
+Cr48 3.2721052e-14
+Fe52 7.126143425e-20
+Ni56 1.224820426e-26
 ------------------------------------
 species creation rates: 
 omegadot(He4): -12.39276033
 omegadot(C12): 10.64099566
 omegadot(O16): 0.01403204362
-omegadot(Ne20): 0.008129701233
-omegadot(Mg24): 0.02972369574
+omegadot(Ne20): 0.008129701234
+omegadot(Mg24): 0.02972369573
 omegadot(Si28): 1.113151728
-omegadot(S32): 0.5297014797
-omegadot(Ar36): 0.05641483161
-omegadot(Ca40): 0.0006109847027
-omegadot(Ti44): 2.004306621e-07
-omegadot(Cr48): 3.272104798e-12
-omegadot(Fe52): 7.12614193e-18
-omegadot(Ni56): 1.224720021e-24
+omegadot(S32): 0.5297014805
+omegadot(Ar36): 0.05641483217
+omegadot(Ca40): 0.0006109847214
+omegadot(Ti44): 2.004306753e-07
+omegadot(Cr48): 3.2721052e-12
+omegadot(Fe52): 7.126143424e-18
+omegadot(Ni56): 1.224720426e-24
 number of steps taken: 255
-AMReX (23.06-6-g1e1c433b401c) finalized
+AMReX (23.08-206-g798ca3356d16) finalized