From 5e1de1e96d6dd500c5cd6fc8275475d26b3a1579 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 14 Jan 2026 14:19:49 -0500 Subject: [PATCH 1/2] add a normalization only after successful VODE step --- integration/VODE/vode_dvstep.H | 25 +++++++++++++++++++++++++ integration/integrator_setup_sdc.H | 16 ++++++++-------- 2 files changed, 33 insertions(+), 8 deletions(-) diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index 049c68c51..c0f319758 100644 --- a/integration/VODE/vode_dvstep.H +++ b/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. diff --git a/integration/integrator_setup_sdc.H b/integration/integrator_setup_sdc.H index a045b1c37..eb192b737 100644 --- a/integration/integrator_setup_sdc.H +++ b/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; From b2f4f7898088ab3cff337c845df4213e0e0d4c79 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 18 Jan 2026 19:07:19 -0500 Subject: [PATCH 2/2] update normalization --- integration/BackwardEuler/be_integrator.H | 1 - integration/VODE/vode_dvstep.H | 40 +++++++++++++++++------ integration/_parameters | 4 +++ 3 files changed, 34 insertions(+), 11 deletions(-) diff --git a/integration/BackwardEuler/be_integrator.H b/integration/BackwardEuler/be_integrator.H index d50422b37..8454b0c01 100644 --- a/integration/BackwardEuler/be_integrator.H +++ b/integration/BackwardEuler/be_integrator.H @@ -153,7 +153,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 diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index c0f319758..bad326b19 100644 --- a/integration/VODE/vode_dvstep.H +++ b/integration/VODE/vode_dvstep.H @@ -325,19 +325,39 @@ 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; - } + 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. diff --git a/integration/_parameters b/integration/_parameters index f6fb6a41f..fad1fe451 100644 --- a/integration/_parameters +++ b/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