@@ -52,8 +52,9 @@ void rk4_step(amrex::Real t, amrex::Real dt,
5252 const amrex::Real (&f)[NEQ], const State& state,
5353 void (*rhs)(amrex::Real, const amrex::Real (&)[NEQ],
5454 const State&, amrex::Real (&)[NEQ]),
55- amrex::Real (&fout)[NEQ])
56- {
55+ amrex::Real (&fout)[NEQ]) {
56+ // Single rk4 update.
57+
5758 amrex::Real k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ];
5859 amrex::Real ftmp[NEQ];
5960
@@ -83,8 +84,9 @@ void rk4_adaptive(amrex::Real t0, amrex::Real tmax,
8384 const amrex::Real (&f0)[NEQ], const State& state,
8485 void (*rhs)(amrex::Real, const amrex::Real (&)[NEQ],
8586 const State&, amrex::Real (&)[NEQ]),
86- amrex::Real (&f_new)[NEQ])
87- {
87+ amrex::Real (&f_new)[NEQ]) {
88+ // rk4 with adaptive timestepping
89+
8890 constexpr amrex::Real SAFETY = 0 .9_rt;
8991 constexpr amrex::Real GROW = -0 .2_rt; // 1/(order+1)
9092 constexpr amrex::Real SHRINK = -0 .25_rt; // 1/order
@@ -186,8 +188,7 @@ void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real dt_init
186188 f_new[n] = f0[n];
187189 }
188190
189- amrex::Real k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ];
190- amrex::Real ftmp[NEQ], f[NEQ];
191+ amrex::Real f[NEQ];
191192
192193 while (t < tmax) {
193194
@@ -201,27 +202,7 @@ void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real dt_init
201202 f[n] = f_new[n];
202203 }
203204
204- rhs (t, f, state, k1);
205-
206- for (int n = 0 ; n < NEQ; ++n) {
207- ftmp[n] = f[n] + 0 .5_rt * dt * k1[n];
208- }
209- rhs (t + 0 .5_rt * dt, ftmp, state, k2);
210-
211- for (int n = 0 ; n < NEQ; ++n) {
212- ftmp[n] = f[n] + 0 .5_rt * dt * k2[n];
213- }
214- rhs (t + 0 .5_rt * dt, ftmp, state, k3);
215-
216- for (int n = 0 ; n < NEQ; ++n) {
217- ftmp[n] = f[n] + dt * k3[n];
218- }
219- rhs (t + dt, ftmp, state, k4);
220-
221- for (int n = 0 ; n < NEQ; ++n) {
222- f_new[n] = f[n] + (dt / 6 .0_rt) *
223- (k1[n] + 2 .0_rt*k2[n] + 2 .0_rt*k3[n] + k4[n]);
224- }
205+ rk4_step (t, dt, f, state, rhs, f_new);
225206
226207 // update time
227208 t += dt;
@@ -431,8 +412,9 @@ amrex::Real heating_term(const amrex::Real T0, const eos_t& eos_state) {
431412 // amrex::Real dFscn_dT = scor * dscor2_dt + dscor_dt * scor2;
432413
433414 // eps_3a = 5.3e21 Fscn rho5^2 Y^3 / T8^3 exp(-44 / T8)
434- amrex::Real heat = 5.3e11 * Fscn * std::pow (tmp_state.rho , 2 ) *
435- std::pow (tmp_state.xn [Species::He4-1 ], 3 ) * std::pow (tmp_state.T * 1.0e-8 , -3 ) *
415+ amrex::Real heat = 5.3e11 * Fscn * tmp_state.rho * tmp_state.rho *
416+ amrex::Math::powi<3 >(tmp_state.xn [Species::He4-1 ]) *
417+ amrex::Math::powi<-3 >(tmp_state.T * 1.0e-8 ) *
436418 std::exp (-44.0 * 1.0e8 / tmp_state.T );
437419
438420 // If there is enough hydrogen for the rapid proton capture reactions:
@@ -517,11 +499,11 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) {
517499 // If we were to compute kappa manually, we need eos.
518500 // i.e. we need two inputs. We have p=g * yt, so we need rho?
519501
520- state.Tt = std::pow ( 2650.0 * state.Ft * state.yt , 0.25 );
502+ state.Tt = std::sqrt ( std::sqrt ( 2650.0 * state.Ft * state.yt ) );
521503
522504 // Now set the lower boundary conditions: yb, Tb, and Fb
523505 // Here Fb is assumed to be given. Note that the original version
524- // includes a special term for compressional heating, if thats included
506+ // includes a special term for compressional heating, if that is included
525507 // then we have to solve for Fb. However, it is sufficient to include
526508 // an extra constant heating term to account for compressional heating.
527509 // yb is what we're after, so we need Tb.
@@ -539,8 +521,13 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) {
539521 amrex::Real f_array[NEQ];
540522 amrex::Real f0_array[NEQ] = {state.Tt , state.Ft , 0 .0_rt};
541523
542- // t0=yt, tmax=yb, step_size =yt, f0={Tt, Ft, 0.0}
524+ // t0=yt, tmax=yb, dt =yt, f0={Tt, Ft, 0.0}
543525 // Integrates to find f_array corresponding to yb.
526+
527+ // rk4<NEQ>(state.yt, state.yb, state.yt,
528+ // f0_array, state, derivs,
529+ // f_array);
530+
544531 rk4_adaptive<NEQ>(state.yt , state.yb , state.yt , 1.0e-6 ,
545532 f0_array, state, derivs,
546533 f_array);
0 commit comments