Skip to content

update the script that tests nse network compatibility #1811

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

Merged
zingale merged 1 commit into from
Jun 15, 2025
Merged

update the script that tests nse network compatibility #1811

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
6 changes: 3 additions & 3 deletions nse_solver/nse_compatibility/GNUmakefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,10 +22,12 @@ MICROPHYSICS_HOME ?= ../..
EOS_DIR := helmholtz

# This sets the network directory
NETWORK_DIR := subch_base
NETWORK_DIR := ase

SCREEN_METHOD = chabrier1998

USE_SIMPLIFIED_SDC = TRUE

CONDUCTIVITY_DIR := stellar

INTEGRATOR_DIR = VODE
Expand All @@ -40,5 +42,3 @@ Bpack := ./Make.package
Blocs := .

include $(MICROPHYSICS_HOME)/unit_test/Make.unit_test


13 changes: 4 additions & 9 deletions nse_solver/nse_compatibility/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,8 @@

This is a script for checking whether integrated mass fractions of
a reaction network will eventually match with the mass fractions
calculated via NSE equations. Currently, this is only valid for networks
that do not involve weak rates.
calculated via NSE equations.

* Note that the integration doesn't achieve NSE even though we have
USE_NSE_NET=TRUE is due to the reference size of the cell. By setting
nse.nse_dx_independet=1 will fix this.

Script will print out all burn_state cases when nse_dx_independent is
not enabled. It will only print out burn_state cases that didn't enter NSE
when nse_dx_independent=1
Script will print out the mass fractions computed via direct integration
and the mass fraction obtained from the NSE equation. It also calculates
the absolute and relative differences between the two.
20 changes: 11 additions & 9 deletions nse_solver/nse_compatibility/_parameters
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,22 +1,24 @@
@namespace: unit_test

run_prefix string ""
run_prefix string ""

# the final time to integrate to
tmax real 1.e3
tmax real 1.e-1

# first output time -- we will output in nsteps logarithmically spaced steps between [tfirst, tmax]
tfirst real 0.0

# number of steps (logarithmically spaced)
nsteps int 100
nsteps int 100

rho_min real 1.e7
rho_max real 1.e9
# density
rho real 1.e7

nrho int 4
# temperature
T real 6.e9

T_min real 6.e9
T_max real 8.e9
# Initial guess for proton chemical potential
mu_p real -3.0

nT int 4
# Initial guess for neutron chemical potential
mu_n real -12.0
12 changes: 6 additions & 6 deletions nse_solver/nse_compatibility/inputs_nse_compatibility
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,14 @@ integrator.jacobian = 1

integrator.renormalize_abundances = 0

integrator.rtol_spec = 1.0e-6
integrator.rtol_enuc = 1.0e-6
integrator.atol_spec = 1.0e-10
integrator.atol_enuc = 1.0e-6
integrator.rtol_spec = 1.0e-15
integrator.rtol_enuc = 1.0e-15
integrator.atol_spec = 1.0e-15
integrator.atol_enuc = 1.0e-15

unit_test.tmax = 1.e3
unit_test.tmax = 1.e0
unit_test.tfirst = 1.e-10
unit_test.nsteps = 100
unit_test.nsteps = 10000

integrator.integrate_energy = 0
integrator.call_eos_in_rhs = 0
Expand Down
143 changes: 58 additions & 85 deletions nse_solver/nse_compatibility/nse_network_compatibility.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,29 @@ void burn_cell_c(const eos_t& eos_state)

burn_state.rho = eos_state.rho;
burn_state.T = eos_state.T;
burn_state.e = eos_state.e;

burn_state.y[SRHO] = eos_state.rho;
burn_state.y[SMX] = 0.0;
burn_state.y[SMY] = 0.0;
burn_state.y[SMZ] = 0.0;

burn_state.y[SEINT] = eos_state.rho * eos_state.e;
burn_state.y[SEDEN] = burn_state.y[SEINT];

burn_state.ydot_a[SRHO] = 0.0;
burn_state.ydot_a[SEINT] = 0.0;

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

burn_state.y_e = eos_state.y_e;

burn_state.sdc_iter = 1;
burn_state.num_sdc_iters = 1;

burn_state.i = 0;
burn_state.j = 0;
burn_state.k = 0;
Expand Down Expand Up @@ -68,6 +86,13 @@ void burn_cell_c(const eos_t& eos_state)

amrex::Real energy_initial = burn_state.e;

// Now since we compiled with USE_NSE_NET=TRUE,
// it will try to burn using NSE, to prevent that, set extremely small tolerance

ase_tol = 1.e-100;
burn_state.mu_p = mu_p;
burn_state.mu_n = mu_n;

// loop over steps, burn, and output the current state

for (int n = 0; n < nsteps; n++){
Expand Down Expand Up @@ -96,34 +121,33 @@ void burn_cell_c(const eos_t& eos_state)
}

// Now find the mass fraction via NSE calculation
burn_state.mu_p = -3.0_rt;
burn_state.mu_n = -12.0_rt;
auto nse_state = get_actual_nse_state(burn_state);

