more true SDC fixes

AMReX-Astro · Sep 20, 2023 · fe4fc47 · fe4fc47
1 parent 9f96ca5
commit fe4fc47
Show file tree

Hide file tree

Showing 3 changed files with 3 additions and 64 deletions.
diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H
@@ -44,24 +44,6 @@ void dvhin (BurnT& state, DvodeT& vstate, Real& H0, int& NITER, int& IER)
     // Set an upper bound on h based on vstate.tout-vstate.t and the initial Y and YDOT.
     Real HUB = PT1 * TDIST;
 
-#ifdef TRUE_SDC
-    for (int i = 1; i <= int_neqs; ++i) {
-        Real atol;
-        if (i == 1) {
-            atol = vstate.atol_dens;
-        } else if (i == NumSpec+2) {
-            atol = vstate.atol_enuc;
-        } else {
-            atol = vstate.atol_spec;
-        }
-        Real DELYI = PT1 * std::abs(vstate.yh(i,1)) + atol;
-        Real AFI = std::abs(vstate.yh(i,2));
-        if (AFI * HUB > DELYI) {
-            HUB = DELYI / AFI;
-        }
-    }
-
-#else
     for (int i = 1; i <= int_neqs; ++i) {
         Real atol = i <= NumSpec ? vstate.atol_spec : vstate.atol_enuc;
         Real DELYI = PT1 * std::abs(vstate.yh(i,1)) + atol;
@@ -70,7 +52,6 @@ void dvhin (BurnT& state, DvodeT& vstate, Real& H0, int& NITER, int& IER)
             HUB = DELYI / AFI;
         }
     }
-#endif
 
     // Set initial guess for h as geometric mean of upper and lower bounds.
     int iter = 0;

diff --git a/integration/VODE/vode_dvode.H b/integration/VODE/vode_dvode.H
@@ -72,24 +72,12 @@ int dvode (BurnT& state, DvodeT& vstate)
     vstate.NQ = 1;
     vstate.H = 1.0_rt;
 
-#ifdef TRUE_SDC
-    vstate.ewt(1) = vstate.rtol_dens * std::abs(vstate.yh(1,1)) + vstate.atol_dens;
-    vstate.ewt(1) = 1.0_rt / vstate.ewt(1);
-
-    for (int i = 2; i <= NumSpec+1; ++i) {
-        vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec;
-        vstate.ewt(i) = 1.0_rt / vstate.ewt(i);
-    }
-    vstate.ewt(NumSpec+2) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+2,1)) + vstate.atol_enuc;
-    vstate.ewt(NumSpec+2) = 1.0_rt / vstate.ewt(NumSpec+2);
-#else
     for (int i = 1; i <= NumSpec; ++i) {
         vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec;
         vstate.ewt(i) = 1.0_rt / vstate.ewt(i);
     }
     vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc;
     vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1);
-#endif
 
     // Call DVHIN to set initial step size H0 to be attempted.
     H0 = 0.0_rt;
@@ -157,24 +145,12 @@ int dvode (BurnT& state, DvodeT& vstate)
 
            }
 
-#ifdef TRUE_SDC
-           vstate.ewt(1) = vstate.rtol_dens * std::abs(vstate.yh(1,1)) + vstate.atol_dens;
-           vstate.ewt(1) = 1.0_rt / vstate.ewt(1);
-
-           for (int i = 2; i <= NumSpec+1; ++i) {
-               vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec;
-               vstate.ewt(i) = 1.0_rt / vstate.ewt(i);
-           }
-           vstate.ewt(NumSpec+2) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+2,1)) + vstate.atol_enuc;
-           vstate.ewt(NumSpec+2) = 1.0_rt / vstate.ewt(NumSpec+2);
-#else
            for (int i = 1; i <= NumSpec; ++i) {
                vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec;
                vstate.ewt(i) = 1.0_rt / vstate.ewt(i);
            }
            vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc;
            vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1);
-#endif
 
        }
        else {

diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H
@@ -149,7 +149,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
         y_save(i) = vstate.y(i);
     }
 #endif
-#ifdef SIMPLIFIED_SDC
+#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
     Real rho_old = state.rho_orig + TOLD * state.ydot_a[SRHO];
     for (int i = 1; i <= int_neqs; ++i) {
         y_save(i) = vstate.y(i);
@@ -224,7 +224,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
         bool valid_update = true;
 
-#ifdef SIMPLIFIED_SDC
+#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
 	if (vstate.y(SEINT+1) < 0.0_rt) {
 	    valid_update = false;
 	}
@@ -271,7 +271,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
 #endif
 
-#ifdef SIMPLIFIED_SDC
+#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
 
             // these are basically the same checks as with Strang
             // above, except now we are evolving rhoX instead of X.
@@ -303,24 +303,6 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
 #endif
 
-#ifdef TRUE_SDC
-            // as above, we want to make sure that the mass fractions stay bounded, but we are
-            // evolving:
-            //
-            //  y(1) = density
-            //  y(2:2+NumSpec-1) = rho X
-            //  y(NumSpec+2) = rho e
-
-            if (vstate.y(1+i) < -species_failure_tolerance * vstate.y(1)) {
-                valid_update = false;
-            }
-
-            if (vstate.y(1+i) > (1.0_rt + species_failure_tolerance) * vstate.y(1)) {
-                valid_update = false;
-            }
-#endif
-
-
         }
 
         // The corrector has converged (NFLAG = 0). The local error test is