[WIP] implement scale_system for simplified-SDC

integration/BackwardEuler/actual_integrator.H

@@ -43,8 +43,9 @@ void actual_integrator (BurnT& state, const Real dt)
 
     eos(eos_input_rt, state);
 
-    // set the scaling for energy if we integrate it dimensionlessly
+    // set the scaling for density / energy if we integrate it dimensionlessly
     state.e_scale = state.e;
+    state.rho_scale = state.rho;
 
     if (scale_system) {
         // the absolute tol for energy needs to reflect the scaled

integration/RKC/actual_integrator.H

@@ -47,8 +47,9 @@ void actual_integrator (BurnT& state, Real dt)
 
     eos(eos_input_rt, state);
 
-    // set the scaling for energy if we integrate it dimensionlessly
+    // set the scaling for energy and density if we integrate it dimensionlessly
     state.e_scale = state.e;
+    scale.rho_scale = state.rho;
 
     if (scale_system) {
         // the absolute tol for energy needs to reflect the scaled

integration/VODE/actual_integrator.H

@@ -59,8 +59,9 @@ void actual_integrator (BurnT& state, Real dt)
 
     eos(eos_input_rt, state);
 
-    // set the scaling for energy if we integrate it dimensionlessly
+    // set the scaling for energy and density if we integrate it dimensionlessly
     state.e_scale = state.e;
+    state.rho_scale = state.rho;
 
     if (scale_system) {
         // the absolute tol for energy needs to reflect the scaled

integration/VODE/actual_integrator_simplified_sdc.H

@@ -82,6 +82,12 @@ void actual_integrator (BurnT& state, Real dt)
     vode_state.atol_spec = sdc_min_density * atol_spec * sdc_tol_fac;
     vode_state.rtol_spec = rtol_spec * sdc_tol_fac;
 
+    // if we are scaling the system, then update the atols accordingly
+    if (scale_system) {
+        vode_state.atol_enuc /= (state.rho_scale * state.e_scale);
+        vode_state.atol_spec /= state.rho_scale;
+    }
+
     // Call the integration routine.
 
     int istate = dvode(state, vode_state);

integration/integrator_rhs_simplified_sdc.H

@@ -74,7 +74,7 @@ void rhs(const Real time, BurnT& state, T& int_state, RArray1D& ydot,
     }
 
     // convert back to the form needed by the integrator -- this will
-    // add the advective terms
+    // add the advective terms and apply any scaling via scale_system
 
     rhs_to_int(time, state, ydot);
 
@@ -176,6 +176,27 @@ void jac (const Real time, BurnT& state, T& int_state, MatrixType& pd)
         }
     }
 
+    // apply scaling via scale_system
+
+    if (scale_system) {
+
+        // scale the derivatives of energy wrt species, d(edot) / dX ...
+
+        for (int j = 1; j <= NumSpec; ++j) {
+            pd(net_ienuc, j) /= state.e_scale;
+        }
+
+        // scale the derivatives of species wrt energy
+
+        for (int i = 1; i <= NumSpec; ++i) {
+            pd(i, net_ienuc) *= state.e_scale;
+        }
+
+
+        // d(edot) / de is unscaled
+
+    }
+
 }
 
 #endif

integration/integrator_type.H

@@ -18,7 +18,7 @@ void update_density_in_time(const Real time, BurnT& state)
     //
     // we are always integrating from t = 0, so there is no offset
     // time needed here.  The indexing of ydot_a is based on
-    // the indices in burn_t and is 0-based    
+    // the indices in burn_t and is 0-based
     state.y[SRHO] = amrex::max(state.rho_orig + state.ydot_a[SRHO] * time, EOSData::mindens);
 
     // for consistency

integration/integrator_type_simplified_sdc.H

@@ -16,10 +16,13 @@ void clean_state(const Real time, BurnT& state, T& int_state)
     // Ensure that mass fractions always stay positive.
 
     if (do_species_clip) {
+
+        // if we are scaling the system, then the species are (rho / rho_scale) * X_k
+        Real rho_ref = scale_system ? state.rho / state.rho_scale : state.rho;
         for (int n = 1; n <= NumSpec; ++n) {
             // we use 1-based indexing, so we need to offset SFS
-            int_state.y(SFS+n) = amrex::max(amrex::min(int_state.y(SFS+n), state.rho),
-                                            state.rho * SMALL_X_SAFE);
+            int_state.y(SFS+n) = amrex::max(amrex::min(int_state.y(SFS+n), rho_ref),
+                                            rho_ref * SMALL_X_SAFE);
         }
     }
 