if (nse_dx_independent != 1 || (nse_dx_independent == 1 && !burn_state.nse)) {
std::cout << std::scientific << std::setprecision(3);
std::cout << "------------------------------------" << std::endl;
std::cout << " - added e = " << burn_state.e - burn_state_in.e << std::endl;
std::cout << " Initial vs. Final T = " << burn_state_in.T << " "
<< burn_state.T << std::endl;
std::cout << " Initial vs. Final density = " << burn_state_in.rho << " "
<< burn_state.rho << std::endl;
std::cout << " Initial vs. Final vs. NSE y_e = " << burn_state_in.y_e << " "
<< burn_state.y_e << " " << nse_state.y_e << std::endl;
std::cout << "------------------------------------" << std::endl;
std::cout << "Element" << std::setw(14) << "Burn Xs" << std::setw(20) << "NSE Xs "
<< std::setw(20) << "abs err" << std::setw(20) << "rel err" << std::endl;
for (int n = 0; n < NumSpec; ++n) {
const std::string& element = short_spec_names_cxx[n];
amrex::Real abs_err = std::abs(burn_state.xn[n] - nse_state.xn[n]);
amrex::Real rel_err = abs_err / burn_state.xn[n];
std::cout << " " << element << std::setw(20-element.size())
<< burn_state.xn[n] << std::setw(20)
<< nse_state.xn[n] << std::setw(20) << abs_err << std::setw(20)
<< rel_err << std::endl;
}
std::cout << "------------------------------------" << std::endl;
auto nse_state = get_actual_nse_state(burn_state, 1.e-12);

// Now print out the result.
std::cout << std::scientific << std::setprecision(8);
std::cout << "------------------------------------" << std::endl;
std::cout << " - added e = " << burn_state.e - burn_state_in.e << std::endl;
std::cout << " Initial vs. Final T = " << burn_state_in.T << " "
<< burn_state.T << std::endl;
std::cout << " Initial vs. Final density = " << burn_state_in.rho << " "
<< burn_state.rho << std::endl;
std::cout << " Initial vs. Final vs. NSE y_e = " << burn_state_in.y_e << " "
<< burn_state.y_e << " " << nse_state.y_e << std::endl;
std::cout << " Solved Chemical Potentials are: " << burn_state.mu_p << " and "
<< burn_state.mu_n << std::endl;
std::cout << "------------------------------------" << std::endl;
std::cout << "Element" << std::setw(14) << "Burn Xs" << std::setw(20) << "NSE Xs "
<< std::setw(20) << "abs err" << std::setw(20) << "rel err" << std::endl;
for (int n = 0; n < NumSpec; ++n) {
const std::string& element = short_spec_names_cxx[n];
amrex::Real abs_err = std::abs(burn_state.xn[n] - nse_state.xn[n]);
amrex::Real rel_err = abs_err / burn_state.xn[n];
std::cout << " " << element << std::setw(20-element.size())
<< burn_state.xn[n] << std::setw(20)
<< nse_state.xn[n] << std::setw(20) << abs_err << std::setw(20)
<< rel_err << std::endl;
}
std::cout << "------------------------------------" << std::endl;
}

void nse_network_compatibility()
Expand All @@ -132,65 +156,14 @@ void nse_network_compatibility()
// with the results from NSE calculations given any network.
// Note that this network should not have any weak rates.

amrex::Real dlogrho = (std::log10(rho_max) - std::log10(rho_min))/static_cast<amrex::Real>(nrho-1);
amrex::Real dlogT = (std::log10(T_max) - std::log10(T_min))/static_cast<amrex::Real>(nT-1);

// Create massfraction 2d array to hold 3 different sets of initial mass fractions
// Corresponding to ye < 0.5, ye = 0.5, and ye > 0.5

amrex::Array2D<amrex::Real, 1, 3, 1, NumSpec+1, Order::C> massfractions;

bool condition1_met = false;
bool condition2_met = false;
bool condition3_met = false;

for (int n = 0; n < NumSpec; ++n) {
if (aion[n] == 2 * zion[n] && !condition1_met) {

// Set the first encountered nuclei with Z/A = 0.5
massfractions(1, n+1) = 0.7_rt;
massfractions(2, n+1) = 1.0_rt;
massfractions(3, n+1) = 0.7_rt;
condition1_met = true;
continue;
}
if (aion[n] > 2 * zion[n] && !condition2_met) {

// This is for achieving ye < 0.5
massfractions(1, n+1) = 0.3_rt;
condition2_met = true;
continue;
}
if (aion[n] < 2 * zion[n] && !condition3_met) {

// This is for achieving ye > 0.5
massfractions(3, n+1) = 0.3_rt;
condition3_met = true;
continue;
}
massfractions(1, n+1) = 0.0_rt;
massfractions(2, n+1) = 0.0_rt;
massfractions(3, n+1) = 0.0_rt;
}

eos_t eos_state;

for (int iye = 0; iye < 3; ++iye) {
for (int irho = 0; irho < nrho; ++irho) {
for (int itemp = 0; itemp < nT; ++itemp) {

amrex::Real T = std::pow(10.0, std::log10(T_min) + itemp * dlogT);
amrex::Real rho = std::pow(10.0, std::log10(rho_min) + irho * dlogrho);

eos_state.T = T;
eos_state.rho = rho;
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = massfractions(iye+1, n+1);
}
eos(eos_input_rt, eos_state);
burn_cell_c(eos_state);
}
}
eos_state.T = T;
eos_state.rho = rho;
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = 1.0_rt;
}
normalize_abundances_burn(eos_state);
eos(eos_input_rt, eos_state);
burn_cell_c(eos_state);
}
#endif
Loading