Skip to content

Commit

Permalink
start of a new NSE table interface
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale committed Dec 7, 2023
1 parent 49c9712 commit 07fa505
Show file tree
Hide file tree
Showing 3 changed files with 129 additions and 48 deletions.
50 changes: 25 additions & 25 deletions integration/nse_update_simplified_sdc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#include <eos.H>

#ifdef NSE_TABLE
#include <nse_table_type.H>
#include <nse_table.H>
#endif
#ifdef NSE_NET
Expand Down Expand Up @@ -85,18 +86,14 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
// get the current NSE state from the table -- this will be used
// to compute the NSE evolution sources

Real abar_out;
Real dq_out;
Real dyedt;
Real dabardt;
Real dbeadt;
Real enu;
Real X[NumSpec];
nse_table_t nse_state;
nse_state.T = T_in;
nse_state.rho = state.rho;
nse_state.Ye = ye_inl

nse_interp(T_in, state.rho, ye_in,
abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X);
nse_interp(nse_state);

Real dyedt_old = dyedt;
Real dyedt_old = nse_state.dyedt;


// predict the U^{n+1,*} state with only estimates of the aux
Expand All @@ -107,7 +104,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) {

// initial aux_sources
Real aux_source[NumAux] = {0.0_rt};
aux_source[iye] = rho_old * dyedt;
aux_source[iye] = rho_old * nse_state.dyedt;

Real rho_aux_new[NumAux];

Expand Down Expand Up @@ -147,36 +144,39 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
// call the NSE table using the * state to get the t^{n+1}
// source estimates

nse_interp(T_new, eos_state.rho, eos_state.aux[iye],
abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X);
nse_state.T = T_new;
nse_state.rho = eos_state.rho;
nse_state.Ye = eos_state.aux[iye];

nse_interp(nse_state);

// note that the abar, dq (B/A), and X here have all already
// seen advection implicitly


// compute the energy release -- we need to remove the advective part

Real rho_bea_tilde = state.y[SRHO] * dq_out - dt * state.ydot_a[SFX+ibea];
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

// convert the energy to erg / cm**3
rho_enucdot = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt;

// now get the updated neutrino term
zbar = abar_out * eos_state.aux[iye];
zbar = nse_state.abar * eos_state.aux[iye];
#ifdef NEUTRINOS
sneut5<do_derivatives>(T_new, eos_state.rho, abar_out, zbar,
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);

// update the new state for the next pass

aux_source[iye] = 0.5_rt * rho_half * (dyedt_old + dyedt);
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];

rho_aux_new[iabar] = state.y[SRHO] * abar_out;
rho_aux_new[ibea] = state.y[SRHO] * dq_out;
rho_aux_new[iabar] = state.y[SRHO] * nse_state.abar;
rho_aux_new[ibea] = state.y[SRHO] * nse_state.bea;

}

Expand All @@ -185,24 +185,24 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
// 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 : X) {
xn = amrex::max(small_x, amrex::min(1.0_rt, xn));
for (auto & xn : nse_state.X) {
xn = std::clamp(xn, small_x, 1.0_rt);
sum_X += xn;
}

for (auto & xn : X) {
for (auto & xn : nse_state.X) {
xn /= sum_X;
}

for (int n = 0; n < NumSpec; ++n) {
state.y[SFS+n] = state.y[SRHO] * X[n];
state.y[SFS+n] = state.y[SRHO] * 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 * abar_out;
state.y[SFX+ibea] = eos_state.rho * dq_out;
state.y[SFX+iabar] = eos_state.rho * nse_state.abar;
state.y[SFX+ibea] = eos_state.rho * nse_state.bea;


// density and momenta have already been updated
Expand Down
1 change: 1 addition & 0 deletions nse_tabular/Make.package
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,4 @@ CEXE_headers += nse_table.H
CEXE_headers += nse_table_check.H
CEXE_headers += nse_table_data.H
CEXE_sources += nse_table_data.cpp
CEXE_headers += nse_table_type.H
126 changes: 103 additions & 23 deletions nse_tabular/nse_table.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,13 @@
#include <AMReX_Array.H>
#include <AMReX_REAL.H>

#include <eos.H>

#include <extern_parameters.H>
#include <nse_table_data.H>
#include <nse_table_size.H>
#include <nse_table_type.H>


using namespace amrex;
using namespace network_rp;
Expand Down Expand Up @@ -235,37 +239,35 @@ Real tricubic(const int ir0, const int it0, const int ic0,
}

AMREX_GPU_HOST_DEVICE AMREX_INLINE
void nse_interp(const Real T, const Real rho, const Real ye,
Real& abar, Real& bea, Real& dyedt, Real& dabardt, Real& dbeadt, Real& e_nu,
Real* X, bool skip_X_fill=false) {
void nse_interp(nse_table_t& nse_state, bool skip_X_fill=false) {

// if skip_X_fill = true then we don't fill X[] with the mass fractions.

using namespace nse_table;
using namespace AuxZero;

Real rholog = std::log10(rho);
Real rholog = std::log10(nse_state.rho);
{
Real rmin = nse_table_size::logrho_min;
Real rmax = nse_table_size::logrho_max;

rholog = std::max(rmin, std::min(rmax, rholog));
rholog = std::clamp(rholog, rmin, rmax);
}

Real tlog = std::log10(T);
Real tlog = std::log10(nse_state.T);
{
Real tmin = nse_table_size::logT_min;
Real tmax = nse_table_size::logT_max;

tlog = std::max(tmin, std::min(tmax, tlog));
tlog = std::clamp(tlog, tmin, tmax);
}

Real yet = ye;
Real yet = nse_state.Ye;
{
Real yemin = nse_table_size::ye_min;
Real yemax = nse_table_size::ye_max;

yet = std::max(yemin, std::min(yemax, yet));
yet = std::clamp(yet, yemin, yemax);
}

if (nse_table_interp_linear) {
Expand All @@ -274,12 +276,12 @@ void nse_interp(const Real T, const Real rho, const Real ye,
int it1 = nse_get_logT_index(tlog);
int ic1 = nse_get_ye_index(yet);

abar = trilinear(ir1, it1, ic1, rholog, tlog, yet, abartab);
bea = trilinear(ir1, it1, ic1, rholog, tlog, yet, beatab);
dyedt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dyedttab);
dabardt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dabardttab);
dbeadt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dbeadttab);
e_nu = trilinear(ir1, it1, ic1, rholog, tlog, yet, enutab);
nse_state.abar = trilinear(ir1, it1, ic1, rholog, tlog, yet, abartab);
nse_state.bea = trilinear(ir1, it1, ic1, rholog, tlog, yet, beatab);
nse_state.dyedt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dyedttab);
nse_state.dabardt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dabardttab);
nse_state.dbeadt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dbeadttab);
nse_state.e_nu = trilinear(ir1, it1, ic1, rholog, tlog, yet, enutab);

