fix the CJ detonation utility (#1699)

util/cj_detonation/README.md

@@ -1,9 +1,23 @@
-Compute the Chapman-Jouguet detonation speed.  At the moment you set
-the state in the source directly (`main.f90`) and then build and
-execute as:
+Compute the Chapman-Jouguet detonation speed.
+
+The thermodynamic state for the fuel and ash are hardcoded
+into main.cpp right now.
+
+You can set the network and EOS in the `GNUmakefile`.
+
+At the moment, this seems to work well with the `aprox13` network.
+It does not seem as robust with larger nets, so some work needs to
+be done there.
+
 
 ```
-main.Linux.gfortran.exe inputs
+main3d.gnu.ex inputs
 ```
 
-You can set the network and EOS in the `GNUmakefile`.
+This will solve for the CJ speed and then try to compute points on
+the Hugoniot curve -- it will eventually give an error saying it
+doesn't converge.
+
+The output (`hugoniot.txt`) can then be plotted using the `cj_plot.py`
+script.
+

util/cj_detonation/_parameters

@@ -1,3 +1,5 @@
+@namespace: cj
+
 smallx     real   1.e-10
 
 small_temp    real       1.e5

util/cj_detonation/cj_det.H

@@ -90,7 +90,7 @@ cj_cond(const eos_t eos_state_fuel, eos_t& eos_state_ash, const Real q) {
             (1.0_rt + (eos_state_ash.p - eos_state_fuel.p) /
              (eos_state_ash.gam1 * eos_state_ash.p));
 
-        Real drho = eos_state_ash.rho - rho_old;
+        drho = eos_state_ash.rho - rho_old;
 
         if (std::abs(drho) < tol * eos_state_ash.rho) {
             converged = true;

util/cj_detonation/inputs

@@ -1 +1,2 @@
-amr.probin_file = probin
+cj.smallx = 1.e-10
+

util/cj_detonation/main.cpp

@@ -11,7 +11,10 @@ using namespace amrex;
 #include <network.H>
 #include <cj_det.H>
 #include <unit_test.H>
-#ifndef NEW_NETWORK_IMPLEMENTATION
+#include <actual_network.H>
+#ifdef NEW_NETWORK_IMPLEMENTATION
+#include <rhs.H>
+#else
 #include <actual_rhs.H>
 #endif
 
@@ -23,25 +26,11 @@ int main(int argc, char *argv[]) {
 
   std::cout << "starting the CJ Det solve..." << std::endl;
 
-  ParmParse ppa("amr");
-
-  std::string probin_file = "probin";
-
-  ppa.query("probin_file", probin_file);
-
-  std::cout << "probin = " << probin_file << std::endl;
-
-  const int probin_file_length = probin_file.length();
-  Vector<int> probin_file_name(probin_file_length);
-
-  for (int i = 0; i < probin_file_length; i++)
-    probin_file_name[i] = probin_file[i];
-
-  init_unit_test(probin_file_name.dataPtr(), &probin_file_length);
+  init_unit_test();
 
   // C++ EOS initialization (must be done after Fortran eos_init and
   // init_extern_parameters)
-  eos_init(small_temp, small_dens);
+  eos_init(cj_rp::small_temp, cj_rp::small_dens);
 
   // C++ Network, RHS, screening, rates initialization
 
@@ -60,9 +49,9 @@ int main(int argc, char *argv[]) {
   eos_state_fuel.rho = 1.e7_rt;
   eos_state_fuel.T = 1.e8_rt;
   for (int n = 0; n < NumSpec; n++) {
-      eos_state_fuel.xn[n] = smallx;
+      eos_state_fuel.xn[n] = cj_rp::smallx;
   }
-  eos_state_fuel.xn[0] = 1.0_rt - (NumSpec - 1) * smallx;
+  eos_state_fuel.xn[0] = 1.0_rt - (NumSpec - 1) * cj_rp::smallx;
 
   eos(eos_input_rt, eos_state_fuel);
 
@@ -71,9 +60,9 @@ int main(int argc, char *argv[]) {
 
   eos_t eos_state_ash = eos_state_fuel;
   for (int n = 0; n < NumSpec; n++) {
-      eos_state_ash.xn[n] = smallx;
+      eos_state_ash.xn[n] = cj_rp::smallx;
   }
-  eos_state_ash.xn[NumSpec-1] = 1.0_rt - (NumSpec - 1) * smallx;
+  eos_state_ash.xn[NumSpec-1] = 1.0_rt - (NumSpec - 1) * cj_rp::smallx;
 
   // get the q value -- we need the change in molar fractions
   Array1D<Real, 1, NumSpec> dymol;
@@ -83,8 +72,23 @@ int main(int argc, char *argv[]) {
                  eos_state_fuel.xn[n-1] * aion_inv[n-1];
   }
 
-  Real q_burn;
-  ener_gener_rate(dymol, q_burn);
+  Real q_burn{};
+
+#ifdef NEW_NETWORK_IMPLEMENTATION
+
+  // note: we are assuming that the network's ener_gener_rate does not
+  // use rhs_state -- this is true, e.g., for aprox13
+  RHS::rhs_state_t<amrex::Real> state;
+  amrex::constexpr_for<1, NumSpec+1>([&] (auto n)
+  {
+      constexpr int species = n;
+      q_burn += RHS::ener_gener_rate<species>(state, dymol(species));
+  });
+#else
+      ener_gener_rate(dymol, q_burn);
+#endif
+
+  std::cout << "q_burn = " << q_burn << std::endl;
 
   // store the shock adiabat and the detonation adiabat