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;