Skip to content

Skip radiation number density during timestep control #2004

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
James471 wants to merge 6 commits into
base: development
Choose a base branch
from
+53 −25
Draft

Skip radiation number density during timestep control #2004

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
a668f7d
use Rosenbroke
chongchonghe Jun 25, 2026
c3e4a01
bump Microphysics
chongchonghe Jun 25, 2026
b7ffff0
Add flag to skip flux integration in microphysics
James471 Jun 25, 2026
acf2e53
Update microphysics to skip radiation convergence check
James471 Jun 26, 2026
f57ca28
Skip radiation number density in microphysics timestep control
James471 Jun 26, 2026
8e57b6a
Update microphysics
James471 Jun 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion extern/Microphysics
Open in desktop
Submodule Microphysics updated 9 files
+1 −1 CMakeLists.txt
+52 −4 integration/Rosenbrock/rosenbrock_integrator.H
+12 −0 integration/Rosenbrock/rosenbrock_type.H
+6 −0 integration/VODE/vode_dvhin.H
+1 −1 integration/VODE/vode_dvnlsd.H
+14 −0 integration/VODE/vode_dvode.H
+2 −0 integration/VODE/vode_type.H
+1 −1 integration/integrator_setup_sdc.H
+10 −4 integration/integrator_setup_strang.H
2 changes: 2 additions & 0 deletions inputs/DTypeFront.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,3 +66,5 @@ integrator.radiation_failure_tolerance = 0.5 # cm^-3; must exceed max N_g
integrator.rtol_spec = 1e-2
integrator.rtol_rad_num = 1e-2
integrator.rtol_enuc = 1e-2

integrator.rosenbrock_tableau = 3
1 change: 1 addition & 0 deletions inputs/StromgrenRSLA.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ integrator.jacobian = 1
integrator.renormalize_abundances = 0
integrator.subtract_internal_energy = 0
integrator.X_reject_buffer = 1e100
integrator.rosenbrock_tableau = 3

plotfile_interval = 500
plotfile_prefix = "StrmConstRecRSLA"
Expand Down
2 changes: 2 additions & 0 deletions inputs/StromgrenTempIndependentRecombination.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@ integrator.jacobian = 1
integrator.renormalize_abundances = 0
integrator.subtract_internal_energy = 0
integrator.X_reject_buffer = 1e100
# ROS2S is too slow for this problem
integrator.rosenbrock_tableau = 1

plotfile_interval = 50
plotfile_prefix = "StrmConstRec"
Expand Down
4 changes: 4 additions & 0 deletions src/networks/network_header.template
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,11 @@ constexpr int NumSpecTotal = NumSpec + NumSpecExtra;

#ifdef PHOTOCHEMISTRY
constexpr int NumChemBands = @NUM_CHEM_BANDS@;
#ifdef SKIP_PHOTOCHEMFLUX
constexpr int MicrophysicsNumRadVarsPerGroup = 1;
#else
constexpr int MicrophysicsNumRadVarsPerGroup = 2;
#endif
constexpr int NumRadEqs = NumChemBands * MicrophysicsNumRadVarsPerGroup;
constexpr amrex::GpuArray<double, NumChemBands + 1> ChemBandsHeader_{ @CHEM_BANDS@ };
#endif
Expand Down
18 changes: 13 additions & 5 deletions src/networks/photoionization/actual_rhs.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,9 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual_rhs(const burn_t &state, Ar
const amrex::Real n_HII = state.xn[2];
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
#ifndef SKIP_PHOTOCHEMFLUX
const amrex::Real n_gamma_flux = state.rn[1];
#endif

int molecular_process_switch = 1;
const amrex::Real n_H_total = n_HI + n_HII;
Expand Down Expand Up @@ -142,7 +144,9 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual_rhs(const burn_t &state, Ar
ydot(net_ienuc) = energy_source / state.rho;
}
ydot(5) = -photoionization_term;
#ifndef SKIP_PHOTOCHEMFLUX
ydot(6) = -chat_sigma_photo * n_HI * n_gamma_flux;
#endif
}

