update self-consistent NSE SDC update with Tabular NSE

networks/CNO_extras/actual_rhs.H

@@ -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,
@@ -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;
 
 }
 

networks/ECSN/actual_rhs.H

@@ -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);
@@ -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,

networks/He-C-Fe-group/actual_rhs.H

@@ -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);
@@ -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,