Update rhoX for nse_solver for SDC and discard strang (#1607)

This is to address #1606 and partially addresses #1542 Also some clean-ups to nse-solver. It mainly stores state.y[SFS+n] at the end.
AMReX-Astro · Jul 11, 2024 · 4142d4f · 4142d4f
1 parent 0b1a231
commit 4142d4f
Show file tree

Hide file tree

Showing 11 changed files with 410 additions and 407 deletions.
diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H
@@ -20,7 +20,9 @@
 // First check to see if we're in the ballpark of nse state
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) {
+void check_nse_molar(const amrex::Array1D<amrex::Real, 1, NumSpec>& Y,
+                     const amrex::Array1D<amrex::Real, 1, NumSpec>& Y_nse,
+                     bool& nse_check) {
 
     // This function gives the first estimate whether we're in the nse or not
     // it checks whether the molar fractions of n,p,a are approximately in NSE
@@ -44,13 +46,13 @@ void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) {
     for (int n = 0; n < NumSpec; ++n) {
 
         if (n == NSE_INDEX::H1_index || n == NSE_INDEX::N_index) {
-            r /= state.xn[n] * state.xn[n] * aion_inv[n] * aion_inv[n];
-            r_nse /= nse_state.xn[n] * nse_state.xn[n] * aion_inv[n] * aion_inv[n];
+            r /= Y(n+1) * Y(n+1);
+            r_nse /= Y_nse(n+1) * Y_nse(n+1);
         }
 
         else if (n == NSE_INDEX::He4_index) {
-            r *= state.xn[n] * aion_inv[n];
-            r_nse *= nse_state.xn[n] * aion_inv[n];
+            r *= Y(n+1);
+            r_nse *= Y_nse(n+1);
         }
     }
 
@@ -59,15 +61,15 @@ void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) {
 
     // if there is neutron in the network
 
-    if ((std::abs(r - r_nse) < 0.5_rt*r_nse)
+    if ((std::abs(r - r_nse) < 0.5_rt * r_nse)
         && (NSE_INDEX::N_index != -1)) {
         nse_check = true;
         return;
     }
 
     // if there is no neutron in the network
 
-    if ((std::abs(r - r_nse) < 0.25_rt*r_nse)
+    if ((std::abs(r - r_nse) < 0.25_rt * r_nse)
         && (NSE_INDEX::N_index == -1)) {
         nse_check = true;
         return;
@@ -76,9 +78,8 @@ void check_nse_molar(burn_t& state, const burn_t& nse_state, bool& nse_check) {
     // Overall molar fraction check
 
     for (int n = 0; n < NumSpec; ++n) {
-        amrex::Real abs_diff = std::abs(state.xn[n] - nse_state.xn[n]) / aion[n];
-        amrex::Real rel_diff = abs_diff / (state.xn[n] / aion[n]);
-        if (abs_diff > nse_abs_tol && rel_diff > nse_rel_tol) {
+        amrex::Real abs_diff = std::abs(Y(n+1) - Y_nse(n+1));
+        if (abs_diff > nse_abs_tol && abs_diff > nse_rel_tol * Y(n+1)) {
             return;
         }
     }
@@ -228,7 +229,7 @@ bool in_single_group(const amrex::Array1D<int, 1, NumSpec>& group_ind) {
 template <typename T>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void fill_reaction_timescale(amrex::Array1D<T, 1, Rates::NumRates>& reaction_timescales,
-                             const int current_rate_index, const burn_t& state,
+                             const int current_rate_index, const amrex::Real rho,
                              const amrex::Array1D<amrex::Real, 1, NumSpec>& Y,
                              const amrex::Array1D<amrex::Real, 1, Rates::NumRates>& screened_rates,
                              const amrex::Real t_s) {
@@ -301,14 +302,14 @@ void fill_reaction_timescale(amrex::Array1D<T, 1, Rates::NumRates>& reaction_tim
         if (NSE_INDEX::rate_indices(current_rate_index, 2) == NSE_INDEX::rate_indices(current_rate_index, 3)) {
             b_f *= 0.5_rt;
         }
-        b_f *= Y(NSE_INDEX::rate_indices(current_rate_index, 2) + 1) * state.rho;
+        b_f *= Y(NSE_INDEX::rate_indices(current_rate_index, 2) + 1) * rho;
     }
 
     if (NSE_INDEX::rate_indices(current_rate_index, 5) != -1) {
         if (NSE_INDEX::rate_indices(current_rate_index, 5) == NSE_INDEX::rate_indices(current_rate_index, 6)) {
             b_r *= 0.5_rt;
         }
-        b_r *= Y(NSE_INDEX::rate_indices(current_rate_index, 5) + 1) * state.rho;
+        b_r *= Y(NSE_INDEX::rate_indices(current_rate_index, 5) + 1) * rho;
     }
 
     // Find the timescale of the rate, See Equation 11 in Kushnir
@@ -437,7 +438,7 @@ void fill_merge_indices(amrex::Array1D<int, 1, 2>& merge_indices,
 
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
-void nse_grouping(amrex::Array1D<int, 1, NumSpec>& group_ind, const burn_t& state,
+void nse_grouping(amrex::Array1D<int, 1, NumSpec>& group_ind, const amrex::Real rho,
                   const amrex::Array1D<amrex::Real, 1, NumSpec>& Y,
                   const amrex::Array1D<amrex::Real, 1, Rates::NumRates>& screened_rates,
                   const amrex::Real t_s) {
@@ -468,7 +469,7 @@ void nse_grouping(amrex::Array1D<int, 1, NumSpec>& group_ind, const burn_t& stat
     amrex::Array1D<int, 1, Rates::NumRates> rate_indices;
 
     for (int n = 1; n <= Rates::NumRates; ++n) {
-        fill_reaction_timescale(reaction_timescales, n, state, Y,
+        fill_reaction_timescale(reaction_timescales, n, rho, Y,
                                 screened_rates, t_s);
         rate_indices(n) = n;
     }
@@ -543,19 +544,22 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) {
 
     const auto nse_state = get_actual_nse_state(current_state);
 
-    burn_t state = current_state;
+    auto state = current_state;
+
+    // Convert to molar fractions
+
+    amrex::Array1D<amrex::Real, 1, NumSpec> Y;
+    amrex::Array1D<amrex::Real, 1, NumSpec> Y_nse;
 
-#ifndef STRANG
-    // if not strang, store mass fractions
     for (int n = 0; n < NumSpec; ++n) {
-        state.xn[n] = current_state.y[SFS+n] / current_state.rho;
+        Y(n+1) = current_state.y[SFS+n] * aion_inv[n] / current_state.rho;
+        Y_nse(n+1) = nse_state.xn[n] * aion_inv[n];
     }
-#endif
 
     // Check whether state is in the ballpark of NSE
 
     if (!skip_molar_check) {
-        check_nse_molar(state, nse_state, current_state.nse);
+        check_nse_molar(Y, Y_nse, current_state.nse);
         if (!current_state.nse) {
             return current_state.nse;
         }
@@ -577,23 +581,16 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) {
     // our current mass fractions are in the ballpark of NSE mass fractions.
 
     if (nse_molar_independent || skip_molar_check) {
-        state = nse_state;
-        state.dx = current_state.dx;
-#ifndef STRANG
-        for (int n = 0; n < NumSpec; ++n) {
-            state.y[SFS+n] = nse_state.xn[n] * nse_state.rho;
+        for (int n = 1; n <= NumSpec; ++n) {
+            Y(n) = Y_nse(n);
         }
-#endif
     }
 
-    // set molar fractions
-
-    amrex::Array1D<amrex::Real, 1, NumSpec> Y;
+    // store xn to find the rates
 
-    for (int n = 1; n <= NumSpec; ++n) {
-        Y(n) = state.xn[n-1] * aion_inv[n-1];
+    for (int n = 0; n < NumSpec; ++n) {
+        state.xn[n] = Y(n+1) * aion[n];
     }
-
     rate_t rate_eval;
     constexpr int do_T_derivatives = 0;
     evaluate_rates<do_T_derivatives, rate_t>(state, rate_eval);
@@ -629,7 +626,7 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) {
 
     // Now do nse grouping
 
-    nse_grouping(group_ind, state, Y, rate_eval.screened_rates, t_s);
+    nse_grouping(group_ind, state.rho, Y, rate_eval.screened_rates, t_s);
 
     // Check if we result in a single group after grouping
 

diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H
@@ -51,6 +51,7 @@ T get_nonexponent_nse_state(const T& state) {
     // the temperature that comes in through T_fixed
 
     amrex::Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T;
+    constexpr amrex::Real ihplanck2 = 1.0_rt / (C::hplanck * C::hplanck);
 
 #ifndef NEW_NETWORK_IMPLEMENTATION
     auto tfactors = evaluate_tfactors(T_in);
@@ -59,7 +60,7 @@ T get_nonexponent_nse_state(const T& state) {
     for (int n = 0; n < NumSpec; ++n) {
 #ifdef NEW_NETWORK_IMPLEMENTATION
         if (n == NSE_INDEX::H1_index) {
-            nse_state.xn[n] = 0.0;
+            nse_state.xn[n] = 0.0_rt;
             continue;
         }
 #endif
@@ -73,11 +74,11 @@ T get_nonexponent_nse_state(const T& state) {
 
         // find nse mass frac without the exponent term.
 
-        amrex::Real power = 2.0 * M_PI * network::mion(n+1) *
-            C::k_B * T_in / (C::hplanck * C::hplanck);
+        amrex::Real power = 2.0_rt * M_PI * network::mion(n+1) *
+            C::k_B * T_in * ihplanck2;
 
         nse_state.xn[n] = network::mion(n+1) * pf * spin / state.rho *
-            std::sqrt(power * power * power);
+            std::sqrt(amrex::Math::powi<3>(power));
     }
 
     return nse_state;
@@ -103,9 +104,9 @@ void compute_coulomb_contribution(amrex::Array1D<amrex::Real, 1, NumSpec>& u_c,
     // so we just treat u_c as a constant.
     //
 
-    const amrex::Real n_e = state.rho * state.y_e / C::m_u;
+    const amrex::Real n_e = state.rho * state.y_e * C::Legacy::n_A;
     const amrex::Real Gamma_e = C::q_e * C::q_e *
-        std::cbrt(4.0_rt * M_PI * n_e / 3.0_rt) / (C::k_B * T_in);
+        std::cbrt(1.333333333333_rt * M_PI * n_e) / (C::k_B * T_in);
     amrex::Real gamma;
 
     for (int n = 0; n < NumSpec; ++n) {
@@ -116,7 +117,7 @@ void compute_coulomb_contribution(amrex::Array1D<amrex::Real, 1, NumSpec>& u_c,
 #endif
         // term for calculating u_c
 
-        gamma = std::pow(zion[n], 5.0_rt/3.0_rt) * Gamma_e;
+        gamma = std::cbrt(amrex::Math::powi<5>(zion[n])) * Gamma_e;
 
         // chemical potential for coulomb correction
         // see appendix of Calder 2007, doi:10.1086/510709 for more detail
@@ -146,12 +147,13 @@ void apply_nse_exponent(T& nse_state,
     // the temperature that comes in through T_fixed
 
     amrex::Real T_in = nse_state.T_fixed > 0.0_rt ? nse_state.T_fixed : nse_state.T;
+    amrex::Real ikTMeV = C::Legacy::MeV2erg / (C::k_B * T_in);
     amrex::Real exponent;
 
     for (int n = 0; n < NumSpec; ++n) {
 #ifdef NEW_NETWORK_IMPLEMENTATION
         if (n == NSE_INDEX::H1_index) {
-            nse_state.xn[n] = 0.0;
+            nse_state.xn[n] = 0.0_rt;
             continue;
         }
 #endif
@@ -163,8 +165,8 @@ void apply_nse_exponent(T& nse_state,
 
         exponent = amrex::min(500.0_rt,
                               (zion[n] * nse_state.mu_p + (aion[n] - zion[n]) *
-                               nse_state.mu_n + network::bion(n+1)) /
-                              C::k_B / T_in * C::Legacy::MeV2erg - u_c(n+1));
+                               nse_state.mu_n + network::bion(n+1)) *
+                              ikTMeV - u_c(n+1));
 
         nse_state.xn[n] *= std::exp(exponent);
     }
@@ -210,6 +212,7 @@ void fcn(amrex::Array1D<amrex::Real, 1, 2>& x, amrex::Array1D<amrex::Real, 1, 2>
         }
 #endif
         // constraint equation 1, mass fraction sum to 1
+        // since we're using rhoX, make it sum to rho
 
         fvec(1) += nse_state.xn[n];
     }
@@ -247,17 +250,18 @@ void jcn(amrex::Array1D<amrex::Real, 1, 2>& x, amrex::Array2D<amrex::Real, 1, 2,
     // the temperature that comes in through T_fixed
 
     amrex::Real T_in = nse_state.T_fixed > 0.0_rt ? nse_state.T_fixed : nse_state.T;
+    amrex::Real ikTMeV = C::Legacy::MeV2erg / (C::k_B * T_in);
 
     for (int n = 0; n < NumSpec; ++n) {
 #ifdef NEW_NETWORK_IMPLEMENTATION
         if (n == NSE_INDEX::H1_index) {
             continue;
         }
 #endif
-        fjac(1, 1) += nse_state.xn[n] * zion[n] / C::k_B / T_in * C::Legacy::MeV2erg ;
-        fjac(1, 2) += nse_state.xn[n] * (aion[n] - zion[n]) / C::k_B / T_in * C::Legacy::MeV2erg;
-        fjac(2, 1) += nse_state.xn[n] * zion[n] * zion[n] * aion_inv[n] / C::k_B / T_in * C::Legacy::MeV2erg;
-        fjac(2, 2) += nse_state.xn[n] * zion[n] * (aion[n] - zion[n]) * aion_inv[n] / C::k_B / T_in * C::Legacy::MeV2erg;
+        fjac(1, 1) += nse_state.xn[n] * zion[n] * ikTMeV;
+        fjac(1, 2) += nse_state.xn[n] * (aion[n] - zion[n]) * ikTMeV;
+        fjac(2, 1) += nse_state.xn[n] * zion[n] * zion[n] * aion_inv[n] * ikTMeV;
+        fjac(2, 2) += nse_state.xn[n] * zion[n] * (aion[n] - zion[n]) * aion_inv[n] * ikTMeV;
     }
 
 }
@@ -499,15 +503,13 @@ T get_actual_nse_state(T& state, amrex::Real eps=1.0e-10_rt,
 
     if (!input_ye_is_valid) {
         // ensure Ye is valid
-#ifdef STRANG
-        composition(state);
-#else
-        state.y_e = 0.0_rt;
+        // We assume the input state has rhoX
 
+        state.y_e = 0.0_rt;
         for (int n = 0; n < NumSpec; ++n) {
-            state.y_e += state.y[SFS+n] * zion[n] * aion_inv[n] / state.rho;
+            state.y_e += state.y[SFS+n] * zion[n] * aion_inv[n];
         }
-#endif
+        state.y_e /= state.rho;
     }
 
     nse_solver_data<T> state_data = {state, {0.0_rt}};
@@ -557,6 +559,12 @@ T get_actual_nse_state(T& state, amrex::Real eps=1.0e-10_rt,
     state.mu_p = state_data.state.mu_p;
     state.mu_n = state_data.state.mu_n;
 
+    // update rhoX in the output nse_state
+
+    for (int n = 0; n < NumSpec; ++n) {
+        state_data.state.y[SFS+n] = state.rho * state_data.state.xn[n];
+    }
+
     return state_data.state;
 }
 #endif
diff --git a/unit_test/test_ase/GNUmakefile b/unit_test/test_ase/GNUmakefile
@@ -28,10 +28,11 @@ NETWORK_DIR := ase
 CONDUCTIVITY_DIR := stellar
 
 SCREEN_METHOD := chabrier1998
-
 #SCREEN_METHOD := null
 #SCREEN_METHOD := chugunov2007
 
+USE_SIMPLIFIED_SDC = TRUE
+
 INTEGRATOR_DIR =  VODE
 
 EXTERN_SEARCH += .
@@ -40,5 +41,3 @@ Bpack   := ./Make.package
 Blocs   := .
 
 include $(MICROPHYSICS_HOME)/unit_test/Make.unit_test
-
-
diff --git a/unit_test/test_ase/burn_cell.H b/unit_test/test_ase/burn_cell.H
@@ -50,12 +50,12 @@ void burn_cell_c()
       std::cout << short_spec_names_cxx[n] << " : " << NSE_STATE.xn[n] << std::endl;
     }
 
-    // Let state.xn equal to nse_state.xn to make sure its in nse state
+    // Let state.y equal to nse_state.y to make sure its in nse state
+
     for (int n = 0; n < NumSpec; ++n){
-      state.xn[n] = NSE_STATE.xn[n];
+      state.y[SFS+n] = NSE_STATE.y[SFS+n];
     }
 
-
     // get eos
     //eos(eos_input_rt, state);
 

diff --git a/unit_test/test_ase/ci-benchmarks/ase_nse_net_unit_test.out b/unit_test/test_ase/ci-benchmarks/ase_nse_net_unit_test.out
@@ -8,28 +8,28 @@ State Density (g/cm^3): 10000000
 State Temperature (K): 6000000000
 electron fraction is 0.5
 NSE state: 
-n : 0.0006075693495
-H1 : 0.001992432534
-He4 : 0.519007877
-C12 : 1.326102e-05
-N13 : 9.477376691e-11
-N14 : 2.565486854e-09
-O16 : 3.225455246e-05
-F18 : 5.422064179e-11
-Ne20 : 7.427509432e-07
-Ne21 : 4.998912611e-08
-Na22 : 2.608008186e-09
-Na23 : 1.00884295e-06
-Mg24 : 0.0001022222509
-Al27 : 0.0005548826609
-Si28 : 0.03680966846
-P31 : 0.04229223791
-S32 : 0.03865771985
-Ar36 : 0.02464742443
-Ca40 : 0.02885152622
-Ti44 : 0.00166123105
-Cr48 : 0.009967017598
-Fe52 : 0.05968005236
-Ni56 : 0.2351208159
+n : 0.000607569351
+H1 : 0.001992432526
+He4 : 0.5190078752
+C12 : 1.326101988e-05
+N13 : 9.477376579e-11
+N14 : 2.56548683e-09
+O16 : 3.225455214e-05
+F18 : 5.422064122e-11
+Ne20 : 7.427509352e-07
+Ne21 : 4.99891257e-08
+Na22 : 2.608008157e-09
+Na23 : 1.008842941e-06
+Mg24 : 0.0001022222498
+Al27 : 0.0005548826564
+Si28 : 0.03680966808
+P31 : 0.0422922376
+S32 : 0.0386577195
+Ar36 : 0.02464742425
+Ca40 : 0.0288515261
+Ti44 : 0.001661231048
+Cr48 : 0.009967017622
+Fe52 : 0.05968005276
+Ni56 : 0.2351208186
 We're in NSE. 
 AMReX (23.07-395-ge6c93bf22695) finalized