Skip to content

Commit

Permalink
update self-consistent NSE SDC update with Tabular NSE (#1569)
Browse files Browse the repository at this point in the history 
the self-consistent and tabular code paths are now consistent -- both do an update that should be second-order accurate.
  • Loading branch information
@zhichen3
zhichen3 authored Jul 14, 2024
1 parent b413290 commit a5408ba
Show file tree
Hide file tree
Showing 16 changed files with 1,734 additions and 143 deletions.
363 changes: 234 additions & 129 deletions integration/nse_update_sdc.H
View file

Large diffs are not rendered by default.

80 changes: 78 additions & 2 deletions networks/CNO_extras/actual_rhs.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -693,11 +693,87 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {

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

rate_eval.enuc_weak = 0.0;
rate_eval.enuc_weak = 0.0_rt;


}

#ifdef NSE_NET
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void get_ydot_weak(const burn_t& state,
amrex::Array1D<amrex::Real, 1, neqs>& ydot_nuc,
amrex::Real& enuc_weak,
[[maybe_unused]] const amrex::Array1D<amrex::Real, 1, NumSpec>& Y) {
///
/// Calculate Ydots contribute only from weak reactions.
/// This is used to calculate dyedt and energy generation from
/// weak reactions for self-consistent NSE
///


// initialize ydot_nuc to 0

for (int i = 1; i <= neqs; ++i) {
ydot_nuc(i) = 0.0_rt;
}

rate_t rate_eval;

[[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma;
[[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e;

rate_eval.enuc_weak = 0.0_rt;

// Calculate tabular rates and get ydot_weak


ydot_nuc(H1) = 0.0_rt;

ydot_nuc(He4) = 0.0_rt;

ydot_nuc(C12) = 0.0_rt;

ydot_nuc(C13) = 0.0_rt;

ydot_nuc(N13) = 0.0_rt;

ydot_nuc(N14) = 0.0_rt;

ydot_nuc(N15) = 0.0_rt;

ydot_nuc(O14) = 0.0_rt;

ydot_nuc(O15) = 0.0_rt;

ydot_nuc(O16) = 0.0_rt;

ydot_nuc(O17) = 0.0_rt;

ydot_nuc(O18) = 0.0_rt;

ydot_nuc(F17) = 0.0_rt;

ydot_nuc(F18) = 0.0_rt;

ydot_nuc(F19) = 0.0_rt;

ydot_nuc(Ne18) = 0.0_rt;

ydot_nuc(Ne19) = 0.0_rt;

ydot_nuc(Ne20) = 0.0_rt;

ydot_nuc(Mg22) = 0.0_rt;

ydot_nuc(Mg24) = 0.0_rt;

ydot_nuc(Fe56) = 0.0_rt;

enuc_weak = rate_eval.enuc_weak;
}
#endif


AMREX_GPU_HOST_DEVICE AMREX_INLINE
void rhs_nuc(const burn_t& state,
amrex::Array1D<amrex::Real, 1, neqs>& ydot_nuc,
Expand Down Expand Up @@ -858,7 +934,7 @@ void rhs_nuc(const burn_t& state,
(screened_rates(k_He4_Ne20_to_Mg24)*Y(He4)*Y(Ne20)*state.rho + -screened_rates(k_Mg24_to_He4_Ne20)*Y(Mg24)) +
(screened_rates(k_C12_O16_to_He4_Mg24)*Y(C12)*Y(O16)*state.rho + -screened_rates(k_He4_Mg24_to_C12_O16)*Y(He4)*Y(Mg24)*state.rho);

ydot_nuc(Fe56) = 0.0;
ydot_nuc(Fe56) = 0.0_rt;

}

Expand Down
83 changes: 82 additions & 1 deletion networks/ECSN/actual_rhs.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {

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

rate_eval.enuc_weak = 0.0;
rate_eval.enuc_weak = 0.0_rt;

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);
Expand Down Expand Up @@ -312,6 +312,87 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {

}

#ifdef NSE_NET
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void get_ydot_weak(const burn_t& state,
amrex::Array1D<amrex::Real, 1, neqs>& ydot_nuc,
amrex::Real& enuc_weak,
[[maybe_unused]] const amrex::Array1D<amrex::Real, 1, NumSpec>& Y) {
///
/// Calculate Ydots contribute only from weak reactions.
/// This is used to calculate dyedt and energy generation from
/// weak reactions for self-consistent NSE
///


// initialize ydot_nuc to 0

for (int i = 1; i <= neqs; ++i) {
ydot_nuc(i) = 0.0_rt;
}

rate_t rate_eval;

[[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma;
[[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e;

rate_eval.enuc_weak = 0.0_rt;

// Calculate tabular rates and get ydot_weak

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;
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;
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;
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;
rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);

auto screened_rates = rate_eval.screened_rates;

ydot_nuc(H1) = 0.0_rt;

ydot_nuc(He4) = 0.0_rt;

ydot_nuc(O16) = 0.0_rt;

ydot_nuc(O20) =
(-screened_rates(k_O20_to_F20)*Y(O20) + screened_rates(k_F20_to_O20)*Y(F20));

ydot_nuc(F20) =
(screened_rates(k_O20_to_F20)*Y(O20) + -screened_rates(k_F20_to_O20)*Y(F20)) +
(-screened_rates(k_F20_to_Ne20)*Y(F20) + screened_rates(k_Ne20_to_F20)*Y(Ne20));

ydot_nuc(Ne20) =
(screened_rates(k_F20_to_Ne20)*Y(F20) + -screened_rates(k_Ne20_to_F20)*Y(Ne20));

ydot_nuc(Mg24) = 0.0_rt;

ydot_nuc(Al27) = 0.0_rt;

ydot_nuc(Si28) = 0.0_rt;

ydot_nuc(P31) = 0.0_rt;

ydot_nuc(S32) = 0.0_rt;

enuc_weak = rate_eval.enuc_weak;
}
#endif


AMREX_GPU_HOST_DEVICE AMREX_INLINE
void rhs_nuc(const burn_t& state,
amrex::Array1D<amrex::Real, 1, neqs>& ydot_nuc,
Expand Down
183 changes: 182 additions & 1 deletion networks/He-C-Fe-group/actual_rhs.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1276,7 +1276,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {

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

rate_eval.enuc_weak = 0.0;
rate_eval.enuc_weak = 0.0_rt;

tabular_evaluate(j_Co55_Fe55_meta, j_Co55_Fe55_rhoy, j_Co55_Fe55_temp, j_Co55_Fe55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
Expand Down Expand Up @@ -1377,6 +1377,187 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {

}

#ifdef NSE_NET
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void get_ydot_weak(const burn_t& state,
amrex::Array1D<amrex::Real, 1, neqs>& ydot_nuc,
amrex::Real& enuc_weak,
[[maybe_unused]] const amrex::Array1D<amrex::Real, 1, NumSpec>& Y) {
///
/// Calculate Ydots contribute only from weak reactions.
/// This is used to calculate dyedt and energy generation from
/// weak reactions for self-consistent NSE
///


// initialize ydot_nuc to 0

for (int i = 1; i <= neqs; ++i) {
ydot_nuc(i) = 0.0_rt;
}

rate_t rate_eval;

[[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma;
[[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e;

rate_eval.enuc_weak = 0.0_rt;

// Calculate tabular rates and get ydot_weak

tabular_evaluate(j_Co55_Fe55_meta, j_Co55_Fe55_rhoy, j_Co55_Fe55_temp, j_Co55_Fe55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Co55_to_Fe55) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Co55) * (edot_nu + edot_gamma);

tabular_evaluate(j_Co56_Fe56_meta, j_Co56_Fe56_rhoy, j_Co56_Fe56_temp, j_Co56_Fe56_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Co56_to_Fe56) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Co56) * (edot_nu + edot_gamma);

tabular_evaluate(j_Co56_Ni56_meta, j_Co56_Ni56_rhoy, j_Co56_Ni56_temp, j_Co56_Ni56_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Co56_to_Ni56) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Co56) * (edot_nu + edot_gamma);

tabular_evaluate(j_Co57_Ni57_meta, j_Co57_Ni57_rhoy, j_Co57_Ni57_temp, j_Co57_Ni57_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Co57_to_Ni57) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Co57) * (edot_nu + edot_gamma);

tabular_evaluate(j_Fe55_Co55_meta, j_Fe55_Co55_rhoy, j_Fe55_Co55_temp, j_Fe55_Co55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Fe55_to_Co55) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe55) * (edot_nu + edot_gamma);

tabular_evaluate(j_Fe55_Mn55_meta, j_Fe55_Mn55_rhoy, j_Fe55_Mn55_temp, j_Fe55_Mn55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Fe55_to_Mn55) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe55) * (edot_nu + edot_gamma);

tabular_evaluate(j_Fe56_Co56_meta, j_Fe56_Co56_rhoy, j_Fe56_Co56_temp, j_Fe56_Co56_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Fe56_to_Co56) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe56) * (edot_nu + edot_gamma);

tabular_evaluate(j_Mn55_Fe55_meta, j_Mn55_Fe55_rhoy, j_Mn55_Fe55_temp, j_Mn55_Fe55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Mn55_to_Fe55) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Mn55) * (edot_nu + edot_gamma);

tabular_evaluate(j_n_p_meta, j_n_p_rhoy, j_n_p_temp, j_n_p_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_n_to_p) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(N) * (edot_nu + edot_gamma);

tabular_evaluate(j_Ni56_Co56_meta, j_Ni56_Co56_rhoy, j_Ni56_Co56_temp, j_Ni56_Co56_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Ni56_to_Co56) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Ni56) * (edot_nu + edot_gamma);

tabular_evaluate(j_Ni57_Co57_meta, j_Ni57_Co57_rhoy, j_Ni57_Co57_temp, j_Ni57_Co57_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Ni57_to_Co57) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Ni57) * (edot_nu + edot_gamma);

tabular_evaluate(j_p_n_meta, j_p_n_rhoy, j_p_n_temp, j_p_n_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_p_to_n) = rate;
rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma);

auto screened_rates = rate_eval.screened_rates;

ydot_nuc(N) =
-screened_rates(k_n_to_p)*Y(N) +
screened_rates(k_p_to_n)*Y(H1);

ydot_nuc(H1) =
screened_rates(k_n_to_p)*Y(N) +
-screened_rates(k_p_to_n)*Y(H1);

ydot_nuc(He4) = 0.0_rt;

ydot_nuc(C12) = 0.0_rt;

ydot_nuc(N13) = 0.0_rt;

ydot_nuc(N14) = 0.0_rt;

ydot_nuc(O16) = 0.0_rt;

ydot_nuc(F18) = 0.0_rt;

ydot_nuc(Ne20) = 0.0_rt;

ydot_nuc(Ne21) = 0.0_rt;

ydot_nuc(Na22) = 0.0_rt;

ydot_nuc(Na23) = 0.0_rt;

ydot_nuc(Mg24) = 0.0_rt;

ydot_nuc(Al27) = 0.0_rt;

ydot_nuc(Si28) = 0.0_rt;

ydot_nuc(P31) = 0.0_rt;

ydot_nuc(S32) = 0.0_rt;

ydot_nuc(Ar36) = 0.0_rt;

ydot_nuc(Ca40) = 0.0_rt;

ydot_nuc(Ti44) = 0.0_rt;

ydot_nuc(Cr48) = 0.0_rt;

ydot_nuc(Mn51) = 0.0_rt;

ydot_nuc(Mn55) =
(screened_rates(k_Fe55_to_Mn55)*Y(Fe55) + -screened_rates(k_Mn55_to_Fe55)*Y(Mn55));

ydot_nuc(Fe52) = 0.0_rt;

ydot_nuc(Fe53) = 0.0_rt;

ydot_nuc(Fe54) = 0.0_rt;

ydot_nuc(Fe55) =
(screened_rates(k_Co55_to_Fe55)*Y(Co55) + -screened_rates(k_Fe55_to_Co55)*Y(Fe55)) +
(-screened_rates(k_Fe55_to_Mn55)*Y(Fe55) + screened_rates(k_Mn55_to_Fe55)*Y(Mn55));

ydot_nuc(Fe56) =
(screened_rates(k_Co56_to_Fe56)*Y(Co56) + -screened_rates(k_Fe56_to_Co56)*Y(Fe56));

ydot_nuc(Co55) =
(-screened_rates(k_Co55_to_Fe55)*Y(Co55) + screened_rates(k_Fe55_to_Co55)*Y(Fe55));

ydot_nuc(Co56) =
(-screened_rates(k_Co56_to_Fe56)*Y(Co56) + screened_rates(k_Fe56_to_Co56)*Y(Fe56)) +
(screened_rates(k_Ni56_to_Co56)*Y(Ni56) + -screened_rates(k_Co56_to_Ni56)*Y(Co56));

ydot_nuc(Co57) =
(screened_rates(k_Ni57_to_Co57)*Y(Ni57) + -screened_rates(k_Co57_to_Ni57)*Y(Co57));

ydot_nuc(Ni56) =
(-screened_rates(k_Ni56_to_Co56)*Y(Ni56) + screened_rates(k_Co56_to_Ni56)*Y(Co56));

ydot_nuc(Ni57) =
(-screened_rates(k_Ni57_to_Co57)*Y(Ni57) + screened_rates(k_Co57_to_Ni57)*Y(Co57));

ydot_nuc(Ni58) = 0.0_rt;

ydot_nuc(Cu59) = 0.0_rt;

ydot_nuc(Zn60) = 0.0_rt;

enuc_weak = rate_eval.enuc_weak;
}
#endif


AMREX_GPU_HOST_DEVICE AMREX_INLINE
void rhs_nuc(const burn_t& state,
amrex::Array1D<amrex::Real, 1, neqs>& ydot_nuc,
Expand Down
Loading

0 comments on commit a5408ba

Please sign in to comment.