Skip to content

Commit

Permalink
start of a new NSE table interface (#1404)
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale authored Dec 8, 2023
1 parent 49c9712 commit bdf6a18
Show file tree
Hide file tree
Showing 6 changed files with 100 additions and 81 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_in;

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
28 changes: 13 additions & 15 deletions integration/nse_update_strang.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 All @@ -37,13 +38,6 @@ void nse_burn(BurnT& state, const Real dt) {
using namespace AuxZero;

// use the NSE table
Real abar_out;
Real dq_out;
Real dyedt;
Real dabardt;
Real dbeadt;
Real enu;
Real X[NumSpec];

// get the energy -- we are assuming that rho, T are valid on input

Expand All @@ -58,25 +52,29 @@ void nse_burn(BurnT& state, const Real dt) {
// source estimates. The thermodynamnics here is specified
// in terms of the auxiliary composition, Ye, abar, and B/A

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

nse_interp(nse_state);

// update Ye

state.aux[iye] += dt * dyedt;
state.aux[iye] += dt * nse_state.dyedt;

// now get the composition from the table using the updated Ye

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


// this is MeV / nucleon -- here aux has not yet been updated, so we
// access the old binding energy
Real deltaq = dq_out - state.aux[ibea];
Real deltaq = nse_state.bea - state.aux[ibea];

state.aux[ibea] += deltaq;
state.aux[iabar] = abar_out;
state.aux[iabar] = nse_state.abar;

#elif defined(NSE_NET)

Expand Down Expand Up @@ -116,7 +114,7 @@ void nse_burn(BurnT& state, const Real dt) {

#if defined(NSE_TABLE)
for (int n = 0; n < NumSpec; ++n) {
state.xn[n] = X[n];
state.xn[n] = nse_state.X[n];
}
#elif defined(NSE_NET)
for (int n = 0; n < NumSpec; ++n) {
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
46 changes: 23 additions & 23 deletions nse_tabular/nse_table.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@
#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 +237,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 +274,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 +288,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 +308,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,7 +322,7 @@ 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);
}
}
}
Expand Down
23 changes: 23 additions & 0 deletions nse_tabular/nse_table_type.H
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
#ifndef NSE_TABLE_TYPE_H
#define NSE_TABLE_TYPE_H

#include <AMReX_REAL.H>
#include <network.H>

struct nse_table_t {

amrex::Real T;
amrex::Real rho;
amrex::Real Ye;
amrex::Real abar;
amrex::Real bea;
amrex::Real dyedt;
amrex::Real dabardt;
amrex::Real dbeadt;
amrex::Real e_nu;
amrex::Real X[NumSpec];

};


#endif
33 changes: 15 additions & 18 deletions unit_test/test_nse_interp/nse_cell.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include <iostream>

#include <nse_table.H>
#include <nse_table_type.H>

using namespace unit_test_rp;

Expand Down Expand Up @@ -120,25 +121,21 @@ void nse_cell_c()

std::cout << "tricubic interpolated values: " << std::endl;

Real abar;
Real bea;
Real dyedt;
Real dabardt;
Real dbeadt;
Real e_nu;
Real X[NumSpec];

nse_interp(temperature, density, ye,
abar, bea, dyedt, dabardt, dbeadt, e_nu, X);

std::cout << "abar = " << abar << std::endl;
std::cout << "bea = " << bea << std::endl;
std::cout << "dyedt = " << dyedt << std::endl;
std::cout << "dabardt = " << dabardt << std::endl;
std::cout << "dbeadt = " << dbeadt << std::endl;
std::cout << "e_nu = " << e_nu << std::endl;
nse_table_t nse_state;
nse_state.T = temperature;
nse_state.rho = density;
nse_state.Ye = ye;

nse_interp(nse_state);

std::cout << "abar = " << nse_state.abar << std::endl;
std::cout << "bea = " << nse_state.bea << std::endl;
std::cout << "dyedt = " << nse_state.dyedt << std::endl;
std::cout << "dabardt = " << nse_state.dabardt << std::endl;
std::cout << "dbeadt = " << nse_state.dbeadt << std::endl;
std::cout << "e_nu = " << nse_state.e_nu << std::endl;
for (int n = 0; n < NumSpec; ++n) {
std::cout << "X(" << short_spec_names_cxx[n] << ") = " << X[n] << std::endl;
std::cout << "X(" << short_spec_names_cxx[n] << ") = " << nse_state.X[n] << std::endl;
}


Expand Down

0 comments on commit bdf6a18

Please sign in to comment.