From ffae4da80aff8b391b0b934151d4723388ba2fb5 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 21 Sep 2023 08:19:44 -0400 Subject: [PATCH] update true-SDC to be in sync with upcoming Castro cleaning (#1338) this will make it more in line with how simplified-SDC works --- integration/VODE/vode_dvhin.H | 19 ------------------- integration/VODE/vode_dvode.H | 24 ------------------------ integration/VODE/vode_dvstep.H | 24 +++--------------------- integration/integrator_data.H | 4 ++-- interfaces/burn_type.H | 6 ------ 5 files changed, 5 insertions(+), 72 deletions(-) diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H index 74700e64dd..1ef7779926 100644 --- 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 index 9466d2676b..924c010325 100644 --- 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 index db28b936f1..5b850225eb 100644 --- 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 diff --git a/integration/integrator_data.H b/integration/integrator_data.H index 62330ff03e..5f64ca7e07 100644 --- a/integration/integrator_data.H +++ b/integration/integrator_data.H @@ -14,7 +14,7 @@ const int INT_NEQS = SVAR_EVOLVE; #endif #ifdef TRUE_SDC -const int INT_NEQS = NumSpec + 2; +const int INT_NEQS = NumSpec + 1; #endif @@ -56,7 +56,7 @@ constexpr int integrator_neqs () #endif #ifdef TRUE_SDC - int_neqs = NumSpec + 2; + int_neqs = NumSpec + 1; #endif return int_neqs; diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H index 77d14f3231..9e7edbd4df 100644 --- a/interfaces/burn_type.H +++ b/interfaces/burn_type.H @@ -131,12 +131,6 @@ struct burn_t // dx is useful for estimating timescales for equilibriation Real dx; -#ifdef TRUE_SDC - Real E_var; - Real mom[3]; - Real f_source[NumSpec+2]; -#endif - Real mu; Real mu_e; Real y_e;