add support for self-consistent NSE to test_sdc

unit_test/test_sdc/_parameters

@@ -6,3 +6,6 @@ temp_min      real       1.e6
 temp_max      real       1.e15
 
 tmax          real       0.1e0
+
+mu_p          real       -6.0
+mu_n          real      -11.0

unit_test/test_sdc/inputs_aprox19.no_nse

@@ -1,7 +1,5 @@
 n_cell = 16
 
-do_cxx = 1
-
 prefix = react_aprox19_sdc_
 
 unit_test.dens_min   = 1.e4

unit_test/test_sdc/inputs_aprox19.nse

@@ -1,7 +1,5 @@
 n_cell = 16
 
-do_cxx = 1
-
 prefix = react_aprox19_sdc_
 
 unit_test.dens_min   = 1.e5

unit_test/test_sdc/inputs_ase

@@ -0,0 +1,24 @@
+n_cell = 16
+
+prefix = react_aprox19_sdc_
+
+unit_test.dens_min   = 1.e7
+unit_test.dens_max   = 5.e8
+unit_test.temp_min   = 1.e8
+unit_test.temp_max   = 6.e9
+
+unit_test.small_dens = 1.e3
+
+unit_test.tmax = 1.e-6
+
+unit_test.primary_species_1 = helium-4
+unit_test.primary_species_2 = nickel-56
+
+integrator.jacobian = 1
+
+nse.nse_molar_independent = 0
+nse.nse_dx_independent = 1
+nse.ase_tol = 0.1
+nse.nse_skip_molar = 0
+nse.solve_nse_e_mode = 1
+

unit_test/test_sdc/main.cpp

@@ -176,8 +176,9 @@ void main_main ()
     }
 
     // allocate a multifab for the number of RHS calls
+    // and whether we are in NSE
     // so we can manually do the reductions (for GPU)
-    iMultiFab integrator_n_rhs(ba, dm, 1, Nghost);
+    iMultiFab integrator_n_rhs(ba, dm, 2, Nghost);
 
     // What time is it now?  We'll use this to compute total react time.
     Real strt_time = ParallelDescriptor::second();
@@ -226,6 +227,10 @@ void main_main ()
     int n_rhs_max = integrator_n_rhs.max(0);
     long n_rhs_sum = integrator_n_rhs.sum(0, 0, true);
 
+#ifdef NSE
+    long n_nse_sum = integrator_n_rhs.sum(1, 0, true);
+#endif
+
     // get the name of the integrator from the build info functions
     // written at compile time.  We will append the name of the
     // integrator to the output file name
@@ -250,10 +255,16 @@ void main_main ()
     // Tell the I/O Processor to write out the "run time"
     amrex::Print() << "Run time = " << stop_time << std::endl;
 
+    long n_cell_cubed = static_cast<long>(n_cell) * n_cell * n_cell;
+
     // print statistics
     std::cout << "min number of rhs calls: " << n_rhs_min << std::endl;
     std::cout << "avg number of rhs calls: "
-              << static_cast<Real>(n_rhs_sum) / static_cast<Real>(n_cell*n_cell*n_cell) << std::endl;
+              << n_rhs_sum / n_cell_cubed << std::endl;
     std::cout << "max number of rhs calls: " << n_rhs_max << std::endl;
+    std::cout << std::endl;
+#ifdef NSE
+    std::cout << "fraction of zones in NSE: " << static_cast<amrex::Real>(n_nse_sum) / n_cell_cubed << std::endl;
+#endif
 
 }

unit_test/test_sdc/react_zones.H

@@ -31,6 +31,13 @@ bool do_react (const plot_t& p, const int i, const int j, const int k,
 #ifdef AUX_THERMO
     set_aux_comp_from_X(eos_state);
 #endif
+#ifdef NSE_NET
+    // if we are doing self-consistent NSE, then set Ye manually
+    eos_state.y_e = 0.0_rt;
+    for (int n = 0; n < NumSpec; ++n) {
+        eos_state.y_e = eos_state.xn[n] * zion[n] * aion_inv[n];
+    }
+#endif
 
     eos(eos_input_rt, eos_state);
 
@@ -49,6 +56,12 @@ bool do_react (const plot_t& p, const int i, const int j, const int k,
         burn_state.y[SFX+n] = burn_state.rho * eos_state.aux[n];
     }
 #endif
+#ifdef NSE_NET
+    burn_state.y_e = eos_state.y_e;
+
+    burn_state.mu_p = unit_test_rp::mu_p;
+    burn_state.mu_n = unit_test_rp::mu_n;
+#endif
 
     burn_state.y[SRHO] = state(i, j, k, p.irho);
 
@@ -109,9 +122,14 @@ bool do_react (const plot_t& p, const int i, const int j, const int k,
     }
 #endif
 
+#ifdef NSE
+    state(i, j, k, p.inse) = burn_state.nse;
+#endif
+
     state(i, j, k, p.irho_Hnuc) = (burn_state.y[SEINT] - rhoe_old) / unit_test_rp::tmax;
 
     n_rhs(i, j, k, 0) = burn_state.n_rhs;
+    n_rhs(i, j, k, 1) = static_cast<int>(burn_state.nse);
 
     return burn_state.success;
 

unit_test/test_sdc/variables.H

@@ -10,25 +10,26 @@ class plot_t {
 
 public:
 
-
-  int irho = -1;
-  int itemp = -1;
-  int ispec = -1;
-  int ispec_old = -1;
+    int irho = -1;
+    int itemp = -1;
+    int ispec = -1;
+    int ispec_old = -1;
 #if NAUX_NET > 0
-  int iaux = -1;
-  int iaux_old = -1;
+    int iaux = -1;
+    int iaux_old = -1;
 #endif
-  int irodot = -1;
-  int irho_Hnuc = -1;
-
-  int n_plot_comps = 0;
+    int irodot = -1;
+    int irho_Hnuc = -1;
+#ifdef NSE
+    int inse = -1;
+#endif
+    int n_plot_comps = 0;
 
-  int next_index(const int num) {
-    int next = n_plot_comps;
-    n_plot_comps += num;
-    return next;
-  }
+    int next_index(const int num) {
+        int next = n_plot_comps;
+        n_plot_comps += num;
+        return next;
+    }
 
 };
 

unit_test/test_sdc/variables.cpp

@@ -19,6 +19,10 @@ plot_t init_variables() {
   p.irodot = p.next_index(NumSpec);
   p.irho_Hnuc = p.next_index(1);
 
+#ifdef NSE
+  p.inse = p.next_index(1);
+#endif
+
   return p;
 }
 
@@ -41,4 +45,9 @@ void get_varnames(const plot_t& p, amrex::Vector<std::string>& names) {
       names[p.iaux_old + n] = "old_" + aux_names_cxx[n];
   }
 #endif
+
+#ifdef NSE
+  names[p.inse] = "in_nse";
+#endif
+
 }