Skip to content

Commit

Permalink
a version of SDC+NSE RK2 with the NSE EOS (#1415)
Browse files Browse the repository at this point in the history
This implements the SDC+NSE table update as a 2nd order Runge-Kutta.
It uses the new NSE EOS interface to get an Abar / T consistent with the updated e
  • Loading branch information
zingale authored Dec 29, 2023
1 parent a38a1ee commit 2ea84b9
Showing 2 changed files with 195 additions and 113 deletions.
8 changes: 8 additions & 0 deletions integration/_parameters
Original file line number Diff line number Diff line change
@@ -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
300 changes: 187 additions & 113 deletions integration/nse_update_simplified_sdc.H
Original file line number Diff line number Diff line change
@@ -14,10 +14,9 @@
#include <extern_parameters.H>

#include <burn_type.H>
#include <eos.H>
#include <nse_eos.H>

#ifdef NSE_TABLE
#include <nse_table_type.H>
#include <nse_table.H>
#endif
#ifdef NSE_NET
@@ -27,163 +26,236 @@
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 <typename BurnT>
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 <B/A> 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};
Real dsnudd{0.0};
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<do_derivatives>(T_in, rho_old, abar, zbar,
sneut5<do_derivatives>(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 <typename BurnT>
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<do_derivatives>(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;

}

0 comments on commit 2ea84b9

Please sign in to comment.