From eb5189fa546da7b644850bb11d5639af0eb5f20d Mon Sep 17 00:00:00 2001 From: Max Katz Date: Sun, 23 Jul 2023 18:36:06 -0400 Subject: [PATCH] Convert ignition_chamulak to new reactions scheme (#1296) --- networks/ignition_chamulak/Make.package | 3 - networks/ignition_chamulak/actual_network.H | 131 +++++++++- .../ignition_chamulak/actual_network_data.cpp | 6 - networks/ignition_chamulak/actual_rhs.H | 226 ------------------ .../ignition_chamulak/actual_rhs_data.cpp | 8 - .../ci-benchmarks/chamulak_VODE_unit_test.out | 20 +- unit_test/test_rhs/ci-benchmarks/chamulak.out | 4 +- 7 files changed, 134 insertions(+), 264 deletions(-) delete mode 100644 networks/ignition_chamulak/actual_network_data.cpp delete mode 100644 networks/ignition_chamulak/actual_rhs_data.cpp diff --git a/networks/ignition_chamulak/Make.package b/networks/ignition_chamulak/Make.package index f2ee53d848..80074414e2 100644 --- a/networks/ignition_chamulak/Make.package +++ b/networks/ignition_chamulak/Make.package @@ -1,10 +1,7 @@ CEXE_headers += network_properties.H ifeq ($(USE_REACT),TRUE) - CEXE_sources += actual_network_data.cpp CEXE_headers += actual_network.H - - CEXE_sources += actual_rhs_data.cpp CEXE_headers += actual_rhs.H USE_SCREENING = TRUE diff --git a/networks/ignition_chamulak/actual_network.H b/networks/ignition_chamulak/actual_network.H index f8f339eb1a..344917ab6b 100644 --- a/networks/ignition_chamulak/actual_network.H +++ b/networks/ignition_chamulak/actual_network.H @@ -1,31 +1,144 @@ #ifndef actual_network_H #define actual_network_H +#define NEW_NETWORK_IMPLEMENTATION + #include #include #include #include #include +#include using namespace amrex; -void actual_network_init(); +AMREX_INLINE +void actual_network_init() {} const std::string network_name = "ignition_chamulak"; -namespace network -{ - // M12_chamulak is the effective number of C12 nuclei destroyed per - // reaction - const Real M12_chamulak = 2.93e0_rt; -} - namespace Rates { enum NetworkRates { - NumRates = 1 + chamulak_rate = 1, + NumRates = chamulak_rate }; }; +namespace RHS { + + AMREX_GPU_HOST_DEVICE AMREX_INLINE + constexpr rhs_t rhs_data (int rate) + { + using namespace Species; + using namespace Rates; + + rhs_t data; + + switch (rate) { + + case chamulak_rate: + data.species_A = C12; + data.species_D = Ash; + + // The composite "ash" particle we are creating has a A = 18, + // while we are consuming 2 C12 nucleons to create it. In order + // to conserve mass, we need to scale the number of ash particles + // created by 24 / 18 = 4 / 3. + + data.number_A = 2.0_rt; + data.number_D = 4.0_rt / 3.0_rt; + + data.exponent_A = 2; + data.exponent_D = 1; + + // For each reaction, we lose 2.93 C12 nuclei (for a single rate, + // C12+C12, we would lose 2. In Chamulak et al. (2008), they say a + // value of 2.93 captures the energetics from a larger network. + + data.forward_branching_ratio = 2.93e0_rt / 2.0_rt; + break; + + } + + return data; + } + + template + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void evaluate_analytical_rate (const rhs_state_t& state, rate_t& rates) + { + using namespace Species; + using namespace Rates; + + if constexpr (rate == chamulak_rate) { + // compute some often used temperature constants + Real T9 = state.tf.t9; + Real dT9dt = 1.0_rt / 1.0e9_rt; + Real T9a = T9 / (1.0_rt + 0.0396e0_rt * T9); + Real dT9adt = (T9a / T9 - (T9a / (1.0_rt + 0.0396e0_rt * T9)) * 0.0396e0_rt) * dT9dt; + + // compute the CF88 rate + Real scratch = std::pow(T9a, 1.0_rt / 3.0_rt); + Real dscratchdt = (1.0_rt / 3.0_rt) * std::pow(T9a, -2.0_rt / 3.0_rt) * dT9adt; + + Real a = 4.27e26_rt * std::pow(T9a, 5.0_rt / 6.0_rt) * std::pow(T9, -1.5e0_rt); + Real dadt = (5.0_rt / 6.0_rt) * (a / T9a) * dT9adt - 1.5e0_rt * (a / T9) * dT9dt; + + Real b = std::exp(-84.165e0_rt / scratch - 2.12e-3_rt * T9 * T9 * T9); + Real dbdt = (84.165e0_rt * dscratchdt / (scratch * scratch) - 3.0_rt * 2.12e-3_rt * T9 * T9 * dT9dt) * b; + + rates.fr = a * b; + rates.frdt = dadt * b + a * dbdt; + + rates.rr = 0.0_rt; + rates.rrdt = 0.0_rt; + } + } + + template + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void postprocess_rate (const rhs_state_t& state, rate_t& rates, + rate_t& rates1, rate_t& rates2, rate_t& rates3) + { + // Nothing to do for this network. + } + + template + AMREX_GPU_HOST_DEVICE AMREX_INLINE + Real ener_gener_rate (const rhs_state_t& rhs_state, Real const& dydt) + { + // Chamulak et al. provide the q-value resulting from C12 burning, + // given as 3 different values (corresponding to 3 different densities). + // Here we do a simple quadratic fit to the 3 values provided (see + // Chamulak et al., p. 164, column 2). + + // our convention is that the binding energies are negative. We convert + // from the MeV values that are traditionally written in astrophysics + // papers by multiplying by 1.e6 eV/MeV * 1.60217646e-12 erg/eV. The + // MeV values are per nucleus, so we divide by aion to make it per + // nucleon and we multiple by Avogardo's # (6.0221415e23) to get the + // value in erg/g + Real rho9 = rhs_state.rho / 1.0e9_rt; + + // q_eff is effective heat evolved per reaction (given in MeV) + Real q_eff = 0.06e0_rt * rho9 * rho9 + 0.02e0_rt * rho9 + 8.83e0_rt; + + // convert from MeV to ergs / gram. Here 2.93 is the + // number of C12 nuclei destroyed in a single reaction and 12.0 is + // the mass of a single C12 nuclei. Also note that our convention + // is that binding energies are negative. + Real ebin_C12 = -q_eff * C::MeV2eV * C::ev2erg * C::n_A / (2.93e0_rt * 12.0e0_rt); + + if constexpr (spec == Species::C12) { + return dydt * NetworkProperties::aion(spec) * ebin_C12; + } + else { + return 0.0_rt; + } + } + +} // namespace RHS + #endif diff --git a/networks/ignition_chamulak/actual_network_data.cpp b/networks/ignition_chamulak/actual_network_data.cpp deleted file mode 100644 index 459fb325bd..0000000000 --- a/networks/ignition_chamulak/actual_network_data.cpp +++ /dev/null @@ -1,6 +0,0 @@ -#include -#include - -void actual_network_init() -{ -} diff --git a/networks/ignition_chamulak/actual_rhs.H b/networks/ignition_chamulak/actual_rhs.H index 75cc31ec00..2cf0ad127e 100644 --- a/networks/ignition_chamulak/actual_rhs.H +++ b/networks/ignition_chamulak/actual_rhs.H @@ -1,232 +1,6 @@ #ifndef actual_rhs_H #define actual_rhs_H -#include -#include -#include -#include -#include - -#include -#include #include -#include -#include -#include -#include - -using namespace amrex; -using namespace ArrayUtil; -using namespace Rates; -using namespace network_rp; - -void actual_rhs_init (); - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -Real get_ebin (const Real& density) -{ - // Chamulak et al. provide the q-value resulting from C12 burning, - // given as 3 different values (corresponding to 3 different densities). - // Here we do a simple quadratic fit to the 3 values provided (see - // Chamulak et al., p. 164, column 2). - - // our convention is that the binding energies are negative. We convert - // from the MeV values that are traditionally written in astrophysics - // papers by multiplying by 1.e6 eV/MeV * 1.60217646e-12 erg/eV. The - // MeV values are per nucleus, so we divide by aion to make it per - // nucleon and we multiple by Avogardo's # (6.0221415e23) to get the - // value in erg/g - Real rho9 = density / 1.0e9_rt; - - // q_eff is effective heat evolved per reaction (given in MeV) - Real q_eff = 0.06e0_rt * rho9 * rho9 + 0.02e0_rt * rho9 + 8.83e0_rt; - - // convert from MeV to ergs / gram. Here M12_chamulak is the - // number of C12 nuclei destroyed in a single reaction and 12.0 is - // the mass of a single C12 nuclei. Also note that our convention - // is that binding energies are negative. - return -q_eff * C::MeV2eV * C::ev2erg * C::n_A / (network::M12_chamulak * 12.0e0_rt); -} - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void get_rates (const burn_t& state, Real& rate, Real& dratedt, Real& sc1212, Real& dsc1212dt) -{ - using namespace Species; - - Real temp = state.T; - Real dens = state.rho; - - // convert mass fractions to molar fractions - Array1D y; - for (int i = 1; i <= NumSpec; ++i) { - y(i) = state.xn[i-1] * aion_inv[i-1]; - } - - // call the screening routine - plasma_state_t pstate; - fill_plasma_state(pstate, temp, dens, y); - - int jscr = 0; - - Real dsc1212dd; - screen(pstate, jscr, - zion[C12-1], aion[C12-1], zion[C12-1], aion[C12-1], - sc1212, dsc1212dt, dsc1212dd); - - // compute some often used temperature constants - Real T9 = temp / 1.e9_rt; - Real dT9dt = 1.0_rt / 1.e9_rt; - Real T9a = T9 / (1.0_rt + 0.0396e0_rt * T9); - Real dT9adt = (T9a / T9 - (T9a / (1.0_rt + 0.0396e0_rt * T9)) * 0.0396e0_rt) * dT9dt; - - // compute the CF88 rate - Real scratch = std::pow(T9a, 1.0_rt / 3.0_rt); - Real dscratchdt = (1.0_rt / 3.0_rt) * std::pow(T9a, -2.0_rt / 3.0_rt) * dT9adt; - - Real a = 4.27e26_rt * std::pow(T9a, 5.0_rt / 6.0_rt) * std::pow(T9, -1.5e0_rt); - Real dadt = (5.0_rt / 6.0_rt) * (a / T9a) * dT9adt - 1.5e0_rt * (a / T9) * dT9dt; - - Real b = std::exp(-84.165e0_rt / scratch - 2.12e-3_rt * T9 * T9 * T9); - Real dbdt = (84.165e0_rt * dscratchdt / (scratch * scratch) - 3.0_rt * 2.12e-3_rt * T9 * T9 * dT9dt) * b; - - rate = a * b; - dratedt = dadt * b + a * dbdt; -} - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -Real ener_gener_rate (const Real& dydt_C12, const Real& ebin_C12) -{ - using namespace Species; - - return dydt_C12 * aion[C12-1] * ebin_C12; -} - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void actual_rhs (const burn_t& state, Array1D& ydot) -{ - using namespace Species; - - for (int i = 1; i <= NumSpec; ++i) { - ydot(i) = 0.0_rt; - } - - Real rate, dratedt, sc1212, dsc1212dt; - get_rates(state, rate, dratedt, sc1212, dsc1212dt); - - Real dens = state.rho; - - // we come in with mass fractions -- convert to molar fractions - Array1D y; - for (int i = 1; i <= NumSpec; ++i) { - y(i) = state.xn[i-1] * aion_inv[i-1]; - } - - - // The change in number density of C12 is - // d(n12)/dt = - M12_chamulak * 1/2 (n12)**2 - // - // where is the average of the relative velocity times the - // cross section for the reaction, and the factor accounting for the - // total number of particle pairs has a 1/2 because we are - // considering a reaction involving identical particles (see Clayton - // p. 293). Finally, the -M12_chamulak means that for each reaction, - // we lose M12_chamulak C12 nuclei (for a single rate, C12+C12, - // M12_chamulak would be 2. In Chamulak et al. (2008), they say a - // value of 2.93 captures the energetics from a larger network - // - // Switching over to mass fractions, using n = rho X N_A/A, where N_A is - // Avogadro's number, and A is the mass number of the nucleon, we get - // - // d(X12)/dt = -M12_chamulak * 1/2 (X12)**2 rho N_A / A12 - // - // The quantity [N_A ] is what is tabulated in Caughlin and Fowler. - - Real xc12tmp = amrex::max(y(C12) * aion[C12-1], 0.e0_rt); - ydot(C12) = -(1.0_rt / 12.0_rt) * 0.5_rt * network::M12_chamulak * dens * sc1212 * rate * xc12tmp * xc12tmp; - ydot(O16) = 0.0_rt; - ydot(Ash) = -ydot(C12); - - // Convert back to molar form - - for (int i = 1; i <= NumSpec; ++i) { - ydot(i) *= aion_inv[i-1]; - } - - Real ebin_C12 = get_ebin(dens); - - ydot(net_ienuc) = ener_gener_rate(ydot(C12), ebin_C12); -} - -template -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void actual_jac (burn_t& state, MatrixType& jac) -{ - using namespace Species; - - Real rate, dratedt, sc1212, dsc1212dt; - get_rates(state, rate, dratedt, sc1212, dsc1212dt); - - // initialize - jac.zero(); - - Real temp = state.T; - Real dens = state.rho; - Real xc12tmp = amrex::max(state.xn[C12-1], 0.0_rt); - - // Carbon jacobian elements - - jac(C12, C12) = -2.0_rt * (1.0_rt / 12.0_rt) * network::M12_chamulak * 0.5_rt * dens * sc1212 * rate * xc12tmp; - jac(Ash, C12) = -jac(C12, C12); - - // Add the temperature derivative df(y_i) / dT (we will convert to d / de later) - - jac(C12, net_ienuc) = -(1.0_rt / 12.0_rt) * network::M12_chamulak * 0.5_rt * - (dens * rate * xc12tmp * xc12tmp * dsc1212dt + - dens * sc1212 * xc12tmp * xc12tmp * dratedt); - - jac(Ash, net_ienuc) = -jac(C12, net_ienuc); - - // Convert back to molar form - - for (int i = 1; i <= neqs; ++i) { - for (int j = 1; j <= NumSpec; ++j) { - jac(i,j) *= aion[j-1]; - jac(j,i) *= aion_inv[j-1]; - } - } - - // Energy generation rate Jacobian elements with respect to species - - Real ebin_C12 = get_ebin(dens); - - for (int j = 1; j <= NumSpec; ++j) { - jac(net_ienuc, j) = ener_gener_rate(jac(C12, j), ebin_C12); - } - - // df(e) / dT - - jac(net_ienuc, net_ienuc) = ener_gener_rate(jac(C12, net_ienuc), ebin_C12); - - // Convert to d / de - - jac(C12, net_ienuc) = temperature_to_energy_jacobian(state, jac(C12, net_ienuc)); - jac(Ash, net_ienuc) = temperature_to_energy_jacobian(state, jac(Ash, net_ienuc)); - jac(net_ienuc, net_ienuc) = temperature_to_energy_jacobian(state, jac(net_ienuc, net_ienuc)); -} - -// Compute and store the more expensive screening factors - -AMREX_INLINE -void set_up_screening_factors () -{ - using namespace Species; - - // note: it is critical that these are called in the exact order - // that the screening calls are done in the RHS routine, since we - // use that order in the screening - - int jscr = 0; - add_screening_factor(jscr++, zion[C12-1], aion[C12-1], zion[C12-1], aion[C12-1]); -} #endif diff --git a/networks/ignition_chamulak/actual_rhs_data.cpp b/networks/ignition_chamulak/actual_rhs_data.cpp deleted file mode 100644 index 01de318328..0000000000 --- a/networks/ignition_chamulak/actual_rhs_data.cpp +++ /dev/null @@ -1,8 +0,0 @@ -#include - -void actual_rhs_init() -{ - set_up_screening_factors(); - - screening_init(); -} diff --git a/unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out b/unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out index ecd2d7288e..cc751a0856 100644 --- a/unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/chamulak_VODE_unit_test.out @@ -12,21 +12,21 @@ RHS at t = 0 ash 0.01230280576 ------------------------------------ successful? 1 - - Hnuc = 5.27728638e+17 - - added e = 8.364498912e+15 - - final T = 1433703794 + - Hnuc = 5.274803093e+17 + - added e = 8.360562902e+15 + - final T = 1433512646 ------------------------------------ e initial = 1.253426044e+18 -e final = 1.261790543e+18 +e final = 1.261786607e+18 ------------------------------------ new mass fractions: -C12 0.9657902582 -O16 1e-30 -ash 0.03420974184 +C12 0.9658063217 +O16 9.999999766e-31 +ash 0.03419367826 ------------------------------------ species creation rates: -omegadot(C12): -2.158343334 -omegadot(O16): 3.315374916e-44 -omegadot(ash): 2.158343334 +omegadot(C12): -2.157329859 +omegadot(O16): -1.476105976e-36 +omegadot(ash): 2.157329859 number of steps taken: 381 AMReX (23.07-7-g88f03408f18a-dirty) finalized diff --git a/unit_test/test_rhs/ci-benchmarks/chamulak.out b/unit_test/test_rhs/ci-benchmarks/chamulak.out index f05f26c3a2..946d240f01 100644 --- a/unit_test/test_rhs/ci-benchmarks/chamulak.out +++ b/unit_test/test_rhs/ci-benchmarks/chamulak.out @@ -17,11 +17,11 @@ J_carbon-12_oxygen-16 0 0 J_oxygen-16_oxygen-16 0 0 J_ash_oxygen-16 0 0 - J_E_oxygen-16 -0 -0 + J_E_oxygen-16 0 0 J_carbon-12_ash 0 0 J_oxygen-16_ash 0 0 J_ash_ash 0 0 - J_E_ash -0 -0 + J_E_ash 0 0 J_carbon-12_E -4.4154660816e-08 -7.3754668622e-85 J_oxygen-16_E 0 0 J_ash_E 4.9169779082e-85 2.9436440544e-08