add a normalization only after successful VODE step

Author: zingale

integration/VODE/vode_dvstep.H

@@ -223,12 +223,18 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
 #ifdef SDC
         if (vstate.y(SEINT+1) < 0.0_rt) {
+#ifdef MICROPHYSICS_DEBUG
+            std::cout << "rejecting due to rho e" << std::endl;
+#endif
             valid_update = false;
         }
 
         amrex::Real rho_current = state.rho_orig + vstate.tn * state.ydot_a[SRHO];
 
         if (rho_current < 0.0_rt) {
+#ifdef MICROPHYSICS_DEBUG
+            std::cout << "rejecting due to rho" << std::endl;
+#endif
             valid_update = false;
         }
 #endif
@@ -294,11 +300,17 @@ int dvstep (BurnT& state, DvodeT& vstate)
             }
 
             if (vstate.y(SFS+i) < -species_failure_tolerance * rho_current) {
+#ifdef MICROPHYSICS_DEBUG
+                std::cout << "rejecting step due to negative species " << i << std::endl;
+#endif
                 valid_update = false;
                 break;
             }
 
             if (vstate.y(SFS+i) > (1.0_rt + species_failure_tolerance) * rho_current) {
+#ifdef MICROPHYSICS_DEBUG
+                std::cout << "rejecting step due to large species " << i << std::endl;
+#endif
                 valid_update = false;
                 break;
             }
@@ -313,6 +325,19 @@ int dvstep (BurnT& state, DvodeT& vstate)
         DSM = ACNRM / vstate.tq(2);
         if (DSM <= 1.0_rt && valid_update) {
 
+#ifdef SDC
+            // do normalization if desired
+            amrex::Real rhoX_sum{};
+            for (int n = 1; n <= NumSpec; ++n) {
+                vstate.y(SFS+n) = amrex::Clamp(vstate.y(SFS+n), 0.0, rho_current);
+                rhoX_sum += vstate.y(SFS+n);
+            }
+
+            amrex::Real fac = rho_current / rhoX_sum;
+            for (int n = 1; n <= NumSpec; ++n) {
+                vstate.y(SFS+n) *= fac;
+            }
+#endif
             // After a successful step, update the yh and TAU arrays and decrement
             // NQWAIT. If NQWAIT is then 1 and NQ .lt. MAXORD, then ACOR is saved
             // for use in a possible order increase on the next step.

integration/integrator_setup_sdc.H

@@ -19,9 +19,9 @@ struct state_backup_t {
     amrex::Real T_in{};
     amrex::Real rhoe_in{};
 #ifndef AMREX_USE_GPU
-    amrex::Real xn_in[NumSpec]{};
+    amrex::Real rhoxn_in[NumSpec]{};
 #ifdef AUX_THERMO
-    amrex::Real aux_in[NumAux]{};
+    amrex::Real rhoaux_in[NumAux]{};
 #endif
 #endif
 };
@@ -115,11 +115,11 @@ state_backup_t integrator_backup (const BurnT& state) {
 
 #ifndef AMREX_USE_GPU
     for (int n = 0; n < NumSpec; ++n) {
-        state_save.xn_in[n] = state.y[SFS+n] / state.y[SRHO];
+        state_save.rhoxn_in[n] = state.y[SFS+n];
     }
 #ifdef AUX_THERMO
     for (int n = 0; n < NumAux; ++n) {
-        state_save.aux_in[n] = state.y[SFX+n] / state.y[SRHO];
+        state_save.rhoaux_in[n] = state.y[SFX+n];
     }
 #endif
 #endif
@@ -223,14 +223,14 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
             std::cout << "dens start = " << std::setprecision(OUTDIGITS) << state.rho_orig << std::endl;
             std::cout << "temp start = " << std::setprecision(OUTDIGITS) << state_save.T_in << std::endl;
             std::cout << "rhoe start = " << std::setprecision(OUTDIGITS) << state_save.rhoe_in << std::endl;
-            std::cout << "xn start = ";
-            for (const auto X : state_save.xn_in) {
+            std::cout << "rho xn start = ";
+            for (const auto X : state_save.rhoxn_in) {
                 std::cout << std::setprecision(OUTDIGITS) << X << " ";
             }
             std::cout << std::endl;
 #ifdef AUX_THERMO
-            std::cout << "aux start = ";
-            for (const auto aux : state_save.aux_in) {
+            std::cout << "rho aux start = ";
+            for (const auto aux : state_save.rhoaux_in) {
                 std::cout << std::setprecision(OUTDIGITS) << aux << " ";
             }
             std::cout << std::endl;