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 ECSN to latest pynucastro #1368

Merged
merged 2 commits into from
Oct 23, 2023
Merged

update ECSN to latest pynucastro #1368

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
89a1b8d
update ECSN to latest pynucastro
zingale Oct 23, 2023
e9d4be4
update
zingale Oct 23, 2023
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
19 changes: 10 additions & 9 deletions networks/ECSN/actual_rhs.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -268,37 +268,39 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {

[[maybe_unused]] Real rate, drate_dt, edot_nu, edot_gamma;

rate_eval.enuc_weak = 0.0;

tabular_evaluate(j_f20_o20_meta, j_f20_o20_rhoy, j_f20_o20_temp, j_f20_o20_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_f20_to_o20) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_f20_to_o20) = drate_dt;
}
rate_eval.add_energy_rate(k_f20_to_o20) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);

tabular_evaluate(j_ne20_f20_meta, j_ne20_f20_rhoy, j_ne20_f20_temp, j_ne20_f20_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_ne20_to_f20) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_ne20_to_f20) = drate_dt;
}
rate_eval.add_energy_rate(k_ne20_to_f20) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne20) * (edot_nu + edot_gamma);

tabular_evaluate(j_o20_f20_meta, j_o20_f20_rhoy, j_o20_f20_temp, j_o20_f20_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_o20_to_f20) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_o20_to_f20) = drate_dt;
}
rate_eval.add_energy_rate(k_o20_to_f20) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(O20) * (edot_nu + edot_gamma);

tabular_evaluate(j_f20_ne20_meta, j_f20_ne20_rhoy, j_f20_ne20_temp, j_f20_ne20_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_f20_to_ne20) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_f20_to_ne20) = drate_dt;
}
rate_eval.add_energy_rate(k_f20_to_ne20) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);


}
Expand Down Expand Up @@ -394,6 +396,7 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
rate_t rate_eval;

constexpr int do_T_derivatives = 0;

evaluate_rates<do_T_derivatives, rate_t>(state, rate_eval);

rhs_nuc(state, ydot, Y, rate_eval.screened_rates);
Expand All @@ -403,11 +406,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
Real enuc;
ener_gener_rate(ydot, enuc);

// include reaction neutrino losses (non-thermal) and gamma heating
enuc += C::Legacy::n_A * Y(F20) * rate_eval.add_energy_rate(k_f20_to_o20);
enuc += C::Legacy::n_A * Y(Ne20) * rate_eval.add_energy_rate(k_ne20_to_f20);
enuc += C::Legacy::n_A * Y(O20) * rate_eval.add_energy_rate(k_o20_to_f20);
enuc += C::Legacy::n_A * Y(F20) * rate_eval.add_energy_rate(k_f20_to_ne20);
// include any weak rate neutrino losses
enuc += rate_eval.enuc_weak;

// Get the thermal neutrino losses

Expand Down Expand Up @@ -613,6 +613,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
rate_derivs_t rate_eval;

constexpr int do_T_derivatives = 1;

evaluate_rates<do_T_derivatives, rate_derivs_t>(state, rate_eval);

// Species Jacobian elements with respect to other species
Expand Down
4 changes: 2 additions & 2 deletions networks/ECSN/reaclib_rates.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,13 +13,13 @@ using namespace Species;

struct rate_t {
Array1D<Real, 1, NumRates> screened_rates;
Array1D<Real, NrateReaclib+1, NrateReaclib+NrateTabular> add_energy_rate;
Real enuc_weak;
};

struct rate_derivs_t {
Array1D<Real, 1, NumRates> screened_rates;
Array1D<Real, 1, NumRates> dscreened_rates_dT;
Array1D<Real, NrateReaclib+1, NrateReaclib+NrateTabular> add_energy_rate;
Real enuc_weak;
};


