diff --git a/integration/_parameters b/integration/_parameters index db40f53ae0..6bae04c821 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -80,3 +80,11 @@ scale_system integer 0 # for the NSE update predictor-corrector, how many iterations # do we take to get the new time NSE state nse_iters integer 3 + +# for SDC+NSE, when estimating the derivatives of the NSE table +# quantities, what fraction of dt do we use for the finite-difference +# estimate +nse_deriv_dt_factor real 0.05 + +# for NSE update, do we include the weak rate neutrino losses? +nse_include_enu_weak integer 1 diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 9c8345e8d2..4ba3132f56 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -14,10 +14,9 @@ #include #include -#include +#include #ifdef NSE_TABLE -#include #include #endif #ifdef NSE_NET @@ -27,45 +26,33 @@ using namespace amrex; #if defined(NSE_TABLE) + /// -/// update the state due to NSE changes for the simplified-SDC case -/// this version works with the tabulated NSE and requires AUX_THERMO +/// this acts as an explicit Euler step for the system (rho e, rho aux) +/// on input, *_source are the reactive sources at time t0 and on output +/// they are the sources at time t0+dt /// -template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void sdc_nse_burn(BurnT& state, const Real dt) { - - using namespace AuxZero; - - state.success = true; - state.n_rhs = 0; - state.n_jac = 0; - - // we need the initial Ye to get the NSE state. We also - // want the initial rho since we may not be entering - // this routine already in NSE. This means that there will - // be an energy release from us instantaneously adjusting - // into NSE (the first call to nse_interp) + an energy - // release from the evolution of NSE over the timestep - - Real ye_in = state.y[SFX+iye] / state.rho; - Real rho_bea_old = state.y[SFX+ibea]; - - // density and momentum have no reactive sources - Real rho_old = state.y[SRHO]; - - state.y[SRHO] += dt * state.ydot_a[SRHO]; - state.y[SMX] += dt * state.ydot_a[SMX]; - state.y[SMY] += dt * state.ydot_a[SMY]; - state.y[SMZ] += dt * state.ydot_a[SMZ]; - - // if we are doing drive_initial_convection, we want to use - // the temperature that comes in through T_fixed - - Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T; +void nse_derivs(const Real rho0, const Real rhoe0, const Real *rhoaux0, + const Real dt, const Real *ydot_a, + Real& drhoedt, Real* drhoauxdt, const Real T_fixed) { + + // start with the current state and compute T + + Real T0; + Real abar; + Real Ye0 = rhoaux0[iye] / rho0; + + if (T_fixed > 0) { + T0 = T_fixed; + abar = rhoaux0[iabar] / rho0; + } else { + Real e0 = rhoe0 / rho0; + T0 = 1.e8; // initial guess + nse_T_abar_from_e(rho0, e0, Ye0, T0, abar); + } - // get the neutrino loss term -- we want to use the state that we - // came in here with, so the original Abar and Zbar + // compute the plasma neutrino losses at t0 Real snu{0.0}; Real dsnudt{0.0}; @@ -73,117 +60,202 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real dsnuda{0.0}; Real dsnudz{0.0}; - Real abar = state.y[SFX+iabar] / rho_old; - Real zbar = abar * ye_in; + Real zbar = abar * Ye0; #ifdef NEUTRINOS constexpr int do_derivatives = 0; - sneut5(T_in, rho_old, abar, zbar, + sneut5(T0, rho0, abar, zbar, snu, dsnudt, dsnudd, dsnuda, dsnudz); #endif - Real snu_old = snu; - // get the current NSE state from the table -- this will be used - // to compute the NSE evolution sources + // call the NSE table at t0 + + constexpr bool skip_X_fill{true}; nse_table_t nse_state; - nse_state.T = T_in; - nse_state.rho = state.rho; - nse_state.Ye = ye_in; + nse_state.T = T0; + nse_state.rho = rho0; + nse_state.Ye = Ye0; + nse_interp(nse_state, skip_X_fill); + + Real abar0_out = nse_state.abar; + Real bea0_out = nse_state.bea; + Real dyedt0 = nse_state.dyedt; + + // construct initial sources at t0 + + Real rhoe_source{}; + if (integrator_rp::nse_include_enu_weak == 1) { + rhoe_source = -rho0 * (nse_state.e_nu + snu); + } else { + rhoe_source = -rho0 * snu; + } - nse_interp(nse_state); + Real rhoaux_source[NumAux]; + rhoaux_source[iye] = rho0 * dyedt0; + rhoaux_source[iabar] = 0.0; + rhoaux_source[ibea] = rho0 * nse_state.dbeadt; - Real dyedt_old = nse_state.dyedt; + // evolve for eps * dt; + Real tau = integrator_rp::nse_deriv_dt_factor * dt; - // predict the U^{n+1,*} state with only estimates of the aux - // reaction terms and advection to dt + Real rho1 = rho0 + tau * ydot_a[SRHO]; + Real rhoe1 = rhoe0 + tau * ydot_a[SEINT] + tau * rhoe_source; + Real rhoaux1[NumAux]; + for (int n = 0; n < NumAux; ++n) { + rhoaux1[n] = rhoaux0[n] + tau * ydot_a[SFX+n] + tau * rhoaux_source[n]; + } - eos_re_t eos_state; - eos_state.T = T_in; // initial guess + // compute the temperature at t0 + tau - // initial aux_sources - Real aux_source[NumAux] = {0.0_rt}; - aux_source[iye] = rho_old * nse_state.dyedt; + Real T1; + Real abar1_out; + Real Ye1 = rhoaux1[iye] / rho1; - Real rho_aux_new[NumAux]; + if (T_fixed > 0) { + T1 = T_fixed; + abar1_out = rhoaux1[iabar] / rho1; + } else { + Real e1 = rhoe1 / rho1; + T1 = T0; + nse_T_abar_from_e(rho1, e1, Ye1, T1, abar1_out); + } - Real rhoe_new; - Real rho_enucdot = -rho_old * snu; + // call NSE at t0 + tau - Real rho_half = 0.5_rt * (rho_old + state.y[SRHO]); + nse_state.T = T1; + nse_state.rho = rho1; + nse_state.Ye = Ye1; - // predict the new aux for the first iteration -- this is really - // just including the advection bits + nse_interp(nse_state, skip_X_fill); + + Real bea1_out = nse_state.bea; + + // construct the finite-difference approximation to the derivatives - for (int n = 0; n < NumAux; n++) { - rho_aux_new[n] = state.y[SFX+n] + dt * state.ydot_a[SFX+n] + dt * aux_source[n]; + // note that abar, (B/A) here have already seen advection + // implicitly + + Real rho_bea_tilde = rho1 * bea1_out - tau * ydot_a[SFX+ibea]; + Real rho_dBEA = rho_bea_tilde - rho0 * bea0_out; // this is MeV / nucleon * g / cm**3 + + Real rho_abar_tilde = rho1 * abar1_out - tau * ydot_a[SFX+iabar]; + Real rho_dabar = rho_abar_tilde - rho0 * abar0_out; // this is MeV / nucleon * g / cm**3 + + drhoedt = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / tau; + if (integrator_rp::nse_include_enu_weak == 1) { + drhoedt -= rho0 * (nse_state.e_nu + snu); + } else { + drhoedt -= rho0 * snu; } + drhoauxdt[iabar] = rho_dabar / tau; + drhoauxdt[iye] = rho0 * dyedt0; + drhoauxdt[ibea] = rho_dBEA / tau; - for (int iter = 0; iter < integrator_rp::nse_iters; iter++) { +} - // update (rho e)^{n+1} based on the new energy generation rate - rhoe_new = state.y[SEINT] + dt * state.ydot_a[SEINT] + dt * rho_enucdot; +/// +/// update the state due to NSE changes for the simplified-SDC case +/// this version works with the tabulated NSE and requires AUX_THERMO +/// +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void sdc_nse_burn(BurnT& state, const Real dt) { - // call the EOS to get the updated T* + using namespace AuxZero; - Real T_new; - eos_state.rho = state.y[SRHO]; - eos_state.e = rhoe_new / state.y[SRHO]; - for (int n = 0; n < NumAux; n++) { - eos_state.aux[n] = rho_aux_new[n] / state.y[SRHO]; - } + state.success = true; + state.n_rhs = 0; + state.n_jac = 0; - if (state.T_fixed > 0) { - T_new = state.T_fixed; - } else { - eos(eos_input_re, eos_state); - T_new = eos_state.T; - } + // store the initial state + + Real rho_old = state.y[SRHO]; + Real rhoe_old = state.y[SEINT]; + Real rhoaux_old[NumAux]; + for (int n = 0; n < NumAux; ++n) { + rhoaux_old[n] = state.y[SFX+n]; + } + + // density and momentum have no reactive sources - // call the NSE table using the * state to get the t^{n+1} - // source estimates + state.y[SRHO] += dt * state.ydot_a[SRHO]; + state.y[SMX] += dt * state.ydot_a[SMX]; + state.y[SMY] += dt * state.ydot_a[SMY]; + state.y[SMZ] += dt * state.ydot_a[SMZ]; - nse_state.T = T_new; - nse_state.rho = eos_state.rho; - nse_state.Ye = eos_state.aux[iye]; + // do an RK2 integration - nse_interp(nse_state); + // get the derivatives at t = t^n - // note that the abar, dq (B/A), and X here have all already - // seen advection implicitly + Real drhoedt; + Real drhoauxdt[NumSpec]; + nse_derivs(rho_old, rhoe_old, rhoaux_old, + dt, state.ydot_a, + drhoedt, drhoauxdt, state.T_fixed); - // compute the energy release -- we need to remove the advective part + // evolve to the midpoint in time - Real rho_bea_tilde = state.y[SRHO] * nse_state.bea - dt * state.ydot_a[SFX+ibea]; - Real rho_dBEA = rho_bea_tilde - rho_bea_old; // this is MeV / nucleon * g / cm**3 + Real rho_tmp = rho_old + 0.5_rt * dt * state.ydot_a[SRHO]; + Real rhoe_tmp = rhoe_old + 0.5_rt * dt * (state.ydot_a[SEINT] + drhoedt); + Real rhoaux_tmp[NumAux]; + for (int n = 0; n < NumAux; ++n) { + rhoaux_tmp[n] = rhoaux_old[n] + + 0.5_rt * dt * (state.ydot_a[SFX+n] + drhoauxdt[n]); + } - // convert the energy to erg / cm**3 - rho_enucdot = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; + // compute the derivatives at the midpoint in time - // now get the updated neutrino term - zbar = nse_state.abar * eos_state.aux[iye]; -#ifdef NEUTRINOS - sneut5(T_new, eos_state.rho, nse_state.abar, zbar, - snu, dsnudt, dsnudd, dsnuda, dsnudz); -#endif - rho_enucdot -= 0.5_rt * rho_half * (snu_old + snu); + nse_derivs(rho_tmp, rhoe_tmp, rhoaux_tmp, + dt, state.ydot_a, + drhoedt, drhoauxdt, state.T_fixed); + + // evolve to the new time - // update the new state for the next pass + Real rho_new = rho_old + dt * state.ydot_a[SRHO]; + Real rhoe_new = rhoe_old + dt * (state.ydot_a[SEINT] + drhoedt); + Real rhoaux_new[NumAux]; + for (int n = 0; n < NumAux; ++n) { + rhoaux_new[n] = rhoaux_old[n] + + dt * (state.ydot_a[SFX+n] + drhoauxdt[n]); + } - aux_source[iye] = 0.5_rt * rho_half * (dyedt_old + nse_state.dyedt); - rho_aux_new[iye] = state.y[SFX+iye] + dt * state.ydot_a[SFX+iye] + dt * aux_source[iye]; + // get the new temperature - rho_aux_new[iabar] = state.y[SRHO] * nse_state.abar; - rho_aux_new[ibea] = state.y[SRHO] * nse_state.bea; + Real T_new; + Real abar_new; + Real Ye_new = rhoaux_new[iye] / rho_new; + if (state.T_fixed > 0) { + T_new = state.T_fixed; + abar_new = rhoaux_new[iabar] / rho_new; + } else { + Real e_new = rhoe_new / rho_new; + T_new = 1.e8; // initial guess + nse_T_abar_from_e(rho_new, e_new, Ye_new, T_new, abar_new); } - // now update the aux quantities + // do a final NSE call -- we want the ending B/A to be consistent + // with our thermodynamic state here, plus we need to get the mass + // fractions + + + constexpr bool skip_X_fill{false}; + + nse_table_t nse_state; + nse_state.T = T_new; + nse_state.rho = rho_new; + nse_state.Ye = Ye_new; + + nse_interp(nse_state, skip_X_fill); + + // store the new state // the new mass fractions are just those that come from the table // make sure they are normalized + Real sum_X{0.0_rt}; for (auto & xn : nse_state.X) { xn = std::clamp(xn, small_x, 1.0_rt); @@ -195,21 +267,23 @@ void sdc_nse_burn(BurnT& state, const Real dt) { } for (int n = 0; n < NumSpec; ++n) { - state.y[SFS+n] = state.y[SRHO] * nse_state.X[n]; + state.y[SFS+n] = rho_new * nse_state.X[n]; } - // aux data comes from the table (except Ye, which we predicted) - - state.y[SFX+iye] = eos_state.rho * eos_state.aux[iye]; - state.y[SFX+iabar] = eos_state.rho * nse_state.abar; - state.y[SFX+ibea] = eos_state.rho * nse_state.bea; + // aux data comes from the integration or the table + state.y[SFX+iye] = rhoaux_new[iye]; + state.y[SFX+iabar] = rho_new * abar_new; + state.y[SFX+ibea] = rho_new * nse_state.bea; // density and momenta have already been updated - // update the total and internal energy now + // get the energy release from the change in rhoe over the timestep + // excluding any advection, and use that to update the total energy - state.y[SEINT] += dt * state.ydot_a[SEINT] + dt * rho_enucdot; + Real rho_enucdot = (rhoe_new - rhoe_old) / dt - state.ydot_a[SEINT]; + + state.y[SEINT] = rhoe_new; state.y[SEDEN] += dt * state.ydot_a[SEDEN] + dt * rho_enucdot; }