Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update rhoX for nse_solver for SDC and discard strang #1607

Merged
merged 15 commits into from
Jul 11, 2024
Merged

Update rhoX for nse_solver for SDC and discard strang #1607

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
efe69f1
update rhoX for nse solver when doing sdc
zhichen3 Jul 9, 2024
ffebad4
update y if use sdc
zhichen3 Jul 9, 2024
3256404
convert to rhoX and discard strang
zhichen3 Jul 9, 2024
b2011bf
update unit_test
zhichen3 Jul 9, 2024
6ad5d9c
Missing comma
zhichen3 Jul 9, 2024
f4d2c36
fix fcn
zhichen3 Jul 9, 2024
983af20
update benchmakr
zhichen3 Jul 9, 2024
11a946a
update benchmakr and initial guess for make_table
zhichen3 Jul 9, 2024
ab6b283
update benchmark
zhichen3 Jul 9, 2024
965642e
add _rt
zhichen3 Jul 9, 2024
431ad46
work with xn during nse calculations, but also store rhoX in the end
zhichen3 Jul 9, 2024
8840b17
change initial guess to avoid overflow
zhichen3 Jul 9, 2024
75e01c3
update benchmark
zhichen3 Jul 9, 2024
d2f4fca
fix comment
zhichen3 Jul 10, 2024
a8a103c
move var outside of loop
zhichen3 Jul 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 31 additions & 34 deletions nse_solver/nse_check.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,9 @@
// First check to see if we're in the ballpark of nse state

