Simple sdc analytic jac memory (#955)

this is similar to the change we did for the jacobian = 3 version in #917
AMReX-Astro · Jul 31, 2022 · 593f02b · 593f02b
1 parent a495221
commit 593f02b
Show file tree

Hide file tree

Showing 4 changed files with 13 additions and 259 deletions.
diff --git a/integration/VODE/Make.package b/integration/VODE/Make.package
@@ -2,8 +2,6 @@ ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
   CEXE_headers += actual_integrator_simplified_sdc.H
   CEXE_headers += vode_rhs_simplified_sdc.H
   CEXE_headers += vode_type_simplified_sdc.H
-  CEXE_headers += vode_jac_util_simplified_sdc.H
-  CEXE_sources += vode_jac_util_simplified_sdc.cpp
 else
   ifneq ($(USE_TRUE_SDC), TRUE)
     CEXE_headers += vode_rhs.H

diff --git a/integration/VODE/vode_jac_util_simplified_sdc.H b/integration/VODE/vode_jac_util_simplified_sdc.H
diff --git a/integration/VODE/vode_jac_util_simplified_sdc.cpp b/integration/VODE/vode_jac_util_simplified_sdc.cpp
diff --git a/integration/VODE/vode_rhs_simplified_sdc.H b/integration/VODE/vode_rhs_simplified_sdc.H
@@ -6,7 +6,6 @@
 #include <extern_parameters.H>
 
 #include <vode_type_simplified_sdc.H>
-#include <vode_jac_util_simplified_sdc.H>
 #include <actual_network.H>
 #include <actual_rhs.H>
 #include <numerical_jacobian.H>
@@ -143,36 +142,40 @@ void jac (burn_t& state, dvode_t& vode_state, MatrixType& pd)
 
         // Call the specific network routine to get the Jacobian.
 
-        JacNetArray2D jac_react;
-
-        actual_jac(state, jac_react);
+        actual_jac(state, 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++) {
-                jac_react(n,m) = jac_react(n,m) * aion[n-1];
+                pd(n,m) = pd(n,m) * aion[n-1];
             }
         }
 
         for (int m = 1; m <= neqs; m++) {
             for (int n = 1; n <= NumSpec; n++) {
-                jac_react(m,n) = jac_react(m,n) * aion_inv[n-1];
+                pd(m,n) = pd(m,n) * aion_inv[n-1];
             }
         }
 
         // apply fudge factor:
 
         if (react_boost > 0.0_rt) {
-            jac_react.mul(react_boost);
+            pd.mul(react_boost);
         }
 
-        // convert to the system we are using
+        // The system we integrate has the form (rho X_k, rho e)
 
-        jac_to_vode(vode_state.tn, state, jac_react, vode_state, pd);
+        // 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
 
     }
 
-
 }
 #endif