[WIP] add a normalization only after successful VODE step

integration/BackwardEuler/be_integrator.H

@@ -174,7 +174,6 @@ int single_step (BurnT& state, BeT& be, const amrex::Real dt)
 
             // if we didn't set another error, then we probably ran
             // out of iterations, so set nonconvergence
-
             ierr = IERR_CORRECTOR_CONVERGENCE;
 
             // reset the solution to the original

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,39 @@ int dvstep (BurnT& state, DvodeT& vstate)
         DSM = ACNRM / vstate.tq(2);
         if (DSM <= 1.0_rt && valid_update) {
 
+            // do normalization if desired
+
+            if (integrator_rp::do_renormalization_each_step &&
+                !integrator_rp::use_number_densities) {
+
+#if defined(SDC)
+                amrex::Real rhoX_sum{};
+                for (int n = 1; n <= NumSpec; ++n) {
+                    vstate.y(SFS+n) = amrex::Clamp(vstate.y(SFS+n),
+                                                   rho_current * integrator_rp::SMALL_X_SAFE, 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;
+                }
+
+#elif defined(STRANG)
+                amrex::Real X_sum{};
+                for (int n = 1; n <= NumSpec; ++n) {
+                    vstate.y(SFS+n) = amrex::Clamp(vstate.y(SFS+n),
+                                                   integrator_rp::SMALL_X_SAFE, 1.0_rt);
+                    X_sum += vstate.y(SFS+n);
+                }
+
+                amrex::Real fac = 1.0_rt / X_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/_parameters

@@ -97,6 +97,10 @@ do_species_clip              bool         0
 # valid (X > 0) before calling the RHS?
 do_corrector_validation      bool         1
 
+# after a successful step, do we clip the mass fractions / partial
+# densities and renormalize?
+do_renormalization_each_step  bool        0
+
 # flag for turning on the use of number densities for all species
 use_number_densities     bool     0
 

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;