template <class MatrixType> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual_jac(const burn_t &state, MatrixType &jac)
Expand Down Expand Up @@ -180,7 +184,9 @@ template <class MatrixType> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual
// Gate for column-5 Jacobian entries: when rn[0] < 0 the clamp makes the RHS
// independent of rn[0], so all derivatives w.r.t. rn[0] are zero.
const amrex::Real n_gamma_active = (state.rn[0] >= 0.0_rt) ? 1.0_rt : 0.0_rt;
#ifndef SKIP_PHOTOCHEMFLUX
const amrex::Real n_gamma_flux = state.rn[1];
#endif

int molecular_process_switch = 1;
const amrex::Real n_H_total = n_HI + n_HII;
Expand All @@ -204,20 +210,17 @@ template <class MatrixType> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual
jac(1, 3) = -recombination_switch * alpha_rec * n_e;
jac(1, net_ienuc) = (d_collisional_ionization_term_dT - d_recombination_term_dT) * dTde;
jac(1, 5) = n_gamma_active * chat_sigma_photo * n_HI;
jac(1, 6) = 0.0_rt;

jac(2, 1) = -jac(1, 1);
jac(2, 2) = -jac(1, 2);
jac(2, 3) = -jac(1, 3);
jac(2, net_ienuc) = -jac(1, net_ienuc);
jac(2, 5) = -jac(1, 5);
jac(2, 6) = -jac(1, 6);
jac(3, 1) = jac(1, 1);
jac(3, 2) = jac(1, 2);
jac(3, 3) = jac(1, 3);
jac(3, net_ienuc) = jac(1, net_ienuc);
jac(3, 5) = jac(1, 5);
jac(3, 6) = jac(1, 6);

jac(net_ienuc, 1) = energy_switch * (-recombination_cooling_switch * Lambda_rec * n_HII - ion_ff_cooling_switch * Lambda_ion_ff * n_HII) / state.rho;
jac(net_ienuc, 2) = energy_switch *
Expand All @@ -227,21 +230,26 @@ template <class MatrixType> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual
jac(net_ienuc, 3) = energy_switch * (-recombination_cooling_switch * Lambda_rec * n_e - ion_ff_cooling_switch * Lambda_ion_ff * n_e) / state.rho;
jac(net_ienuc, net_ienuc) = energy_switch * (-d_recombination_cooling_dT - d_KI_cooling_dT - d_ion_ff_cooling_dT) * dTde / state.rho;
jac(net_ienuc, 5) = energy_switch * (photoheating_switch * Gamma_photo * chat_sigma_photo * n_HI) * n_gamma_active / state.rho;
jac(net_ienuc, 6) = 0.0_rt;

jac(5, 1) = 0.0_rt;
jac(5, 2) = -chat_sigma_photo * n_gamma;
jac(5, 3) = 0.0_rt;
jac(5, net_ienuc) = 0.0_rt;
jac(5, 5) = -n_gamma_active * chat_sigma_photo * n_HI;
jac(5, 6) = 0.0_rt;

#ifndef SKIP_PHOTOCHEMFLUX
jac(1, 6) = 0.0_rt;
jac(2, 6) = -jac(1, 6);
jac(3, 6) = jac(1, 6);
jac(net_ienuc, 6) = 0.0_rt;
jac(5, 6) = 0.0_rt;
jac(6, 1) = 0.0_rt;
jac(6, 2) = -chat_sigma_photo * n_gamma_flux;
jac(6, 3) = 0.0_rt;
jac(6, net_ienuc) = 0.0_rt;
jac(6, 5) = 0.0_rt;
jac(6, 6) = -n_HI * chat_sigma_photo;
#endif
}

#endif
5 changes: 3 additions & 2 deletions src/problems/DTypeFront/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@ if(AMReX_SPACEDIM GREATER_EQUAL 2)
# photoionization for this
# directory only
setup_target_for_microphysics_compilation(${microphysics_network_name}
"${CMAKE_CURRENT_BINARY_DIR}/")
"${CMAKE_CURRENT_BINARY_DIR}/"
"Rosenbrock")

# use the BEFORE keyword so that these files get priority in compilation for
# targets in this directory this is critical to ensure the correct Microphysics
Expand Down Expand Up @@ -91,7 +92,7 @@ if(AMReX_SPACEDIM GREATER_EQUAL 2)
../../radiation/photochemistry.cpp ${photoionization_sources})

