From f208637ab01507bec5f429f4adab56b81f6e5325 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 5 Jan 2025 18:07:38 -0500 Subject: [PATCH 1/2] fix the CJ detonation utility --- util/cj_detonation/GNUmakefile | 2 +- util/cj_detonation/_parameters | 2 ++ util/cj_detonation/inputs | 3 ++- util/cj_detonation/main.cpp | 15 +++++++-------- util/cj_detonation/probin | 5 ----- 5 files changed, 12 insertions(+), 15 deletions(-) delete mode 100644 util/cj_detonation/probin diff --git a/util/cj_detonation/GNUmakefile b/util/cj_detonation/GNUmakefile index 223308e541..ded3cce357 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 := aprox13 +NETWORK_DIR := subch_simple CONDUCTIVITY_DIR := stellar diff --git a/util/cj_detonation/_parameters b/util/cj_detonation/_parameters index b1191dab7f..d5d3e3ed6d 100644 --- a/util/cj_detonation/_parameters +++ b/util/cj_detonation/_parameters @@ -1,3 +1,5 @@ +@namespace: cj + smallx real 1.e-10 small_temp real 1.e5 diff --git a/util/cj_detonation/inputs b/util/cj_detonation/inputs index 9204ae72be..b0edb3d397 100644 --- a/util/cj_detonation/inputs +++ b/util/cj_detonation/inputs @@ -1 +1,2 @@ -amr.probin_file = probin +cj.smallx = 1.e-10 + diff --git a/util/cj_detonation/main.cpp b/util/cj_detonation/main.cpp index d4b520f1b8..632be518dc 100644 --- a/util/cj_detonation/main.cpp +++ b/util/cj_detonation/main.cpp @@ -11,9 +11,8 @@ using namespace amrex; #include #include #include -#ifndef NEW_NETWORK_IMPLEMENTATION +#include #include -#endif using namespace unit_test_rp; @@ -37,11 +36,11 @@ int main(int argc, char *argv[]) { 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 +59,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 +70,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 dymol; diff --git a/util/cj_detonation/probin b/util/cj_detonation/probin deleted file mode 100644 index f4576dcbd8..0000000000 --- a/util/cj_detonation/probin +++ /dev/null @@ -1,5 +0,0 @@ -&probin - - smallx = 1.e-10 - -/ From 644d6adc4896b6e55d96204a39f675842100b890 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 5 Jan 2025 18:56:21 -0500 Subject: [PATCH 2/2] 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 #include #include +#ifdef NEW_NETWORK_IMPLEMENTATION +#include +#else #include +#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 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 state; + amrex::constexpr_for<1, NumSpec+1>([&] (auto n) + { + constexpr int species = n; + q_burn += RHS::ener_gener_rate(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