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