@@ -38,7 +41,13 @@ void clean_state(const Real time, BurnT& state, T& int_state)
             // use 1-based indexing
             nspec_sum += int_state.y(SFS+n);
         }
+
+        // if we are doing scale_system, then we summed (rho / rho_scale) X_k,
+        // so the normalization is rho_scale / rho
         nspec_sum /= state.y[SRHO];
+        if (scale_system) {
+            nspec_sum *= state.rho_scale;
+        }
 
         for (int n = 1; n <= NumSpec; n++) {
             int_state.y(SFS+n) /= nspec_sum;
@@ -62,9 +71,21 @@ void burn_to_int(BurnT& state, T& int_state)
 
     // store the original rho and rho e
     state.rho_orig = state.y[SRHO];
-
     state.rhoe_orig = state.y[SEINT];
 
+    // normalize
+    state.e_scale = state.rhoe_orig / state.rho_orig;
+    state.rho_scale = state.rho_orig;
+
+    if (scale_system) {
+        // we are integrating (rho X, rho e), so normalize by the
+        // initial scale
+        int_state.y(net_ienuc) /= (state.rho_scale * state.e_scale);
+        for (int n = 1; n <= NumSpec; ++n) {
+            int_state.y(n) /= state.rho_scale;
+        }
+    }
+
 }
 
 
@@ -82,6 +103,14 @@ void int_to_burn(const Real time, const T& int_state, BurnT& state)
         state.y[n] = int_state.y(n+1);
     }
 
+    // undo any scaling
+    if (scale_system) {
+        state.y[SEINT] *= state.rho_scale * state.e_scale;
+        for (int n = 0; n < NumSpec; ++n) {
+            state.y[SFS+n] *= state.rho_scale;
+        }
+    }
+
     // update rho in the burn_t state
     // this may be redundant, but better to be safe
 
@@ -90,8 +119,7 @@ void int_to_burn(const Real time, const T& int_state, BurnT& state)
     Real rhoInv = 1.0_rt / state.rho;
 
     for (int n = 0; n < NumSpec; n++) {
-        // int_state uses 1-based indexing
-        state.xn[n] = int_state.y(SFS+1+n) * rhoInv;
+        state.xn[n] = state.y[SFS+n] * rhoInv;
     }
 
 #ifdef AUX_THERMO
@@ -106,7 +134,7 @@ void int_to_burn(const Real time, const T& int_state, BurnT& state)
 
     // set internal energy for EOS call
 
-    state.e = int_state.y(SEINT+1) * rhoInv;
+    state.e = state.y[SEINT] * rhoInv;
 
     state.T = state.T;  // initial guess
 
@@ -134,7 +162,7 @@ void rhs_to_int([[maybe_unused]] const Real time,
     // reaction network.  Note that these are in terms of dY/dt
     // we will convert these in place to the form needed for SDC
 
-    // convert from dY/dt to dX/dt. The species derivatives are the
+    // convert from dY/dt to rho dX/dt. The species derivatives are the
     // first NumSpec components of ydot coming from the reaction network
 
     BL_ASSERT(SFS == 0);
@@ -159,6 +187,15 @@ void rhs_to_int([[maybe_unused]] const Real time,
         ydot(n+1) += state.ydot_a[n];
     }
 
+    // finally scale the ydots if we are doing scale_system
+
+    if (scale_system) {
+        for (int n = 1; n <= NumSpec; ++n) {
+            ydot(n) /=  state.rho_scale;
+        }
+        ydot(SEINT+1) /= (state.rho_scale * state.e_scale);
+    }
+
 }
 
 #endif

interfaces/burn_type.H

@@ -77,8 +77,9 @@ struct burn_t
   Real aux[NumAux];
 #endif
 
-  // scaling used to make e we integrate dimensionless
+  // scaling used to make e or (rho e) we integrate dimensionless
   Real e_scale;
+  Real rho_scale;
 
   // derivatives needed for calling the EOS
   Real dedr;