// massfractab is 2-d, so we wrap the access in a lambda already
// indexing the component
Expand All @@ -288,7 +290,7 @@ void nse_interp(const Real T, const Real rho, const Real ye,
for (int n = 1; n <= NumSpec; n++) {
Real _X = trilinear(ir1, it1, ic1, rholog, tlog, yet,
[=] (const int i) {return massfractab(n, i);});
X[n-1] = amrex::max(0.0_rt, amrex::min(1.0_rt, _X));
nse_state.X[n-1] = std::clamp(_X, 0.0_rt, 1.0_rt);
}
}

Expand All @@ -308,12 +310,12 @@ void nse_interp(const Real T, const Real rho, const Real ye,
int ic0 = nse_get_ye_index(yet) - 1;
ic0 = amrex::max(1, amrex::min(nse_table_size::nye-3, ic0));

abar = tricubic(ir0, it0, ic0, rholog, tlog, yet, abartab);
bea = tricubic(ir0, it0, ic0, rholog, tlog, yet, beatab);
dyedt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dyedttab);
dabardt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dabardttab);
dbeadt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dbeadttab);
e_nu = tricubic(ir0, it0, ic0, rholog, tlog, yet, enutab);
nse_state.abar = tricubic(ir0, it0, ic0, rholog, tlog, yet, abartab);
nse_state.bea = tricubic(ir0, it0, ic0, rholog, tlog, yet, beatab);
nse_state.dyedt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dyedttab);
nse_state.dabardt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dabardttab);
nse_state.dbeadt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dbeadttab);
nse_state.e_nu = tricubic(ir0, it0, ic0, rholog, tlog, yet, enutab);

// massfractab is 2-d, so we wrap the access in a lambda already
// indexing the component
Expand All @@ -322,11 +324,89 @@ void nse_interp(const Real T, const Real rho, const Real ye,
for (int n = 1; n <= NumSpec; n++) {
Real _X = tricubic(ir0, it0, ic0, rholog, tlog, yet,
[=] (const int i) {return massfractab(n, i);});
X[n-1] = amrex::max(0.0_rt, amrex::min(1.0_rt, _X));
nse_state.X[n-1] = std::clamp(_X, 0.0_rt, 1.0_rt);
}
}
}

}


///
/// if we are in NSE, then the entire thermodynamic state is just
/// a function of rho, T, Ye. We can write the energy as:
///
/// e = e(rho, T, Y_e, Abar(rho, T, Ye))
///
/// where we note that Abar is a function of those same inputs.
/// This function inverts this form of the EOS to find the T
/// that satisfies the EOS and NSE.
///
AMREX_GPU_HOST_DEVICE AMREX_INLINE
Real
nse_T_from_e(const Real rho, const Real e_in, const Real Ye, const Real T_guess) {

using namespace AuxZero;

const Real ttol{1.e-8_rt};
const int max_iter{100};
const Real eps{1.e-8_rt};

// we need the full EOS type, since we need de/dA
//eos_extra_t eos_state;

bool converged{false};

Real T{T_guess};
int iter{};

while (not converged && iter < max_iter) {

// call NSE table to get Abar
nse_table_t nse_state;
nse_state.T = T;
nse_state.rho = rho;
nse_state.Ye = Ye;

constexpr bool skip_X_fill{true};
nse_interp(nse_state, skip_X_fill);
Real abar_old = nse_state.abar;

// call the EOS with the initial guess for T

eos_state.rho = rho;
eos_state.T = T;
eos_state.aux[iye] = Ye;
eos_state.aux[iabar] = abar_old;
eos(eos_input_rt, eos_state);

// f is the quantity we want to zero

Real f = eos_state.e - e_in;

// perturb T and call NSE again to get perturbed Abar
nse_state.T *= (1.0_rt + eps);
nse_interp(nse_state,skip_X_fill);

// estimate dAbar/dT
Real dabar_dT = (nse_state.abar - abar_old) / (eps * T);

// compute the correction to our guess
Real dT = f / (eos_state.dedT + eos_state.dedA * dabar_dT);

// update the temperature
T = std::clamp(T + dT, 0.5 * T, 2.0 * T);

// check convergence

if (std::abs(dT) < ttol * T) {
converged = true;
}
iter++;
}

return T;

}

#endif

0 comments on commit 07fa505

Please sign in to comment.