fix(photochem): exclude flux from VODE convergence; hard-code atol for DTypeFront; add docs

inputs/DTypeFront.toml

@@ -13,15 +13,14 @@ amr.v              = 0       # verbosity in Amr
 # *****************************************************************
 # Resolution and refinement
 # *****************************************************************
-amr.n_cell          = [64, 64, 64]
+amr.n_cell          = [8, 8, 8]
 amr.max_level       = 0     # number of levels = max_level + 1
-amr.blocking_factor = 32  # grid size must be divisible by this
+amr.blocking_factor = 8  # grid size must be divisible by this
 amr.max_grid_size = 128  # at least 128 for GPUs
 do_reflux = 0
 do_subcycle = 0
 suppress_output = 0
 
-max_timesteps = 10000
 stop_time = 2.9e12
 
 stromgen.primary_species_1 = 1.0e-10
@@ -39,7 +38,20 @@ integrator.subtract_internal_energy = 0
 integrator.X_reject_buffer = 1e100
 integrator.ode_max_steps = 50000
 
-plotfile_interval = 1000000
+plotfile_interval = -1
 plotfile_prefix = "DType"
 
 network.allow_mixed_cell_cooling = 0
+
+
+max_timesteps = 10
+
+integrator.atol_spec = 1e-3   # cm^-3; 1e-4 × n_H = physical negligibility floor
+integrator.rtol_spec = 1e-2
+integrator.atol_rad_num = 1e-2  # cm^-3; must exceed N_gamma_floor = Erad_floor/E_photon ~ 3.5e-4 so dark cells return in one VODE step
+integrator.radiation_failure_tolerance = 0.1  # N_gamma can dip to -atol_rad_num = -0.01 in cells where photons are fully absorbed; must exceed atol_rad_num
+integrator.rtol_rad_num = 1e-2
+integrator.atol_rad_flux = 1e-2  # normalized flux (dimensionless); crossover y_cross = atol/rtol = 1, always absolute control
+integrator.rtol_rad_flux = 1e-2
+integrator.atol_enuc = 1.24e8  # erg/g; ~2.5% of e at T=37 K neutral gas
+integrator.rtol_enuc = 1e-2

src/networks/photoionization/actual_rhs.H

@@ -100,7 +100,8 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual_rhs(const burn_t &state, Ar
 	const amrex::Real n_e = state.xn[0];
 	const amrex::Real n_HI = state.xn[1];
 	const amrex::Real n_HII = state.xn[2];
-	const amrex::Real n_gamma = state.rn[0];
+	const amrex::Real n_gamma =
+	    amrex::max(state.rn[0], 0.0_rt); // clamp at zero: prevents sign reversal if VODE steps rn[0] slightly negative within tolerance
 	const amrex::Real n_gamma_flux = state.rn[1];
 
 	int molecular_process_switch = 1;
@@ -174,7 +175,8 @@ template <class MatrixType> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual
 	const amrex::Real n_e = state.xn[0];
 	const amrex::Real n_HI = state.xn[1];
 	const amrex::Real n_HII = state.xn[2];
-	const amrex::Real n_gamma = state.rn[0];
+	const amrex::Real n_gamma =
+	    amrex::max(state.rn[0], 0.0_rt); // clamp at zero: prevents sign reversal if VODE steps rn[0] slightly negative within tolerance
 	const amrex::Real n_gamma_flux = state.rn[1];
 
 	int molecular_process_switch = 1;

src/problems/DTypeFront/testDTypeFront.cpp

@@ -60,7 +60,8 @@ template <> struct Physics_Traits<DTypeFront> {
 
 template <> struct RadSystem_Traits<DTypeFront> {
 	static constexpr double c_hat_over_c = c_hat / C::c_light;
-	static constexpr double Erad_floor = 1e-99;
+	static constexpr double Erad_floor =
+	    C::a_rad * 1.0; // a_rad * (1 K)^4; N_gamma_floor ~ 3.5e-4 cm^-3 < atol_rad_num, negligible ionization over 10 steps
 	static constexpr int beta_order = 0;
 	static constexpr auto ChemBands() { return ChemBandsHeader_; }
 };