AMREX_GPU_HOST_DEVICE AMREX_INLINE
void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) {
void check_nse_molar(const amrex::Array1D<amrex::Real, 1, NumSpec>& Y,
const amrex::Array1D<amrex::Real, 1, NumSpec>& Y_nse,
bool& nse_check) {

// This function gives the first estimate whether we're in the nse or not
// it checks whether the molar fractions of n,p,a are approximately in NSE
Expand All @@ -44,13 +46,13 @@ void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) {
for (int n = 0; n < NumSpec; ++n) {

if (n == NSE_INDEX::H1_index || n == NSE_INDEX::N_index) {
r /= state.xn[n] * state.xn[n] * aion_inv[n] * aion_inv[n];
r_nse /= nse_state.xn[n] * nse_state.xn[n] * aion_inv[n] * aion_inv[n];
r /= Y(n+1) * Y(n+1);
r_nse /= Y_nse(n+1) * Y_nse(n+1);
}

else if (n == NSE_INDEX::He4_index) {
r *= state.xn[n] * aion_inv[n];
r_nse *= nse_state.xn[n] * aion_inv[n];
r *= Y(n+1);
r_nse *= Y_nse(n+1);
}
}

Expand All @@ -59,15 +61,15 @@ void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) {

// if there is neutron in the network

if ((std::abs(r - r_nse) < 0.5_rt*r_nse)
if ((std::abs(r - r_nse) < 0.5_rt * r_nse)
&& (NSE_INDEX::N_index != -1)) {
nse_check = true;
return;
}

// if there is no neutron in the network

if ((std::abs(r - r_nse) < 0.25_rt*r_nse)
if ((std::abs(r - r_nse) < 0.25_rt * r_nse)
&& (NSE_INDEX::N_index == -1)) {
nse_check = true;
return;
Expand All @@ -76,9 +78,8 @@ void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) {
// Overall molar fraction check

for (int n = 0; n < NumSpec; ++n) {
amrex::Real abs_diff = std::abs(state.xn[n] - nse_state.xn[n]) / aion[n];
amrex::Real rel_diff = abs_diff / (state.xn[n] / aion[n]);
if (abs_diff > nse_abs_tol && rel_diff > nse_rel_tol) {
amrex::Real abs_diff = std::abs(Y(n+1) - Y_nse(n+1));
if (abs_diff > nse_abs_tol && abs_diff > nse_rel_tol * Y(n+1)) {
return;
}
}
Expand Down Expand Up @@ -228,7 +229,7 @@ bool in_single_group(const amrex::Array1D<int, 1, NumSpec>& group_ind) {
template <typename T>
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void fill_reaction_timescale(amrex::Array1D<T, 1, Rates::NumRates>& reaction_timescales,
const int current_rate_index, const burn_t& state,
const int current_rate_index, const amrex::Real rho,
const amrex::Array1D<amrex::Real, 1, NumSpec>& Y,
const amrex::Array1D<amrex::Real, 1, Rates::NumRates>& screened_rates,
const amrex::Real t_s) {
Expand Down Expand Up @@ -301,14 +302,14 @@ void fill_reaction_timescale(amrex::Array1D<T, 1, Rates::NumRates>& reaction_tim
if (NSE_INDEX::rate_indices(current_rate_index, 2) == NSE_INDEX::rate_indices(current_rate_index, 3)) {
b_f *= 0.5_rt;
}
b_f *= Y(NSE_INDEX::rate_indices(current_rate_index, 2) + 1) * state.rho;
b_f *= Y(NSE_INDEX::rate_indices(current_rate_index, 2) + 1) * rho;
}

if (NSE_INDEX::rate_indices(current_rate_index, 5) != -1) {
if (NSE_INDEX::rate_indices(current_rate_index, 5) == NSE_INDEX::rate_indices(current_rate_index, 6)) {
b_r *= 0.5_rt;
}
b_r *= Y(NSE_INDEX::rate_indices(current_rate_index, 5) + 1) * state.rho;
b_r *= Y(NSE_INDEX::rate_indices(current_rate_index, 5) + 1) * rho;
}

// Find the timescale of the rate, See Equation 11 in Kushnir
Expand Down Expand Up @@ -437,7 +438,7 @@ void fill_merge_indices(amrex::Array1D<int, 1, 2>& merge_indices,


AMREX_GPU_HOST_DEVICE AMREX_INLINE
void nse_grouping(amrex::Array1D<int, 1, NumSpec>& group_ind, const burn_t& state,
void nse_grouping(amrex::Array1D<int, 1, NumSpec>& group_ind, const amrex::Real rho,
const amrex::Array1D<amrex::Real, 1, NumSpec>& Y,
const amrex::Array1D<amrex::Real, 1, Rates::NumRates>& screened_rates,
const amrex::Real t_s) {
Expand Down Expand Up @@ -468,7 +469,7 @@ void nse_grouping(amrex::Array1D<int, 1, NumSpec>& group_ind, const burn_t& stat
amrex::Array1D<int, 1, Rates::NumRates> rate_indices;

for (int n = 1; n <= Rates::NumRates; ++n) {
fill_reaction_timescale(reaction_timescales, n, state, Y,
fill_reaction_timescale(reaction_timescales, n, rho, Y,
screened_rates, t_s);
rate_indices(n) = n;
}
Expand Down Expand Up @@ -543,19 +544,22 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) {

const auto nse_state = get_actual_nse_state(current_state);

burn_t state = current_state;
auto state = current_state;

// Convert to molar fractions

amrex::Array1D<amrex::Real, 1, NumSpec> Y;
amrex::Array1D<amrex::Real, 1, NumSpec> Y_nse;

#ifndef STRANG
// if not strang, store mass fractions
for (int n = 0; n < NumSpec; ++n) {
state.xn[n] = current_state.y[SFS+n] / current_state.rho;
Y(n+1) = current_state.y[SFS+n] * aion_inv[n] / current_state.rho;
Y_nse(n+1) = nse_state.xn[n] * aion_inv[n];
}
#endif

// Check whether state is in the ballpark of NSE

if (!skip_molar_check) {
check_nse_molar(state, nse_state, current_state.nse);
check_nse_molar(Y, Y_nse, current_state.nse);
if (!current_state.nse) {
return current_state.nse;
}
Expand All @@ -577,23 +581,16 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) {
// our current mass fractions are in the ballpark of NSE mass fractions.

if (nse_molar_independent || skip_molar_check) {
state = nse_state;
state.dx = current_state.dx;
#ifndef STRANG
for (int n = 0; n < NumSpec; ++n) {
state.y[SFS+n] = nse_state.xn[n] * nse_state.rho;
for (int n = 1; n <= NumSpec; ++n) {
Y(n) = Y_nse(n);
}
#endif
}

// set molar fractions

amrex::Array1D<amrex::Real, 1, NumSpec> Y;
// store xn to find the rates

for (int n = 1; n <= NumSpec; ++n) {
Y(n) = state.xn[n-1] * aion_inv[n-1];
for (int n = 0; n < NumSpec; ++n) {
state.xn[n] = Y(n+1) * aion[n];
}

rate_t rate_eval;
constexpr int do_T_derivatives = 0;
evaluate_rates<do_T_derivatives, rate_t>(state, rate_eval);
Expand Down Expand Up @@ -629,7 +626,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) {

// Now do nse grouping

nse_grouping(group_ind, state, Y, rate_eval.screened_rates, t_s);
nse_grouping(group_ind, state.rho, Y, rate_eval.screened_rates, t_s);

// Check if we result in a single group after grouping

Expand Down
48 changes: 28 additions & 20 deletions nse_solver/nse_solver.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ T get_nonexponent_nse_state(const T& state) {
for (int n = 0; n < NumSpec; ++n) {
#ifdef NEW_NETWORK_IMPLEMENTATION
if (n == NSE_INDEX::H1_index) {
nse_state.xn[n] = 0.0;
nse_state.xn[n] = 0.0_rt;
continue;
}
#endif
Expand All @@ -73,11 +73,12 @@ T get_nonexponent_nse_state(const T& state) {

// find nse mass frac without the exponent term.

amrex::Real power = 2.0 * M_PI * network::mion(n+1) *
C::k_B * T_in / (C::hplanck * C::hplanck);
constexpr amrex::Real ihplanck2 = 1.0_rt / (C::hplanck * C::hplanck);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this can be moved outside of the loop

amrex::Real power = 2.0_rt * M_PI * network::mion(n+1) *
C::k_B * T_in * ihplanck2;

nse_state.xn[n] = network::mion(n+1) * pf * spin / state.rho *
std::sqrt(power * power * power);
std::sqrt(amrex::Math::powi<3>(power));
}

return nse_state;
Expand All @@ -103,9 +104,9 @@ void compute_coulomb_contribution(amrex::Array1D<amrex::Real, 1, NumSpec>& u_c,
// so we just treat u_c as a constant.
//

const amrex::Real n_e = state.rho * state.y_e / C::m_u;
const amrex::Real n_e = state.rho * state.y_e * C::Legacy::n_A;
const amrex::Real Gamma_e = C::q_e * C::q_e *
std::cbrt(4.0_rt * M_PI * n_e / 3.0_rt) / (C::k_B * T_in);
std::cbrt(1.333333333333_rt * M_PI * n_e) / (C::k_B * T_in);
amrex::Real gamma;

for (int n = 0; n < NumSpec; ++n) {
Expand All @@ -116,7 +117,7 @@ void compute_coulomb_contribution(amrex::Array1D<amrex::Real, 1, NumSpec>& u_c,
#endif
// term for calculating u_c

gamma = std::pow(zion[n], 5.0_rt/3.0_rt) * Gamma_e;
gamma = std::cbrt(amrex::Math::powi<5>(zion[n])) * Gamma_e;

// chemical potential for coulomb correction
// see appendix of Calder 2007, doi:10.1086/510709 for more detail
Expand Down Expand Up @@ -146,12 +147,13 @@ void apply_nse_exponent(T& nse_state,
// the temperature that comes in through T_fixed

amrex::Real T_in = nse_state.T_fixed > 0.0_rt ? nse_state.T_fixed : nse_state.T;
amrex::Real ikTMeV = C::Legacy::MeV2erg / (C::k_B * T_in);
amrex::Real exponent;

for (int n = 0; n < NumSpec; ++n) {
#ifdef NEW_NETWORK_IMPLEMENTATION
if (n == NSE_INDEX::H1_index) {
nse_state.xn[n] = 0.0;
nse_state.xn[n] = 0.0_rt;
continue;
}
#endif
Expand All @@ -163,8 +165,8 @@ void apply_nse_exponent(T& nse_state,

exponent = amrex::min(500.0_rt,
(zion[n] * nse_state.mu_p + (aion[n] - zion[n]) *
nse_state.mu_n + network::bion(n+1)) /
C::k_B / T_in * C::Legacy::MeV2erg - u_c(n+1));
nse_state.mu_n + network::bion(n+1)) *
ikTMeV - u_c(n+1));

nse_state.xn[n] *= std::exp(exponent);
}
Expand Down Expand Up @@ -210,6 +212,7 @@ void fcn(amrex::Array1D<amrex::Real, 1, 2>& x, amrex::Array1D<amrex::Real, 1, 2>
}
#endif
// constraint equation 1, mass fraction sum to 1
// since we're using rhoX, make it sum to rho

fvec(1) += nse_state.xn[n];
}
Expand Down Expand Up @@ -247,17 +250,18 @@ void jcn(amrex::Array1D<amrex::Real, 1, 2>& x, amrex::Array2D<amrex::Real, 1, 2,
// the temperature that comes in through T_fixed

amrex::Real T_in = nse_state.T_fixed > 0.0_rt ? nse_state.T_fixed : nse_state.T;
amrex::Real ikTMeV = C::Legacy::MeV2erg / (C::k_B * T_in);

for (int n = 0; n < NumSpec; ++n) {
#ifdef NEW_NETWORK_IMPLEMENTATION
if (n == NSE_INDEX::H1_index) {
continue;
}
#endif
fjac(1, 1) += nse_state.xn[n] * zion[n] / C::k_B / T_in * C::Legacy::MeV2erg ;
fjac(1, 2) += nse_state.xn[n] * (aion[n] - zion[n]) / C::k_B / T_in * C::Legacy::MeV2erg;
fjac(2, 1) += nse_state.xn[n] * zion[n] * zion[n] * aion_inv[n] / C::k_B / T_in * C::Legacy::MeV2erg;
fjac(2, 2) += nse_state.xn[n] * zion[n] * (aion[n] - zion[n]) * aion_inv[n] / C::k_B / T_in * C::Legacy::MeV2erg;
fjac(1, 1) += nse_state.xn[n] * zion[n] * ikTMeV;
fjac(1, 2) += nse_state.xn[n] * (aion[n] - zion[n]) * ikTMeV;
fjac(2, 1) += nse_state.xn[n] * zion[n] * zion[n] * aion_inv[n] * ikTMeV;
fjac(2, 2) += nse_state.xn[n] * zion[n] * (aion[n] - zion[n]) * aion_inv[n] * ikTMeV;
}

}
Expand Down Expand Up @@ -499,15 +503,13 @@ T get_actual_nse_state(T& state, amrex::Real eps=1.0e-10_rt,

if (!input_ye_is_valid) {
// ensure Ye is valid
#ifdef STRANG
composition(state);
#else
state.y_e = 0.0_rt;
// We assume the input state has rhoX

state.y_e = 0.0_rt;
for (int n = 0; n < NumSpec; ++n) {
state.y_e += state.y[SFS+n] * zion[n] * aion_inv[n] / state.rho;
state.y_e += state.y[SFS+n] * zion[n] * aion_inv[n];
}
#endif
state.y_e /= state.rho;
}

nse_solver_data<T> state_data = {state, {0.0_rt}};
Expand Down Expand Up @@ -557,6 +559,12 @@ T get_actual_nse_state(T& state, amrex::Real eps=1.0e-10_rt,
state.mu_p = state_data.state.mu_p;
state.mu_n = state_data.state.mu_n;

// update rhoX in the output nse_state

for (int n = 0; n < NumSpec; ++n) {
state_data.state.y[SFS+n] = state.rho * state_data.state.xn[n];
}

return state_data.state;
}
#endif
5 changes: 2 additions & 3 deletions unit_test/test_ase/GNUmakefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,11 @@ NETWORK_DIR := ase
CONDUCTIVITY_DIR := stellar

SCREEN_METHOD := chabrier1998

#SCREEN_METHOD := null
#SCREEN_METHOD := chugunov2007

USE_SIMPLIFIED_SDC = TRUE

INTEGRATOR_DIR = VODE

EXTERN_SEARCH += .
Expand All @@ -40,5 +41,3 @@ Bpack := ./Make.package
Blocs := .

include $(MICROPHYSICS_HOME)/unit_test/Make.unit_test


6 changes: 3 additions & 3 deletions unit_test/test_ase/burn_cell.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,12 +50,12 @@ void burn_cell_c()
std::cout << short_spec_names_cxx[n] << " : " << NSE_STATE.xn[n] << std::endl;
}

// Let state.xn equal to nse_state.xn to make sure its in nse state
// Let state.y equal to nse_state.y to make sure its in nse state

for (int n = 0; n < NumSpec; ++n){
state.xn[n] = NSE_STATE.xn[n];
state.y[SFS+n] = NSE_STATE.y[SFS+n];
}


// get eos
//eos(eos_input_rt, state);

Expand Down
46 changes: 23 additions & 23 deletions unit_test/test_ase/ci-benchmarks/ase_nse_net_unit_test.out
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,28 +8,28 @@ State Density (g/cm^3): 10000000
State Temperature (K): 6000000000
electron fraction is 0.5
NSE state:
n : 0.0006075693495
H1 : 0.001992432534
He4 : 0.519007877
C12 : 1.326102e-05
N13 : 9.477376691e-11
N14 : 2.565486854e-09
O16 : 3.225455246e-05
F18 : 5.422064179e-11
Ne20 : 7.427509432e-07
Ne21 : 4.998912611e-08
Na22 : 2.608008186e-09
Na23 : 1.00884295e-06
Mg24 : 0.0001022222509
Al27 : 0.0005548826609
Si28 : 0.03680966846
P31 : 0.04229223791
S32 : 0.03865771985
Ar36 : 0.02464742443
Ca40 : 0.02885152622
Ti44 : 0.00166123105
Cr48 : 0.009967017598
Fe52 : 0.05968005236
Ni56 : 0.2351208159
n : 0.000607569351
H1 : 0.001992432526
He4 : 0.5190078752
C12 : 1.326101988e-05
N13 : 9.477376579e-11
N14 : 2.56548683e-09
O16 : 3.225455214e-05
F18 : 5.422064122e-11
Ne20 : 7.427509352e-07
Ne21 : 4.99891257e-08
Na22 : 2.608008157e-09
Na23 : 1.008842941e-06
Mg24 : 0.0001022222498
Al27 : 0.0005548826564
Si28 : 0.03680966808
P31 : 0.0422922376
S32 : 0.0386577195
Ar36 : 0.02464742425
Ca40 : 0.0288515261
Ti44 : 0.001661231048
Cr48 : 0.009967017622
Fe52 : 0.05968005276
Ni56 : 0.2351208186
We're in NSE.
AMReX (23.07-395-ge6c93bf22695) finalized
Loading
Loading