Expand Down
13 changes: 11 additions & 2 deletions networks/ECSN/table_rates.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,11 @@ void init_tab_info(const table_t& tf, const std::string& file, R& log_rhoy_table
std::ifstream table;
table.open(file);

if (!table.is_open()) {
// the table was not present or we could not open it; abort
amrex::Error("table could not be opened");
}

std::string line;

// read and skip over the header
Expand All @@ -99,6 +104,10 @@ void init_tab_info(const table_t& tf, const std::string& file, R& log_rhoy_table
for (int j = 1; j <= tf.nrhoy; ++j) {
for (int i = 1; i <= tf.ntemp; ++i) {
std::getline(table, line);
if (line.empty()) {
amrex::Error("Error reading table data");
}

std::istringstream sdata(line);

sdata >> log_rhoy_table(j) >> log_temp_table(i);
Expand Down Expand Up @@ -265,8 +274,8 @@ evaluate_dr_dtemp(const table_t& table_meta, const R& log_rhoy_table, const T& l
// Finally, we perform a 1d-linear interpolation between dlogr_dlogt_ip1
// and dlogr_dlogt_i to compute dlogr_dlogt

Real log_rhoy_lo = log_temp_table(irhoy_lo);
Real log_rhoy_hi = log_temp_table(irhoy_hi);
Real log_rhoy_lo = log_rhoy_table(irhoy_lo);
Real log_rhoy_hi = log_rhoy_table(irhoy_hi);

Real log_temp_lo = log_temp_table(jtemp_lo);
Real log_temp_hi = log_temp_table(jtemp_hi);
Expand Down
50 changes: 25 additions & 25 deletions unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
AMReX (23.07-106-gde7c6189623b-dirty) initialized
AMReX (23.10-13-gecaa14456e3f) initialized
starting the single zone burn...
reading in network electron-capture / beta-decay tables...
Maximum Time (s): 0.01
Expand Down Expand Up @@ -29,37 +29,37 @@ RHS at t = 0
s32 1323519.174
------------------------------------
successful? 1
- Hnuc = 4.716248406e+19
- added e = 4.716248406e+17
- final T = 6804058463
- Hnuc = 4.4691188e+19
- added e = 4.4691188e+17
- final T = 6711517341
------------------------------------
e initial = 5.995956082e+17
e final = 1.071220449e+18
e final = 1.046507488e+18
------------------------------------
new mass fractions:
h1 0.009032124504
he4 4.003293528e-13
o16 5.888714278e-07
o20 0.00985678636
f20 0.009855432543
ne20 7.501479981e-14
mg24 1.21710929e-09
al27 1.951815457e-21
si28 0.26794899
p31 8.245884489e-14
s32 0.7033060765
h1 0.007036728874
he4 4.569221711e-13
o16 6.11671909e-07
o20 0.009363850335
f20 0.009362566168
ne20 8.595667728e-14
mg24 8.576669062e-06
al27 9.999999996e-31
si28 0.2286832772
p31 9.999999996e-31
s32 0.745544389
------------------------------------
species creation rates:
omegadot(h1): -0.0967875496
omegadot(h1): -0.2963271126
omegadot(he4): -3
omegadot(o16): -49.99994111
omegadot(o20): -0.01432136399
omegadot(f20): -0.0144567457
omegadot(o16): -49.99993883
omegadot(o20): -0.0636149665
omegadot(f20): -0.06374338321
omegadot(ne20): -30
omegadot(mg24): -9.999999878
omegadot(mg24): -9.999142333
omegadot(al27): -1
omegadot(si28): 25.794899
omegadot(si28): 21.86832772
omegadot(p31): -1
omegadot(s32): 69.33060765
number of steps taken: 31502
AMReX (23.07-106-gde7c6189623b-dirty) finalized
omegadot(s32): 73.5544389
number of steps taken: 14714
AMReX (23.10-13-gecaa14456e3f) finalized
2 changes: 1 addition & 1 deletion unit_test/test_rhs/ci-benchmarks/ecsn.out
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,7 +161,7 @@
J_hydrogen-1_E -0.12580793235 1.2089803642e+12
J_helium-4_E -0.0078065301852 1.4905283338e+12
J_oxygen-16_E -5.3990173959e+12 1.2739583515e-06
J_oxygen-20_E -2.079874023e-18 3.5544574816e-18
J_oxygen-20_E -2.8880340474e-19 3.5544574816e-18
J_fluorine-20_E -2.0716943987e-18 7.6531304664e-18
J_neon-20_E -8.8836179369e-05 1.9862184663e-05
J_magnesium-24_E -0.0040169125096 0.018563573855
Expand Down