# this will add #define CHEMISTRY
target_compile_definitions(${JOB_NAME} PUBLIC PHOTOCHEMISTRY)
target_compile_definitions(${JOB_NAME} PUBLIC PHOTOCHEMISTRY SKIP_PHOTOCHEMFLUX SKIP_RAD_NUM_CONV)

if(AMReX_GPU_BACKEND MATCHES "CUDA")
setup_target_for_cuda_compilation(${JOB_NAME})
Expand Down
2 changes: 1 addition & 1 deletion src/problems/OneZonePhotoionization/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ add_executable(
../../radiation/photochemistry.cpp ${photoionization_sources})

# this will add #define PHOTOCHEMISTRY
target_compile_definitions(${JOB_NAME} PUBLIC PHOTOCHEMISTRY)
target_compile_definitions(${JOB_NAME} PUBLIC PHOTOCHEMISTRY SKIP_PHOTOCHEMFLUX SKIP_RAD_NUM_CONV)

if(AMReX_GPU_BACKEND MATCHES "CUDA")
setup_target_for_cuda_compilation(${JOB_NAME})
Expand Down
5 changes: 3 additions & 2 deletions src/problems/StromgrenSphere/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@ if(AMReX_SPACEDIM EQUAL 3)
# photoionization for this
# directory only
setup_target_for_microphysics_compilation(${microphysics_network_name}
"${CMAKE_CURRENT_BINARY_DIR}/")
"${CMAKE_CURRENT_BINARY_DIR}/"
"Rosenbrock")

# use the BEFORE keyword so that these files get priority in compilation for
# targets in this directory this is critical to ensure the correct Microphysics
Expand Down Expand Up @@ -91,7 +92,7 @@ if(AMReX_SPACEDIM EQUAL 3)
../../radiation/photochemistry.cpp ${photoionization_sources})

# this will add #define PHOTOCHEMISTRY
target_compile_definitions(${JOB_NAME} PUBLIC PHOTOCHEMISTRY)
target_compile_definitions(${JOB_NAME} PUBLIC PHOTOCHEMISTRY SKIP_PHOTOCHEMFLUX SKIP_RAD_NUM_CONV)

if(AMReX_GPU_BACKEND MATCHES "CUDA")
setup_target_for_cuda_compilation(${JOB_NAME})
Expand Down
2 changes: 1 addition & 1 deletion src/problems/StromgrenSphereRSLA/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ if(AMReX_SPACEDIM GREATER_EQUAL 2)
../../radiation/photochemistry.cpp ${photoionization_sources})

# this will add #define CHEMISTRY
target_compile_definitions(${JOB_NAME} PUBLIC PHOTOCHEMISTRY)
target_compile_definitions(${JOB_NAME} PUBLIC PHOTOCHEMISTRY SKIP_PHOTOCHEMFLUX SKIP_RAD_NUM_CONV)

@gemini-code-assist gemini-code-assist Bot Jun 26, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

Since StromgrenRSLA.toml has been updated to use integrator.rosenbrock_tableau = 3, the setup_target_for_microphysics_compilation call at the beginning of this file (around line 6) should also be updated to pass "Rosenbrock" as the third argument, similar to the changes made in DTypeFront and StromgrenSphere. Otherwise, the Rosenbrock integrator will not be compiled for this target, leading to build or runtime failures.


