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/_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/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/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..34add27baa 100644 --- a/util/cj_detonation/main.cpp +++ b/util/cj_detonation/main.cpp @@ -11,7 +11,10 @@ using namespace amrex; #include #include #include -#ifndef NEW_NETWORK_IMPLEMENTATION +#include +#ifdef NEW_NETWORK_IMPLEMENTATION +#include +#else #include #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 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 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 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 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 - -/