From 644d6adc4896b6e55d96204a39f675842100b890 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sun, 5 Jan 2025 18:56:21 -0500
Subject: [PATCH] this now works

---
 util/cj_detonation/GNUmakefile |  2 +-
 util/cj_detonation/README.md   | 24 +++++++++++++++++-----
 util/cj_detonation/cj_det.H    |  2 +-
 util/cj_detonation/main.cpp    | 37 +++++++++++++++++++---------------
 4 files changed, 42 insertions(+), 23 deletions(-)

diff --git a/util/cj_detonation/GNUmakefile b/util/cj_detonation/GNUmakefile
index ded3cce357..223308e541 100644
--- a/util/cj_detonation/GNUmakefile
+++ b/util/cj_detonation/GNUmakefile
@@ -21,7 +21,7 @@ MICROPHYSICS_HOME  := ../..
 EOS_DIR     := helmholtz
 
 # This sets the network directory
-NETWORK_DIR := subch_simple
+NETWORK_DIR := aprox13
 
 CONDUCTIVITY_DIR := stellar
 
diff --git a/util/cj_detonation/README.md b/util/cj_detonation/README.md
index aa3d40b0c6..862ee418bb 100644
--- a/util/cj_detonation/README.md
+++ b/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.
+
diff --git a/util/cj_detonation/cj_det.H b/util/cj_detonation/cj_det.H
index fcba430bb7..4653794dd2 100644
--- a/util/cj_detonation/cj_det.H
+++ b/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;
diff --git a/util/cj_detonation/main.cpp b/util/cj_detonation/main.cpp
index 632be518dc..34add27baa 100644
--- a/util/cj_detonation/main.cpp
+++ b/util/cj_detonation/main.cpp
@@ -12,7 +12,11 @@ using namespace amrex;
 #include <cj_det.H>
 #include <unit_test.H>
 #include <actual_network.H>
+#ifdef NEW_NETWORK_IMPLEMENTATION
+#include <rhs.H>
+#else
 #include <actual_rhs.H>
+#endif
 
 using namespace unit_test_rp;
 
@@ -22,20 +26,6 @@ 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();
 
   // C++ EOS initialization (must be done after Fortran eos_init and
@@ -82,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