if(AMReX_GPU_BACKEND MATCHES "CUDA")
setup_target_for_cuda_compilation(${JOB_NAME})
Expand Down
35 changes: 22 additions & 13 deletions src/radiation/photochemistry.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,10 +80,17 @@ auto computePhotoChemistry(amrex::MultiFab &mf, const Real dt, const int stage,
for (int nn = 0; nn < NumSpec; ++nn) {
photochemstate.xn[nn] = state(i, j, k, RadSystem<problem_t>::scalar0_index + nn) / spmasses[nn];
}
#ifdef SKIP_PHOTOCHEMFLUX
amrex::GpuArray<Real, NumChemBands> n_gamma_initial{};
#endif
for (int nn = 0; nn < NumChemBands; ++nn) {
photochemstate.rn[0 + MicrophysicsNumRadVarsPerGroup * nn] =
state(i, j, k, firstChemIndex + Physics_NumVars::numRadVarsPerGroup * nn) * invChemBandQuanta[nn];
const Real n_gamma = state(i, j, k, firstChemIndex + Physics_NumVars::numRadVarsPerGroup * nn) * invChemBandQuanta[nn];
photochemstate.rn[0 + MicrophysicsNumRadVarsPerGroup * nn] = n_gamma;
#ifdef SKIP_PHOTOCHEMFLUX
n_gamma_initial[nn] = n_gamma;
#else
photochemstate.rn[1 + MicrophysicsNumRadVarsPerGroup * nn] = 1.0_rt;
#endif
}
photochemstate.rho = rho;
photochemstate.e = Eint / rho;
Expand Down Expand Up @@ -125,17 +132,19 @@ auto computePhotoChemistry(amrex::MultiFab &mf, const Real dt, const int stage,
state(i, j, k, RadSystem<problem_t>::scalar0_index + nn) = photochemstate.xn[nn] * spmasses[nn];
}
for (int nn = 0; nn < NumChemBands; ++nn) {
state(i, j, k, firstChemIndex + Physics_NumVars::numRadVarsPerGroup * nn) =
photochemstate.rn[0 + MicrophysicsNumRadVarsPerGroup * nn] * chemBandQuanta[nn];
state(i, j, k, firstChemFxIndex + Physics_NumVars::numRadVarsPerGroup * nn) =
photochemstate.rn[1 + MicrophysicsNumRadVarsPerGroup * nn] *
state(i, j, k, firstChemFxIndex + Physics_NumVars::numRadVarsPerGroup * nn);
state(i, j, k, firstChemFyIndex + Physics_NumVars::numRadVarsPerGroup * nn) =
photochemstate.rn[1 + MicrophysicsNumRadVarsPerGroup * nn] *
state(i, j, k, firstChemFyIndex + Physics_NumVars::numRadVarsPerGroup * nn);
state(i, j, k, firstChemFzIndex + Physics_NumVars::numRadVarsPerGroup * nn) =
photochemstate.rn[1 + MicrophysicsNumRadVarsPerGroup * nn] *
state(i, j, k, firstChemFzIndex + Physics_NumVars::numRadVarsPerGroup * nn);
const Real n_gamma_final = photochemstate.rn[0 + MicrophysicsNumRadVarsPerGroup * nn];
state(i, j, k, firstChemIndex + Physics_NumVars::numRadVarsPerGroup * nn) = n_gamma_final * chemBandQuanta[nn];
#ifdef SKIP_PHOTOCHEMFLUX
const Real flux_ratio = (n_gamma_initial[nn] > 0.0_rt) ? (n_gamma_final / n_gamma_initial[nn]) : 0.0_rt;

@gemini-code-assist gemini-code-assist Bot Jun 26, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

Using 0.0_rt as the threshold for n_gamma_initial[nn] can lead to floating-point overflow (to infinity) when n_gamma_initial[nn] is extremely small (e.g., subnormal values close to underflow) and n_gamma_final is clamped to small_x. If the initial flux is 0.0, multiplying it by infinity produces NaN, which will propagate and corrupt the simulation state. Using a safe, physically negligible threshold like 1e-100_rt prevents this overflow and subsequent NaN propagation.

const Real flux_ratio = (n_gamma_initial[nn] > 1e-100_rt) ? (n_gamma_final / n_gamma_initial[nn]) : 0.0_rt;

state(i, j, k, firstChemFxIndex + Physics_NumVars::numRadVarsPerGroup * nn) *= flux_ratio;
state(i, j, k, firstChemFyIndex + Physics_NumVars::numRadVarsPerGroup * nn) *= flux_ratio;
state(i, j, k, firstChemFzIndex + Physics_NumVars::numRadVarsPerGroup * nn) *= flux_ratio;
#else
const Real flux_norm = photochemstate.rn[1 + MicrophysicsNumRadVarsPerGroup * nn];
state(i, j, k, firstChemFxIndex + Physics_NumVars::numRadVarsPerGroup * nn) *= flux_norm;
state(i, j, k, firstChemFyIndex + Physics_NumVars::numRadVarsPerGroup * nn) *= flux_norm;
state(i, j, k, firstChemFzIndex + Physics_NumVars::numRadVarsPerGroup * nn) *= flux_norm;
#endif
}
// Quokka uses rho*eint
const Real dEint = (photochemstate.e * photochemstate.rho) - Eint;
Expand Down
Loading