From 9e2281237cac5e059f75d77cd9b1cbd0077d8176 Mon Sep 17 00:00:00 2001 From: zhi Date: Sun, 4 Jan 2026 21:06:04 -0500 Subject: [PATCH 1/5] this compiles --- networks/general_null/xrb_ignition.net | 8 + util/xrb_ignition/GNUmakefile | 40 ++ util/xrb_ignition/Make.package | 2 + util/xrb_ignition/README.md | 14 + util/xrb_ignition/_parameters | 9 + util/xrb_ignition/ignition_condition.H | 515 +++++++++++++++++++++++++ util/xrb_ignition/inputs_ignition | 5 + util/xrb_ignition/main.cpp | 33 ++ 8 files changed, 626 insertions(+) create mode 100644 networks/general_null/xrb_ignition.net create mode 100644 util/xrb_ignition/GNUmakefile create mode 100644 util/xrb_ignition/Make.package create mode 100644 util/xrb_ignition/README.md create mode 100644 util/xrb_ignition/_parameters create mode 100644 util/xrb_ignition/ignition_condition.H create mode 100644 util/xrb_ignition/inputs_ignition create mode 100644 util/xrb_ignition/main.cpp diff --git a/networks/general_null/xrb_ignition.net b/networks/general_null/xrb_ignition.net new file mode 100644 index 000000000..880cd9081 --- /dev/null +++ b/networks/general_null/xrb_ignition.net @@ -0,0 +1,8 @@ +# this is the null description for the X-ray burst ignition condition +# + +# name short-name aion zion + hydrogen-1 H1 1.0 1.0 + helium-4 He4 4.0 2.0 + oxygen-14 O14 14.0 8.0 + oxygen-15 O15 15.0 8.0 diff --git a/util/xrb_ignition/GNUmakefile b/util/xrb_ignition/GNUmakefile new file mode 100644 index 000000000..7b6f6b451 --- /dev/null +++ b/util/xrb_ignition/GNUmakefile @@ -0,0 +1,40 @@ +PRECISION = DOUBLE +PROFILE = FALSE + +DEBUG = FALSE + +DIM = 3 + +COMP = gnu + +USE_MPI = FALSE +USE_OMP = FALSE + +USE_REACT = TRUE + +EBASE = main + +# define the location of the Microphysics top directory +MICROPHYSICS_HOME := ../.. + +# This sets the EOS directory +EOS_DIR := helmholtz + +# This sets the network directory +NETWORK_DIR := general_null +NETWORK_INPUTS := xrb_ignition.net + +USE_CONDUCTIVITY := TRUE +CONDUCTIVITY_DIR := stellar + +USE_SCREENING := TRUE +SCREEN_METHOD := screen5 + +USE_SIMPLIFIED_SDC = TRUE + +EXTERN_SEARCH += . + +Bpack := ./Make.package +Blocs := . + +include $(MICROPHYSICS_HOME)/unit_test/Make.unit_test diff --git a/util/xrb_ignition/Make.package b/util/xrb_ignition/Make.package new file mode 100644 index 000000000..68346547e --- /dev/null +++ b/util/xrb_ignition/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += main.cpp +CEXE_headers += ignition_condition.H diff --git a/util/xrb_ignition/README.md b/util/xrb_ignition/README.md new file mode 100644 index 000000000..7dbb9c9d1 --- /dev/null +++ b/util/xrb_ignition/README.md @@ -0,0 +1,14 @@ +# `xrb_ignition` + +`xrb_ignition` finds the ignition condition of the X-ray burst +based on https://github.com/andrewcumming/settle +See DOI: 10.1086/317191 for details. + +Usage: +./main.ex + +X: Hydrogen massfraction +Z: CNO massfraction +Fb: Base Flux in MeV/nucleon/mdot +mdot: Local accretion rate in Eddington unit +COMPRESS: To include compressional heating or not. diff --git a/util/xrb_ignition/_parameters b/util/xrb_ignition/_parameters new file mode 100644 index 000000000..ef1dbc7c2 --- /dev/null +++ b/util/xrb_ignition/_parameters @@ -0,0 +1,9 @@ +@namespace: unit_test + +run_prefix string "" + +X real 0.71 +Z real 0.01 +mdot real 0.1 +Fb real 150.0 +COMPRESS int 0 diff --git a/util/xrb_ignition/ignition_condition.H b/util/xrb_ignition/ignition_condition.H new file mode 100644 index 000000000..f20699d0d --- /dev/null +++ b/util/xrb_ignition/ignition_condition.H @@ -0,0 +1,515 @@ +#ifndef IGNITION_CONDITION_H +#define IGNITION_CONDITION_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* mass fractions of CNO elements (metal) based on beta-decay timescale + O14 -> N14 (70.641 s) and O15 -> N15 (122.24 s).*/ + +constexpr amrex::Real O14 = 0.352; +constexpr amrex::Real O15 = 0.648; + +struct State +{ + // Radius and gravitational acceleration + amrex::Real R, g; + + // Initial abundance for hydrogen, helium, and metal + amrex::Real X0, Y0, Z0; + + // Local Accretion rate + amrex::Real mdot; + + // Depletion Column Depth + amrex::Real yd; + + // Top boundary column depth, temperature and flux + // This is used as the initial condition from integration. + amrex::Real yt, Tt, Ft; + + // Bot boundary column depth, temperature and flux + // This is assumed to be the ignition condition + amrex::Real yb, Tb, Fb; +}; + + +template +void rk4(amrex::Real t0, amrex::Real tmax, amrex::Real stepsize, + const amrex::Real (&f0)[NEQ], const State& state, + void (*rhs)(const amrex::Real t, const amrex::Real (&f)[NEQ], + const State& state, amrex::Real (&dfdt)[NEQ]), + amrex::Real* t_array, amrex::Real (*f_array)[NEQ], int nsteps) { + /* Runge-Kutta 4 IVP solver for multiple variables: df1/dt, df2dt, ... + where NEQ is the number of variables. + + Parameters + ---------- + t0: Initial time + tmax: Final time + stepsize: Stepsize used for integration + f0: Initial value corresponding to t0, f(t0) = f0. + rhs: Function that computes the derivative, df / dt + */ + + // Set initial condition + amrex::Real dt = stepsize; + + t_array[0] = t0; + amrex::Real t = t0; + int i =0; + + for (int n = 0; n < NEQ; ++n) { + f_array[0][n] = f0[n]; + } + + amrex::Real k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ]; + amrex::Real ytmp[NEQ]; + + while (t < tmax) { + + // Change stepsize so that we don't overshoot + if (t + stepsize > tmax) { + dt = tmax - t; + } + + // Array of y's of size NEQ, number of variables. + const amrex::Real (&y)[NEQ] = f_array[i]; + + rhs(t, y, state, k1); + + for (int n = 0; n < NEQ; ++n) { + ytmp[n] = y[n] + 0.5_rt * dt * k1[n]; + } + rhs(t + 0.5_rt * dt, ytmp, state, k2); + + for (int n = 0; n < NEQ; ++n) { + ytmp[n] = y[n] + 0.5_rt * dt * k2[n]; + } + rhs(t + 0.5_rt * dt, ytmp, state, k3); + + for (int n = 0; n < NEQ; ++n) { + ytmp[n] = y[n] + dt * k3[n]; + } + rhs(t + dt, ytmp, state, k4); + + for (int n = 0; n < NEQ; ++n) { + f_array[i+1][n] = y[n] + (dt / 6.0_rt) * + (k1[n] + 2.0_rt*k2[n] + 2.0_rt*k3[n] + k4[n]); + } + + t += dt; + t_array[i+1] = t; + i++; + } +} + + +amrex::Real zbrent(amrex::Real a, + amrex::Real b, + amrex::Real tol, + int max_iter, + amrex::Real (*fcn)(amrex::Real, State&), + State& state) { + /* + Brent's method for finding root without Jacobian within interval [a, b] + */ + + // Compute the constraint equation at the boundary + amrex::Real fa = fcn(a, state); + amrex::Real fb = fcn(b, state); + + // Makre sure roots are bracketed + // i.e. they have opposite sign + if (fa * fb > 0.0_rt) { + amrex::Abort("Brent: root not bracketed"); + } + + // Make large value on the left? + if (std::abs(fa) < std::abs(fb)) { + std::swap(a, b); + std::swap(fa, fb); + } + + amrex::Real c = a; + amrex::Real fc = fa; + + // Flag to see if we went through bisection method previously + bool mflag = true; + + // s is the new solution + amrex::Real s = b; + amrex::Real d = 0.0_rt; + + for (int iter = 0; iter < max_iter; ++iter) { + + if (fa != fc && fb != fc) { + // Inverse quadratic interpolation if we have 3 distinct points + // a=x_n-2, b=x_n-1, c=x_n + s = a * fb * fc / ((fa - fb) * (fa - fc)) + + b * fa * fc / ((fb - fa) * (fb - fc)) + + c * fa * fb / ((fc - fa) * (fc - fb)); + } else { + // Secant method otherwise + s = b - fb * (b - a) / (fb - fa); + } + + // Conditions to fall back to bisection + if ((s < (3*a + b) * 0.25_rt || s > b) || + (mflag && std::abs(s - b) >= 0.5_rt * std::abs(b - c)) || + (!mflag && std::abs(s - b) >= 0.5_rt * std::abs(c - d)) || + (mflag && std::abs(b - c) < tol) || + (!mflag && std::abs(c - d) < tol)) { + + s = 0.5_rt * (a + b); + mflag = true; + } else { + mflag = false; + } + + amrex::Real fs = fcn(s, state); + + // d is previous value of c + d = c; + + // c is previous value of b + c = b; + fc = fb; + + if (fa * fs < 0.0_rt) { + b = s; + fb = fs; + } else { + a = s; + fa = fs; + } + + if (std::abs(fa) < std::abs(fb)) { + std::swap(a, b); + std::swap(fa, fb); + } + + // Convergence test + if (std::abs(b - a) < tol || std::abs(fb) < tol) { + return b; + } + } + + amrex::Abort("Brent: did not converge"); + return b; // unreachable +} + + +void derivs(const amrex::Real y, const amrex::Real (&f)[3], + const State& state, amrex::Real (&dfdy)[3]) { + /* + derivatives for integration: + dfdy = (dT/dy, dF/dy, dz/dy) and f = (T, F, z) + + dT / dy = 3 κ / (4 a c) F / T^3 + dT / dy = F / (ϱ k_th) + where k_th = (4 a c T^3) / (3 κ ϱ) + + dF/dy = -epsilon_H if X > 0, X = X0(1 - y/yd) + + dz/dy = -1 / ϱ + */ + eos_t eos_state; + eos_state.T = std::max(f[0], 1.0e7); + eos_state.p = state.g * y; + + amrex::Real X = std::max(state.X0 * (1.0_rt - y / state.yd), 0.0_rt); + eos_state.xn[Species::H1-1] = X; + eos_state.xn[Species::He4-1] = 1.0_rt - X - state.Z0; + eos_state.xn[Species::O14-1] = O14 * state.Z0; + eos_state.xn[Species::O15-1] = O15 * state.Z0; + + // Find density via Temperature and pressure + eos(eos_input_tp, eos_state); + + // Find conductivity + actual_conductivity(eos_state); + + amrex::Real dTdy = f[1] / (eos_state.rho * eos_state.conductivity); + + dfdy[0] = dTdy; + + // Hot CNO burning + amrex::Real dFdy = 0.0_rt; + amrex::Real epsilon_H = 5.8e13 * (state.Z0 / 0.01); + if (X > 0.0_rt) { + dFdy = -epsilon_H; + } + dfdy[1] = dFdy; + + // Hydrostatic balance + dfdy[2] = - 1.0_rt / eos_state.rho; +} + + +amrex::Real heating_rate(const eos_t& eos_state) { + /* + Evaluates the heating term -- triple alpha process + */ + + // Q-value of C12 + 2p -> O14 ( C12(p,g)N13(p,g)O14 ) = 6.57 MeV + // Q-value 3 He4 -> C12 = 7.274 MeV + constexpr amrex::Real Q12 = 6.57_rt; + constexpr amrex::Real Q3a = 7.274_rt; + constexpr amrex::Real enhance_fac = 1.0 + Q12 / Q3a; + + // Find screening factor + constexpr int do_T_derivatives = 1; + using number_t = std::conditional_t; + number_t temp = eos_state.T; + autodiff::seed(temp); + + plasma_state_t pstate{}; + amrex::Array1D Y; + for (int n = 0; n < NumSpec; ++n) { + Y(n+1) = eos_state.xn[n] * aion_inv[n]; + } + + amrex::Real scor, scor2, dscor_dt, dscor2_dt; + + fill_plasma_state(pstate, temp, eos_state.rho, Y); + + constexpr auto scn_fac = scrn::calculate_screen_factor(2.0_rt, 4.0_rt, 2.0_rt, 4.0_rt); + actual_screen(pstate, scn_fac, scor, dscor_dt); + + constexpr auto scn_fac2 = scrn::calculate_screen_factor(2.0_rt, 4.0_rt, 4.0_rt, 8.0_rt); + actual_screen(pstate, scn_fac2, scor2, dscor2_dt); + + amrex::Real Fscn = scor * scor2; + amrex::Real dFscn_dT = scor * dscor2_dt + dscor_dt * dscor2_dt; + + // Maybe use actual triple alpha from pynucastro? But then + // eps_3a = 5.3e21 F_scn rho5^2 Y^3 / T8^3 exp(-44 / T8) + amrex::Real epsilon_3a = 5.3e21 * std::pow(eos_state.rho * 1.0e-5, 2) * + std::pow(eos_state.xn[Species::He4-1], 3) * std::pow(eos_state.T * 1.0e-8, -3) * + std::exp(-44.0 * 1.0e8 / eos_state.T); + + // If there is enough hydrogen for the rapid proton capture reactions: + // C12(p,g)N13(p,g)O14, i.e. when X(H1) > 1.0e-3, + // then apply an additional enhancement factor due to enhancement, ~ 1.9 + if (eos_state.xn[Species::H1-1] > 1.0e-3_rt) { + epsilon_3a *= enhance_fac; + } + + // Heating term is temperature derivative of epsilon_3a including screening enhancement. + amrex::Real heat = (-3.0_rt * Fscn / eos_state.T + + 44.0_rt * Fscn / (eos_state.T * eos_state.T * 1.0e-8_rt) + + dFscn_dT) * epsilon_3a; + + return heat; +} + +amrex::Real cooling_rate(const amrex::Real yb, const eos_t& eos_state) { + /* + Evaluates the cooling term. + Get cooling rate by combining: + + dT/dy = 3 κ / (4 a c) F / T^3 and ε = -dF/dy + + ε ≅ a c / (3 κ) 1 / (y^2) T^4 + + with k_th = (4 a c T^3) / (3 κ ϱ) + + ε ≅ 0.25 k_th ϱ T / y^2 + Note that since opacity is also temperature dependent, + we should do a finite difference here. + */ + + // Find conductivity + auto tmp_state = eos_state; + actual_conductivity(tmp_state); + amrex::Real k_th = tmp_state.conductivity; + + // Find conductivity temperature dependence via finite difference + + tmp_state.T = eos_state.T * 1.001_rt; + actual_conductivity(tmp_state); + amrex::Real k_hi = tmp_state.conductivity; + + tmp_state.T = eos_state.T / 1.001_rt; + actual_conductivity(tmp_state); + amrex::Real k_lo = tmp_state.conductivity; + + amrex::Real d_kth_dT = (k_hi - k_lo) / (eos_state.T * (1.001_rt - 1.0_rt / 1.001_rt)); + + amrex::Real cool = 0.25_rt * eos_state.rho * + (k_th + eos_state.T * d_kth_dT) / (yb * yb); + + return cool; +} + + +amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) { + /* Constraint equation to find ignition column depth + given the guess ignition column depth, yb_0. + Note the guess value is given in log10(yb) */ + + state.yb = std::pow(log10yb_0, 10); + + // Set upper boundary conditions, i.e. initial conditions: yt, Tt, and Ft + + // Pick arbitrary top boundary column depth + // Follow "settle" version on Github + state.yt = 1.e3_rt; + + // Compute the top boundary flux + // Ft = F_b + epsilon_H * (yb or yd) + // whichever is smaller. + // If yb < yd, the ignition happens during hydrogen-helium layer + // If yb > yd, the ignition happens during pure helium layer + + amrex::Real epsilon_H = 5.8e13 * (state.Z0 / 0.01); + state.Ft = state.Fb + epsilon_H * std::min(state.yb, state.yd); + + // Compute the top boundary temperature + // Given by the analytic radiative zero solution for constant flux + // Tt = std::pow(3 κ Ft yt/ [a c], 1/4) + + // Here 2650 is 3 κ / (a c), which is precomputed. + // where they assumed the κ at top is 0.22 + + // If we were to compute kappa manually, we need eos. + // i.e. we need two inputs. We have p=g * yt, so we need rho? + + state.Tt = std::pow(2650.0 * state.Ft * state.yt, 0.25); + + // Now set the lower boundary conditions: yb, Tb, and Fb + // Here Fb is assumed to be given. Note that the original version + // includes a special term for compressional heating, if thats included + // then we have to solve for Fb. However, it is sufficient to include + // an extra constant heating term to account for compressional heating. + // yb is what we're after, so we need Tb. + + // The original version seems to use an ODE solver for this. + // The idea is that opacity is also temperature dependent, + // so we can't integrate dT/dy = 3 κ / (4 a c) F / T^3 easily + // Also note that dT/dy depends on F, so we need to evolve F as well with dF/dy + // dF/dy = -epsilon_H, where + // Use RK4 to integrate from Tt and Ft all the way down to find T(yb) and F(yb) + // How do I make sure F(yb) = Fb? + // Note that we used Fb to find Ft: Ft = epsilon_H + + // Use state.yt as stepsize and find array_size + int nsteps = static_cast(std::ceil((state.yb - state.yt) / state.yt)) + 1; + + constexpr int NEQ = 3; + amrex::Real y_array[nsteps]; + amrex::Real f_array[nsteps][NEQ]; + amrex::Real f0_array[NEQ] = {state.Tt, state.Ft, 0.0_rt}; + + // t0=yt, tmax=yb, step_size=yt, f0=Tt, + rk4(state.yt, state.yb, state.yt, + f0_array, state, derivs, + y_array, f_array, nsteps); + + // Update solution Tb to State + state.Tb = f_array[nsteps-1][0]; + + // Make sure solution to F(yb) is the same as Fb + + assert(state.Fb - f_array[nsteps-1][-1] < 1.0e-3); + + // Since we need density and opacity, get eos first. + + eos_t eos_state; + eos_state.T = state.Tb; + eos_state.p = state.g * state.yb; + + amrex::Real X = std::max(state.X0 * (1.0_rt - state.yb / state.yd), 0.0_rt); + eos_state.xn[Species::H1-1] = X; + eos_state.xn[Species::He4-1] = 1 - X - state.Z0; + eos_state.xn[Species::O14-1] = O14 * state.Z0; + eos_state.xn[Species::O15-1] = O15 * state.Z0; + + // Find density via Temperature and pressure + eos(eos_input_tp, eos_state); + + // After Tb is solved, assume we have runaway when + // heating rate is larger than cooling rate, i.e. + // d eps_3a / dT = d eps_cool / dT + + amrex::Real heat = heating_rate(eos_state); + amrex::Real cool = cooling_rate(state.yb, eos_state); + return (heat - cool) / cool; +} + + +void find_yb(State& state) { + // Root-find using the constraint equation to find yb + // Use zbrent method so that Jacobian is not required. + // Assumes that yb is between 10^8 ~ 10^10 + + zbrent(8.0_rt, 10.0_rt, 1.0e-3_rt, 1000, yb_constraint_eq, state); +} + + +AMREX_INLINE +void find_ignition_condition() { + /* + Main driver to determine the ignition condition + */ + + // Gravitational Acceleration for 1.4 solar mass and 11 km NS. + constexpr amrex::Real g = 1.5e14_rt; + constexpr amrex::Real R = 1.1e6_rt; + + State state; + + // Read in input parameters: X, Z, F_b, and mdot + + state.X0 = unit_test_rp::X; + state.Z0 = unit_test_rp::Z; + + // Compute the Eddington Accretion Rate + // mdot_Edd = 2 m_p c / [(1 + X) R σ_th ] + // They used a value of 8.8e4 g cm^-2 s^-1 + // We can compute it at runtime in the future + amrex::Real mdot_Edd = 8.8e4_rt; // Come back later + + // Input mdot should be in eddington unit. + amrex::Real mdot = unit_test_rp::mdot * mdot_Edd; + + // Input Fb is MeV/nucleon/mdot. + // Brown & Bildstein showed that electron capture rates can + // release ~ 1 MeV/nucleon deep in the crust + // The compressional term can give additional c_p T ~ 20 T8 keV/nucleon + // If don't include compressional term, use higher constant Fb + // to account for compressional heating. They use ~1.5 MeV/nucleon + // Convert to CGS. + state.Fb = unit_test_rp::Fb * state.mdot * 14.62_rt * 5.8e21_rt; // Not sure why multiplied by these constants, supposedly it is to convert to CGS... + + // Get the initial helium mass fraction + state.Y0 = 1.0_rt - state.X0 - state.Z0; + + // Energy release per gram from hot CNO. + // This is mass difference between helium and 4 hydrogen + constexpr amrex::Real E_H = 6.4e18_rt; + + // Specific nuclear energy generation rate for hot CNO + amrex::Real epsilon_H = 5.8e13 * (state.Z0 / 0.01); + + // Depletion Column Depth: yd = X mdot E_H / epsilon_H + state.yd = state.X0 * state.mdot * E_H / epsilon_H; + + // Now compute the base column depth: yb via root find. + // This is also the ignition column depth + // Assume that we work with log10(yb) + + find_yb(state); +} + +#endif diff --git a/util/xrb_ignition/inputs_ignition b/util/xrb_ignition/inputs_ignition new file mode 100644 index 000000000..d18036b90 --- /dev/null +++ b/util/xrb_ignition/inputs_ignition @@ -0,0 +1,5 @@ +unit_test.X = 0.71 +unit_test.Z = 0.01 +unit_test.mdot = 0.1 +unit_test.Fb = 150.0 +unit_test.COMPRESS = 0 diff --git a/util/xrb_ignition/main.cpp b/util/xrb_ignition/main.cpp new file mode 100644 index 000000000..23530071e --- /dev/null +++ b/util/xrb_ignition/main.cpp @@ -0,0 +1,33 @@ +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +int main(int argc, char *argv[]) { + + amrex::Initialize(argc, argv); + + std::cout << "starting the single zone burn..." << std::endl; + + amrex::ParmParse ppa("amr"); + + init_unit_test(); + + // C++ EOS initialization (must be done after init_extern_parameters) + eos_init(unit_test_rp::small_temp, unit_test_rp::small_dens); + + // C++ Network, RHS, screening, rates initialization + network_init(); + + find_ignition_condition(); + + amrex::Finalize(); +} From d98812714784bef5839bb110c0ac840b2da13bf3 Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 7 Jan 2026 12:38:35 -0500 Subject: [PATCH 2/5] this runs --- util/xrb_ignition/ignition_condition.H | 118 +++++++++++++++++-------- 1 file changed, 81 insertions(+), 37 deletions(-) diff --git a/util/xrb_ignition/ignition_condition.H b/util/xrb_ignition/ignition_condition.H index f20699d0d..14174cf5f 100644 --- a/util/xrb_ignition/ignition_condition.H +++ b/util/xrb_ignition/ignition_condition.H @@ -40,17 +40,20 @@ struct State // Bot boundary column depth, temperature and flux // This is assumed to be the ignition condition amrex::Real yb, Tb, Fb; + + // Height below the top layer corresponding to yt, + amrex::Real z; }; template -void rk4(amrex::Real t0, amrex::Real tmax, amrex::Real stepsize, +void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real stepsize, const amrex::Real (&f0)[NEQ], const State& state, - void (*rhs)(const amrex::Real t, const amrex::Real (&f)[NEQ], - const State& state, amrex::Real (&dfdt)[NEQ]), - amrex::Real* t_array, amrex::Real (*f_array)[NEQ], int nsteps) { + void (*rhs)(const amrex::Real, const amrex::Real (&)[NEQ], + const State&, amrex::Real (&)[NEQ]), + amrex::Real (&f_new)[NEQ]) { /* Runge-Kutta 4 IVP solver for multiple variables: df1/dt, df2dt, ... - where NEQ is the number of variables. + where NEQ is the number of variables. It just solves to tmax. Parameters ---------- @@ -62,18 +65,16 @@ void rk4(amrex::Real t0, amrex::Real tmax, amrex::Real stepsize, */ // Set initial condition - amrex::Real dt = stepsize; - t_array[0] = t0; + amrex::Real dt = stepsize; amrex::Real t = t0; - int i =0; for (int n = 0; n < NEQ; ++n) { - f_array[0][n] = f0[n]; + f_new[n] = f0[n]; } amrex::Real k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ]; - amrex::Real ytmp[NEQ]; + amrex::Real ftmp[NEQ], f[NEQ]; while (t < tmax) { @@ -82,34 +83,35 @@ void rk4(amrex::Real t0, amrex::Real tmax, amrex::Real stepsize, dt = tmax - t; } - // Array of y's of size NEQ, number of variables. - const amrex::Real (&y)[NEQ] = f_array[i]; + // Set current f + for (int n = 0; n < NEQ; ++n) { + f[n] = f_new[n]; + } - rhs(t, y, state, k1); + rhs(t, f, state, k1); for (int n = 0; n < NEQ; ++n) { - ytmp[n] = y[n] + 0.5_rt * dt * k1[n]; + ftmp[n] = f[n] + 0.5_rt * dt * k1[n]; } - rhs(t + 0.5_rt * dt, ytmp, state, k2); + rhs(t + 0.5_rt * dt, ftmp, state, k2); for (int n = 0; n < NEQ; ++n) { - ytmp[n] = y[n] + 0.5_rt * dt * k2[n]; + ftmp[n] = f[n] + 0.5_rt * dt * k2[n]; } - rhs(t + 0.5_rt * dt, ytmp, state, k3); + rhs(t + 0.5_rt * dt, ftmp, state, k3); for (int n = 0; n < NEQ; ++n) { - ytmp[n] = y[n] + dt * k3[n]; + ftmp[n] = f[n] + dt * k3[n]; } - rhs(t + dt, ytmp, state, k4); + rhs(t + dt, ftmp, state, k4); for (int n = 0; n < NEQ; ++n) { - f_array[i+1][n] = y[n] + (dt / 6.0_rt) * + f_new[n] = f[n] + (dt / 6.0_rt) * (k1[n] + 2.0_rt*k2[n] + 2.0_rt*k3[n] + k4[n]); } + // update time t += dt; - t_array[i+1] = t; - i++; } } @@ -131,7 +133,7 @@ amrex::Real zbrent(amrex::Real a, // Makre sure roots are bracketed // i.e. they have opposite sign if (fa * fb > 0.0_rt) { - amrex::Abort("Brent: root not bracketed"); + std::cout << "Warning: root not bracketed" << std::endl; } // Make large value on the left? @@ -358,7 +360,7 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) { given the guess ignition column depth, yb_0. Note the guess value is given in log10(yb) */ - state.yb = std::pow(log10yb_0, 10); + state.yb = std::pow(10, log10yb_0); // Set upper boundary conditions, i.e. initial conditions: yt, Tt, and Ft @@ -403,25 +405,25 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) { // How do I make sure F(yb) = Fb? // Note that we used Fb to find Ft: Ft = epsilon_H - // Use state.yt as stepsize and find array_size - int nsteps = static_cast(std::ceil((state.yb - state.yt) / state.yt)) + 1; - constexpr int NEQ = 3; - amrex::Real y_array[nsteps]; - amrex::Real f_array[nsteps][NEQ]; + amrex::Real f_array[NEQ]; amrex::Real f0_array[NEQ] = {state.Tt, state.Ft, 0.0_rt}; - // t0=yt, tmax=yb, step_size=yt, f0=Tt, + // t0=yt, tmax=yb, step_size=yt, f0={Tt, Ft, 0.0} + // Integrates to find f_array corresponding to yb. rk4(state.yt, state.yb, state.yt, f0_array, state, derivs, - y_array, f_array, nsteps); + f_array); // Update solution Tb to State - state.Tb = f_array[nsteps-1][0]; + state.Tb = f_array[0]; // Make sure solution to F(yb) is the same as Fb - assert(state.Fb - f_array[nsteps-1][-1] < 1.0e-3); + assert(std::abs(state.Fb - f_array[1]) < 1.0e-3 * state.Fb); + + // update relative height + state.z = f_array[2]; // Since we need density and opacity, get eos first. @@ -444,7 +446,15 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) { amrex::Real heat = heating_rate(eos_state); amrex::Real cool = cooling_rate(state.yb, eos_state); - return (heat - cool) / cool; + amrex::Real diff = (heat - cool) / cool; + + std::cout << std::scientific; + std::cout << "yb = " << state.yb + << ", Tb = " << state.Tb + << ", (heat-cool) / cool = " << diff + << std::endl; + + return diff; } @@ -464,11 +474,15 @@ void find_ignition_condition() { */ // Gravitational Acceleration for 1.4 solar mass and 11 km NS. - constexpr amrex::Real g = 1.5e14_rt; + // constexpr amrex::Real g = 1.5e14_rt; + constexpr amrex::Real g = 2.45e14_rt; // This is what settle version uses. constexpr amrex::Real R = 1.1e6_rt; State state; + state.g = g; + state.R = R; + // Read in input parameters: X, Z, F_b, and mdot state.X0 = unit_test_rp::X; @@ -481,7 +495,7 @@ void find_ignition_condition() { amrex::Real mdot_Edd = 8.8e4_rt; // Come back later // Input mdot should be in eddington unit. - amrex::Real mdot = unit_test_rp::mdot * mdot_Edd; + state.mdot = unit_test_rp::mdot * mdot_Edd; // Input Fb is MeV/nucleon/mdot. // Brown & Bildstein showed that electron capture rates can @@ -490,7 +504,8 @@ void find_ignition_condition() { // If don't include compressional term, use higher constant Fb // to account for compressional heating. They use ~1.5 MeV/nucleon // Convert to CGS. - state.Fb = unit_test_rp::Fb * state.mdot * 14.62_rt * 5.8e21_rt; // Not sure why multiplied by these constants, supposedly it is to convert to CGS... + // Here they used mdot in Eddington unit + state.Fb = unit_test_rp::Fb * unit_test_rp::mdot * 14.62_rt * 5.8e21_rt; // Not sure why multiplied by these constants, supposedly it is to convert to CGS... // Get the initial helium mass fraction state.Y0 = 1.0_rt - state.X0 - state.Z0; @@ -505,11 +520,40 @@ void find_ignition_condition() { // Depletion Column Depth: yd = X mdot E_H / epsilon_H state.yd = state.X0 * state.mdot * E_H / epsilon_H; + // Display input parameters to stdout + std::cout << "Base flux is " << unit_test_rp::Fb << std::endl; + std::cout << "In CGS, base flux is " << state.Fb << std::endl; + std::cout << "Metallicity is " << unit_test_rp::Z << std::endl; + std::cout << "Accreted H fraction is " << unit_test_rp::X << std::endl; + std::cout << "mdot/Edd is " << unit_test_rp::mdot << std::endl; + std::cout << "Depletion column depth is " << state.yd << std::endl; + // Now compute the base column depth: yb via root find. // This is also the ignition column depth // Assume that we work with log10(yb) + std::cout << "Searching for ignition column depth..." << std::endl; find_yb(state); + + // Use solution to find density via EOS + + eos_t eos_state; + eos_state.T = state.Tb; + eos_state.p = state.g * state.yb; + + amrex::Real X = std::max(state.X0 * (1.0_rt - state.yb / state.yd), 0.0_rt); + eos_state.xn[Species::H1-1] = X; + eos_state.xn[Species::He4-1] = 1 - X - state.Z0; + eos_state.xn[Species::O14-1] = O14 * state.Z0; + eos_state.xn[Species::O15-1] = O15 * state.Z0; + + // Find density via Temperature and pressure + eos(eos_input_tp, eos_state); + + std::cout << "Solved ignition column depth: y = " << state.yb << std::endl; + std::cout << "With ignition temperature: T = " << state.Tb << std::endl; + std::cout << "With ignition pressure: P = " << state.g * state.yb << std::endl; + std::cout << "With ignition density: ϱ = " << eos_state.rho << std::endl; } #endif From 2ef9ad94d1fcfadfe73cfca4a0addc1bde6496ce Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 7 Jan 2026 12:59:59 -0500 Subject: [PATCH 3/5] add adaptive rk --- util/xrb_ignition/ignition_condition.H | 139 +++++++++++++++++++++++-- 1 file changed, 131 insertions(+), 8 deletions(-) diff --git a/util/xrb_ignition/ignition_condition.H b/util/xrb_ignition/ignition_condition.H index 14174cf5f..ce349705c 100644 --- a/util/xrb_ignition/ignition_condition.H +++ b/util/xrb_ignition/ignition_condition.H @@ -47,7 +47,120 @@ struct State template -void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real stepsize, +AMREX_FORCE_INLINE +void rk4_step(amrex::Real t, amrex::Real dt, + const amrex::Real (&f)[NEQ], const State& state, + void (*rhs)(amrex::Real, const amrex::Real (&)[NEQ], + const State&, amrex::Real (&)[NEQ]), + amrex::Real (&fout)[NEQ]) +{ + amrex::Real k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ]; + amrex::Real ftmp[NEQ]; + + rhs(t, f, state, k1); + + for (int n = 0; n < NEQ; ++n) + ftmp[n] = f[n] + 0.5_rt * dt * k1[n]; + rhs(t + 0.5_rt*dt, ftmp, state, k2); + + for (int n = 0; n < NEQ; ++n) + ftmp[n] = f[n] + 0.5_rt * dt * k2[n]; + rhs(t + 0.5_rt*dt, ftmp, state, k3); + + for (int n = 0; n < NEQ; ++n) + ftmp[n] = f[n] + dt * k3[n]; + rhs(t + dt, ftmp, state, k4); + + for (int n = 0; n < NEQ; ++n) + fout[n] = f[n] + (dt / 6.0_rt) * + (k1[n] + 2.0_rt*k2[n] + 2.0_rt*k3[n] + k4[n]); +} + + +template +void rk4_adaptive(amrex::Real t0, amrex::Real tmax, + amrex::Real dt_init, amrex::Real eps, + const amrex::Real (&f0)[NEQ], const State& state, + void (*rhs)(amrex::Real, const amrex::Real (&)[NEQ], + const State&, amrex::Real (&)[NEQ]), + amrex::Real (&f_new)[NEQ]) +{ + constexpr amrex::Real SAFETY = 0.9_rt; + constexpr amrex::Real GROW = -0.2_rt; // 1/(order+1) + constexpr amrex::Real SHRINK = -0.25_rt; // 1/order + constexpr amrex::Real TINY = 1.e-30_rt; + + amrex::Real t = t0; + amrex::Real dt = dt_init; + + for (int n = 0; n < NEQ; ++n) { + f_new[n] = f0[n]; + } + + amrex::Real y_big[NEQ], y_half[NEQ], y_tmp[NEQ], yscal[NEQ]; + + while (t < tmax) { + + if (t + dt > tmax) { + dt = tmax - t; + } + + // scaling for error norm + for (int n = 0; n < NEQ; ++n) { + yscal[n] = std::abs(f_new[n]) + TINY; + } + + amrex::Real errmax; + + for (;;) { + // one big step + rk4_step(t, dt, f_new, state, rhs, y_big); + + // two half steps + rk4_step(t, 0.5_rt*dt, f_new, state, rhs, y_tmp); + rk4_step(t + 0.5_rt*dt, 0.5_rt*dt, y_tmp, state, rhs, y_half); + + // compute error + errmax = 0.0_rt; + for (int n = 0; n < NEQ; ++n) { + amrex::Real err = std::abs(y_half[n] - y_big[n]) / yscal[n]; + errmax = std::max(errmax, err); + } + errmax /= eps; + + if (errmax <= 1.0_rt) { + break; + } + + // reduce timestep + amrex::Real dt_new = SAFETY * dt * std::pow(errmax, SHRINK); + dt = std::max(dt_new, 0.1_rt * dt); + + if (dt <= TINY) { + amrex::Abort("rk4_adaptive: timestep underflow"); + } + } + + // accept step + t += dt; + for (int n = 0; n < NEQ; ++n) + f_new[n] = y_half[n]; + + // increase timestep + amrex::Real dt_new; + if (errmax > 1.e-4_rt) { + dt_new = SAFETY * dt * std::pow(errmax, GROW); + } else { + dt_new = 5.0_rt * dt; + } + dt = std::min(dt_new, tmax - t); + } +} + + + +template +void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real dt_init, const amrex::Real (&f0)[NEQ], const State& state, void (*rhs)(const amrex::Real, const amrex::Real (&)[NEQ], const State&, amrex::Real (&)[NEQ]), @@ -59,14 +172,14 @@ void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real stepsiz ---------- t0: Initial time tmax: Final time - stepsize: Stepsize used for integration + dt_init: initial Stepsize used for integration f0: Initial value corresponding to t0, f(t0) = f0. rhs: Function that computes the derivative, df / dt */ // Set initial condition - amrex::Real dt = stepsize; + amrex::Real dt = dt_init; amrex::Real t = t0; for (int n = 0; n < NEQ; ++n) { @@ -79,7 +192,7 @@ void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real stepsiz while (t < tmax) { // Change stepsize so that we don't overshoot - if (t + stepsize > tmax) { + if (t + dt > tmax) { dt = tmax - t; } @@ -312,6 +425,8 @@ amrex::Real heating_rate(const eos_t& eos_state) { 44.0_rt * Fscn / (eos_state.T * eos_state.T * 1.0e-8_rt) + dFscn_dT) * epsilon_3a; + std::cout << std::scientific; + std::cout << "heat rate is = " << heat << std::endl; return heat; } @@ -351,6 +466,9 @@ amrex::Real cooling_rate(const amrex::Real yb, const eos_t& eos_state) { amrex::Real cool = 0.25_rt * eos_state.rho * (k_th + eos_state.T * d_kth_dT) / (yb * yb); + std::cout << std::scientific; + std::cout << "cool rate is = " << cool << std::endl; + return cool; } @@ -411,9 +529,9 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) { // t0=yt, tmax=yb, step_size=yt, f0={Tt, Ft, 0.0} // Integrates to find f_array corresponding to yb. - rk4(state.yt, state.yb, state.yt, - f0_array, state, derivs, - f_array); + rk4_adaptive(state.yt, state.yb, state.yt, 1.0e-6, + f0_array, state, derivs, + f_array); // Update solution Tb to State state.Tb = f_array[0]; @@ -526,7 +644,7 @@ void find_ignition_condition() { std::cout << "Metallicity is " << unit_test_rp::Z << std::endl; std::cout << "Accreted H fraction is " << unit_test_rp::X << std::endl; std::cout << "mdot/Edd is " << unit_test_rp::mdot << std::endl; - std::cout << "Depletion column depth is " << state.yd << std::endl; + std::cout << std::scientific << "Depletion column depth is " << state.yd << std::endl; // Now compute the base column depth: yb via root find. // This is also the ignition column depth @@ -554,6 +672,11 @@ void find_ignition_condition() { std::cout << "With ignition temperature: T = " << state.Tb << std::endl; std::cout << "With ignition pressure: P = " << state.g * state.yb << std::endl; std::cout << "With ignition density: ϱ = " << eos_state.rho << std::endl; + std::cout << "With relative height: z = " << state.z << std::endl; + + for (int n = 0; n < NumSpec; ++n) { + std::cout << "X(" << short_spec_names_cxx[n] << ") is " << eos_state.xn[n] << std::endl; + } } #endif From 6f48e87a535ba60ebf790ae2a10d00209ad0d1e6 Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 7 Jan 2026 16:45:30 -0500 Subject: [PATCH 4/5] use numerical differencing to catch temperature dependence on rho --- util/xrb_ignition/ignition_condition.H | 122 ++++++++++++++++--------- util/xrb_ignition/main.cpp | 2 - 2 files changed, 78 insertions(+), 46 deletions(-) diff --git a/util/xrb_ignition/ignition_condition.H b/util/xrb_ignition/ignition_condition.H index ce349705c..6630143a9 100644 --- a/util/xrb_ignition/ignition_condition.H +++ b/util/xrb_ignition/ignition_condition.H @@ -371,9 +371,24 @@ void derivs(const amrex::Real y, const amrex::Real (&f)[3], } -amrex::Real heating_rate(const eos_t& eos_state) { +amrex::Real heating_term(const amrex::Real T0, const eos_t& eos_state) { /* Evaluates the heating term -- triple alpha process + Use general fitting: + ϵ_3⍺ = 5.3e21 F_scn rho5^2 Y^3 / T8^3 exp(-44 / T8) + + Also apply optional enhancement factor due to quick + subsequent reaction rate: C12(p,g)N13(p,g)O14, where + + F_enh = 1 + Q12 / Q3a ≅ 1.9 + + only when X(H) > 1.e-3. + + Note that we technically want temperature derivative. + Although we technically have analytic form, but since + the EOS is e(T,p), then ϱ is also temperature dependent. + So to capture the effect, let's do a numerical difference. + Use logarithmnic diff instead of linear diff for better accuracy. */ // Q-value of C12 + 2p -> O14 ( C12(p,g)N13(p,g)O14 ) = 6.57 MeV @@ -382,21 +397,29 @@ amrex::Real heating_rate(const eos_t& eos_state) { constexpr amrex::Real Q3a = 7.274_rt; constexpr amrex::Real enhance_fac = 1.0 + Q12 / Q3a; + // Call EOS to find density. + auto tmp_state = eos_state; + tmp_state.T = T0; + eos(eos_input_tp, tmp_state); + // Find screening factor + constexpr int do_T_derivatives = 1; using number_t = std::conditional_t; - number_t temp = eos_state.T; - autodiff::seed(temp); + number_t temp = tmp_state.T; + if constexpr(do_T_derivatives) { + autodiff::seed(temp); + } plasma_state_t pstate{}; amrex::Array1D Y; for (int n = 0; n < NumSpec; ++n) { - Y(n+1) = eos_state.xn[n] * aion_inv[n]; + Y(n+1) = tmp_state.xn[n] * aion_inv[n]; } amrex::Real scor, scor2, dscor_dt, dscor2_dt; - fill_plasma_state(pstate, temp, eos_state.rho, Y); + fill_plasma_state(pstate, temp, tmp_state.rho, Y); constexpr auto scn_fac = scrn::calculate_screen_factor(2.0_rt, 4.0_rt, 2.0_rt, 4.0_rt); actual_screen(pstate, scn_fac, scor, dscor_dt); @@ -405,32 +428,35 @@ amrex::Real heating_rate(const eos_t& eos_state) { actual_screen(pstate, scn_fac2, scor2, dscor2_dt); amrex::Real Fscn = scor * scor2; - amrex::Real dFscn_dT = scor * dscor2_dt + dscor_dt * dscor2_dt; + // amrex::Real dFscn_dT = scor * dscor2_dt + dscor_dt * scor2; - // Maybe use actual triple alpha from pynucastro? But then - // eps_3a = 5.3e21 F_scn rho5^2 Y^3 / T8^3 exp(-44 / T8) - amrex::Real epsilon_3a = 5.3e21 * std::pow(eos_state.rho * 1.0e-5, 2) * - std::pow(eos_state.xn[Species::He4-1], 3) * std::pow(eos_state.T * 1.0e-8, -3) * - std::exp(-44.0 * 1.0e8 / eos_state.T); + // eps_3a = 5.3e21 Fscn rho5^2 Y^3 / T8^3 exp(-44 / T8) + amrex::Real heat = 5.3e11 * Fscn * std::pow(tmp_state.rho, 2) * + std::pow(tmp_state.xn[Species::He4-1], 3) * std::pow(tmp_state.T * 1.0e-8, -3) * + std::exp(-44.0 * 1.0e8 / tmp_state.T); // If there is enough hydrogen for the rapid proton capture reactions: // C12(p,g)N13(p,g)O14, i.e. when X(H1) > 1.0e-3, // then apply an additional enhancement factor due to enhancement, ~ 1.9 - if (eos_state.xn[Species::H1-1] > 1.0e-3_rt) { - epsilon_3a *= enhance_fac; + if (tmp_state.xn[Species::H1-1] > 1.0e-3_rt) { + heat *= enhance_fac; } // Heating term is temperature derivative of epsilon_3a including screening enhancement. - amrex::Real heat = (-3.0_rt * Fscn / eos_state.T + - 44.0_rt * Fscn / (eos_state.T * eos_state.T * 1.0e-8_rt) + - dFscn_dT) * epsilon_3a; + // Also note that since the EOS is e(T, p), this means ϱ is also temperature dependent. + // Need dϱ/dT, but eos doesn't seem to provide this. + // Let's use numerical jacobian instead of analytic, like settle version. + + // amrex::Real heat = (-3.0_rt * Fscn / tmp_state.T + + // 44.0_rt * Fscn / (tmp_state.T * tmp_state.T * 1.0e-8_rt) + + // dFscn_dT) * epsilon_3a; - std::cout << std::scientific; - std::cout << "heat rate is = " << heat << std::endl; return heat; } -amrex::Real cooling_rate(const amrex::Real yb, const eos_t& eos_state) { + +amrex::Real cooling_term(const amrex::Real T0, const amrex::Real yb, + const eos_t& eos_state) { /* Evaluates the cooling term. Get cooling rate by combining: @@ -446,28 +472,14 @@ amrex::Real cooling_rate(const amrex::Real yb, const eos_t& eos_state) { we should do a finite difference here. */ - // Find conductivity + // Call EOS to find conductivity auto tmp_state = eos_state; - actual_conductivity(tmp_state); - amrex::Real k_th = tmp_state.conductivity; - - // Find conductivity temperature dependence via finite difference + tmp_state.T = T0; + eos(eos_input_tp, tmp_state); - tmp_state.T = eos_state.T * 1.001_rt; actual_conductivity(tmp_state); - amrex::Real k_hi = tmp_state.conductivity; - - tmp_state.T = eos_state.T / 1.001_rt; - actual_conductivity(tmp_state); - amrex::Real k_lo = tmp_state.conductivity; - - amrex::Real d_kth_dT = (k_hi - k_lo) / (eos_state.T * (1.001_rt - 1.0_rt / 1.001_rt)); - - amrex::Real cool = 0.25_rt * eos_state.rho * - (k_th + eos_state.T * d_kth_dT) / (yb * yb); - - std::cout << std::scientific; - std::cout << "cool rate is = " << cool << std::endl; + amrex::Real k_th = tmp_state.conductivity; + amrex::Real cool = 0.25_rt * tmp_state.rho * tmp_state.T * k_th / (yb * yb); return cool; } @@ -555,18 +567,33 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) { eos_state.xn[Species::O14-1] = O14 * state.Z0; eos_state.xn[Species::O15-1] = O15 * state.Z0; - // Find density via Temperature and pressure - eos(eos_input_tp, eos_state); - // After Tb is solved, assume we have runaway when // heating rate is larger than cooling rate, i.e. // d eps_3a / dT = d eps_cool / dT - amrex::Real heat = heating_rate(eos_state); - amrex::Real cool = cooling_rate(state.yb, eos_state); - amrex::Real diff = (heat - cool) / cool; + // To find temperature derivative of heating (triple-alpha) and cooling term + // Do a logarithmic difference: de/dT = (e / T) dln(e)/dln(T) + // dln(e)/dln(T) = [ln(e2) - ln(e1)] / [ln(T2) - ln(T1)] = ln(e2/e1) / ln(T2/T1) + // This is more accurate than linear difference for expressions that involve exponentials. + + constexpr amrex::Real PERTURB_RATIO = 1.0e-3_rt; + amrex::Real dT = state.Tb * PERTURB_RATIO; + + amrex::Real heat = heating_term(state.Tb, eos_state); + amrex::Real heatp = heating_term(state.Tb + dT, eos_state); + amrex::Real dheat_dT = heat * std::log(heatp/heat) / + (state.Tb * std::log(1 + PERTURB_RATIO)); + + amrex::Real cool = cooling_term(state.Tb, state.yb, eos_state); + amrex::Real coolp = cooling_term(state.Tb + dT, state.yb, eos_state); + amrex::Real dcool_dT = cool * std::log(coolp/cool) / + (state.Tb * std::log(1 + PERTURB_RATIO)); + + amrex::Real diff = (dheat_dT - dcool_dT) / dcool_dT; std::cout << std::scientific; + std::cout << "heat rate is = " << dheat_dT << std::endl; + std::cout << "cool rate is = " << dcool_dT << std::endl; std::cout << "yb = " << state.yb << ", Tb = " << state.Tb << ", (heat-cool) / cool = " << diff @@ -650,6 +677,7 @@ void find_ignition_condition() { // This is also the ignition column depth // Assume that we work with log10(yb) + std::cout << "\n" << std::endl; std::cout << "Searching for ignition column depth..." << std::endl; find_yb(state); @@ -668,12 +696,18 @@ void find_ignition_condition() { // Find density via Temperature and pressure eos(eos_input_tp, eos_state); + std::cout << "\n" << std::endl; + std::cout << "------------------------ SOLUTION ------------------------" << std::endl; std::cout << "Solved ignition column depth: y = " << state.yb << std::endl; std::cout << "With ignition temperature: T = " << state.Tb << std::endl; std::cout << "With ignition pressure: P = " << state.g * state.yb << std::endl; std::cout << "With ignition density: ϱ = " << eos_state.rho << std::endl; std::cout << "With relative height: z = " << state.z << std::endl; + // Redshift factor used in settle + constexpr amrex::Real ZZ = 1.31; + + std::cout << "With recurrence time (hours): y/mdot = " << ZZ * state.yb / state.mdot / 3600.0 << std::endl; for (int n = 0; n < NumSpec; ++n) { std::cout << "X(" << short_spec_names_cxx[n] << ") is " << eos_state.xn[n] << std::endl; } diff --git a/util/xrb_ignition/main.cpp b/util/xrb_ignition/main.cpp index 23530071e..81c10a1dc 100644 --- a/util/xrb_ignition/main.cpp +++ b/util/xrb_ignition/main.cpp @@ -15,8 +15,6 @@ int main(int argc, char *argv[]) { amrex::Initialize(argc, argv); - std::cout << "starting the single zone burn..." << std::endl; - amrex::ParmParse ppa("amr"); init_unit_test(); From 701d4c8fa72e697019837bb6e37d06096b909652 Mon Sep 17 00:00:00 2001 From: zhi Date: Fri, 9 Jan 2026 17:26:43 -0500 Subject: [PATCH 5/5] fix codespell and powi --- util/xrb_ignition/ignition_condition.H | 51 ++++++++++---------------- 1 file changed, 19 insertions(+), 32 deletions(-) diff --git a/util/xrb_ignition/ignition_condition.H b/util/xrb_ignition/ignition_condition.H index 6630143a9..f04a9df41 100644 --- a/util/xrb_ignition/ignition_condition.H +++ b/util/xrb_ignition/ignition_condition.H @@ -52,8 +52,9 @@ void rk4_step(amrex::Real t, amrex::Real dt, const amrex::Real (&f)[NEQ], const State& state, void (*rhs)(amrex::Real, const amrex::Real (&)[NEQ], const State&, amrex::Real (&)[NEQ]), - amrex::Real (&fout)[NEQ]) -{ + amrex::Real (&fout)[NEQ]) { + // Single rk4 update. + amrex::Real k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ]; amrex::Real ftmp[NEQ]; @@ -83,8 +84,9 @@ void rk4_adaptive(amrex::Real t0, amrex::Real tmax, const amrex::Real (&f0)[NEQ], const State& state, void (*rhs)(amrex::Real, const amrex::Real (&)[NEQ], const State&, amrex::Real (&)[NEQ]), - amrex::Real (&f_new)[NEQ]) -{ + amrex::Real (&f_new)[NEQ]) { + // rk4 with adaptive timestepping + constexpr amrex::Real SAFETY = 0.9_rt; constexpr amrex::Real GROW = -0.2_rt; // 1/(order+1) constexpr amrex::Real SHRINK = -0.25_rt; // 1/order @@ -186,8 +188,7 @@ void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real dt_init f_new[n] = f0[n]; } - amrex::Real k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ]; - amrex::Real ftmp[NEQ], f[NEQ]; + amrex::Real f[NEQ]; while (t < tmax) { @@ -201,27 +202,7 @@ void rk4(const amrex::Real t0, const amrex::Real tmax, const amrex::Real dt_init f[n] = f_new[n]; } - rhs(t, f, state, k1); - - for (int n = 0; n < NEQ; ++n) { - ftmp[n] = f[n] + 0.5_rt * dt * k1[n]; - } - rhs(t + 0.5_rt * dt, ftmp, state, k2); - - for (int n = 0; n < NEQ; ++n) { - ftmp[n] = f[n] + 0.5_rt * dt * k2[n]; - } - rhs(t + 0.5_rt * dt, ftmp, state, k3); - - for (int n = 0; n < NEQ; ++n) { - ftmp[n] = f[n] + dt * k3[n]; - } - rhs(t + dt, ftmp, state, k4); - - for (int n = 0; n < NEQ; ++n) { - f_new[n] = f[n] + (dt / 6.0_rt) * - (k1[n] + 2.0_rt*k2[n] + 2.0_rt*k3[n] + k4[n]); - } + rk4_step(t, dt, f, state, rhs, f_new); // update time t += dt; @@ -431,8 +412,9 @@ amrex::Real heating_term(const amrex::Real T0, const eos_t& eos_state) { // amrex::Real dFscn_dT = scor * dscor2_dt + dscor_dt * scor2; // eps_3a = 5.3e21 Fscn rho5^2 Y^3 / T8^3 exp(-44 / T8) - amrex::Real heat = 5.3e11 * Fscn * std::pow(tmp_state.rho, 2) * - std::pow(tmp_state.xn[Species::He4-1], 3) * std::pow(tmp_state.T * 1.0e-8, -3) * + amrex::Real heat = 5.3e11 * Fscn * tmp_state.rho * tmp_state.rho * + amrex::Math::powi<3>(tmp_state.xn[Species::He4-1]) * + amrex::Math::powi<-3>(tmp_state.T * 1.0e-8) * std::exp(-44.0 * 1.0e8 / tmp_state.T); // If there is enough hydrogen for the rapid proton capture reactions: @@ -517,11 +499,11 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) { // If we were to compute kappa manually, we need eos. // i.e. we need two inputs. We have p=g * yt, so we need rho? - state.Tt = std::pow(2650.0 * state.Ft * state.yt, 0.25); + state.Tt = std::sqrt(std::sqrt(2650.0 * state.Ft * state.yt)); // Now set the lower boundary conditions: yb, Tb, and Fb // Here Fb is assumed to be given. Note that the original version - // includes a special term for compressional heating, if thats included + // includes a special term for compressional heating, if that is included // then we have to solve for Fb. However, it is sufficient to include // an extra constant heating term to account for compressional heating. // yb is what we're after, so we need Tb. @@ -539,8 +521,13 @@ amrex::Real yb_constraint_eq(amrex::Real log10yb_0, State& state) { amrex::Real f_array[NEQ]; amrex::Real f0_array[NEQ] = {state.Tt, state.Ft, 0.0_rt}; - // t0=yt, tmax=yb, step_size=yt, f0={Tt, Ft, 0.0} + // t0=yt, tmax=yb, dt=yt, f0={Tt, Ft, 0.0} // Integrates to find f_array corresponding to yb. + + // rk4(state.yt, state.yb, state.yt, + // f0_array, state, derivs, + // f_array); + rk4_adaptive(state.yt, state.yb, state.yt, 1.0e-6, f0_array, state, derivs, f_array);