Merge branch 'development' into ROCK4

integration/_parameters

@@ -77,3 +77,6 @@ sdc_burn_tol_factor          real         1.d0
 # for Strang, this simply means scaling e by the initial energy?
 scale_system              integer         0
 
+# for the NSE update predictor-corrector, how many iterations
+# do we take to get the new time NSE state
+nse_iters                 integer         3

integration/nse_update_simplified_sdc.H

@@ -93,7 +93,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
         rho_aux_new[n] = state.y[SFX+n] + dt * state.ydot_a[SFX+n] + dt * aux_source[n];
     }
 
-    for (int iter = 0; iter < 3; iter++) {
+    for (int iter = 0; iter < integrator_rp::nse_iters; iter++) {
 
         // update (rho e)^{n+1} based on the new energy generation rate
         rhoe_new = state.y[SEINT] + dt * state.ydot_a[SEINT] + dt * rho_enucdot;
@@ -258,7 +258,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) {
         rhoX_new[n] = state.y[SFS+n] + dt * state.ydot_a[SFS+n] + dt * rhoX_source[n];
     }
 
-    for (int iter = 0; iter < 3; iter++) {
+    for (int iter = 0; iter < integrator_rp::nse_iters; iter++) {
 
         // update (rho e)^{n+1} based on the new energy generation rate
         rhoe_new = state.y[SEINT] + dt * state.ydot_a[SEINT] + dt * rho_enucdot;

networks/ignition_reaclib/URCA-simple/actual_network.H

@@ -28,19 +28,20 @@ namespace Rates
         k_C12_C12_to_n_Mg23 = 2,
         k_C12_C12_to_p_Na23 = 3,
         k_He4_C12_to_O16 = 4,
-        k_n_to_p_weak_wc12 = 5,
-        k_Na23_to_Ne23 = 6,
-        k_Ne23_to_Na23 = 7,
-        NumRates = k_Ne23_to_Na23
+        k_Na23_to_Ne23 = 5,
+        k_Ne23_to_Na23 = 6,
+        k_n_to_p = 7,
+        k_p_to_n = 8,
+        NumRates = k_p_to_n
     };
 
     // number of reaclib rates
 
-    const int NrateReaclib = 5;
+    const int NrateReaclib = 4;
 
     // number of tabular rates
 
-    const int NrateTabular = 2;
+    const int NrateTabular = 4;
 
     // rate names -- note: the rates are 1-based, not zero-based, so we pad
     // this vector with rate_names[0] = "" so the indices line up with the
@@ -52,9 +53,10 @@ namespace Rates
         "C12_C12_to_n_Mg23",  // 2,
         "C12_C12_to_p_Na23",  // 3,
         "He4_C12_to_O16",  // 4,
-        "n_to_p_weak_wc12",  // 5,
-        "Na23_to_Ne23",  // 6,
-        "Ne23_to_Na23"  // 7,
+        "Na23_to_Ne23",  // 5,
+        "Ne23_to_Na23",  // 6,
+        "n_to_p",  // 7,
+        "p_to_n"  // 8,
     };
 
 }

networks/ignition_reaclib/URCA-simple/actual_network_data.cpp

@@ -15,9 +15,10 @@ namespace NSE_INDEX
         -1, 3, 3, -1, 0, 8, -1,
         -1, 3, 3, -1, 1, 7, -1,
         -1, 2, 3, -1, -1, 4, -1,
-        -1, -1, 0, -1, -1, 1, -1,
         -1, -1, 7, -1, -1, 6, -1,
-        -1, -1, 6, -1, -1, 7, 6
+        -1, -1, 6, -1, -1, 7, 5,
+        -1, -1, 0, -1, -1, 1, -1,
+        -1, -1, 1, -1, -1, 0, -1
     };
 }
 #endif

networks/ignition_reaclib/URCA-simple/actual_rhs.H

@@ -151,6 +151,22 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     }
     rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne23) * (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;
+    if constexpr (std::is_same<T, rate_derivs_t>::value) {
+        rate_eval.dscreened_rates_dT(k_n_to_p) = drate_dt;
+    }
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(N) * (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;
+    if constexpr (std::is_same<T, rate_derivs_t>::value) {
+        rate_eval.dscreened_rates_dT(k_p_to_n) = drate_dt;
+    }
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma);
+
 
 }
 
@@ -163,12 +179,14 @@ void rhs_nuc(const burn_t& state,
     using namespace Rates;
 
     ydot_nuc(N) =
-        -screened_rates(k_n_to_p_weak_wc12)*Y(N) +
+        -screened_rates(k_n_to_p)*Y(N) +
+        screened_rates(k_p_to_n)*Y(H1) +
         0.5*screened_rates(k_C12_C12_to_n_Mg23)*std::pow(Y(C12), 2)*state.rho;
 
     ydot_nuc(H1) =
         0.5*screened_rates(k_C12_C12_to_p_Na23)*std::pow(Y(C12), 2)*state.rho +
-        screened_rates(k_n_to_p_weak_wc12)*Y(N);
+        screened_rates(k_n_to_p)*Y(N) +
+        -screened_rates(k_p_to_n)*Y(H1);
 
     ydot_nuc(He4) =
         0.5*screened_rates(k_C12_C12_to_He4_Ne20)*std::pow(Y(C12), 2)*state.rho +
@@ -254,15 +272,21 @@ void jac_nuc(const burn_t& state,
 
     Real scratch;
 
-    scratch = -screened_rates(k_n_to_p_weak_wc12);
+    scratch = -screened_rates(k_n_to_p);
     jac.set(N, N, scratch);
 
+    scratch = screened_rates(k_p_to_n);
+    jac.set(N, H1, scratch);
+
     scratch = 1.0*screened_rates(k_C12_C12_to_n_Mg23)*Y(C12)*state.rho;
     jac.set(N, C12, scratch);
 
-    scratch = screened_rates(k_n_to_p_weak_wc12);
+    scratch = screened_rates(k_n_to_p);
     jac.set(H1, N, scratch);
 
+    scratch = -screened_rates(k_p_to_n);
+    jac.set(H1, H1, scratch);
+
     scratch = 1.0*screened_rates(k_C12_C12_to_p_Na23)*Y(C12)*state.rho;
     jac.set(H1, C12, scratch);