diff --git a/extern/Microphysics b/extern/Microphysics index ab8f7f5419..72847c101c 160000 --- a/extern/Microphysics +++ b/extern/Microphysics @@ -1 +1 @@ -Subproject commit ab8f7f54191f271654a8c772e6aeeab4155a82bb +Subproject commit 72847c101c54ff5f17b143b066d6ba1871c70278 diff --git a/inputs/DTypeFront.toml b/inputs/DTypeFront.toml index f7d152a8bc..0d4fc32003 100644 --- a/inputs/DTypeFront.toml +++ b/inputs/DTypeFront.toml @@ -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 diff --git a/inputs/StromgrenRSLA.toml b/inputs/StromgrenRSLA.toml index c98072241e..90e6a90c97 100644 --- a/inputs/StromgrenRSLA.toml +++ b/inputs/StromgrenRSLA.toml @@ -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" diff --git a/inputs/StromgrenTempIndependentRecombination.toml b/inputs/StromgrenTempIndependentRecombination.toml index 3fb75d09ce..468ca60cd8 100644 --- a/inputs/StromgrenTempIndependentRecombination.toml +++ b/inputs/StromgrenTempIndependentRecombination.toml @@ -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" diff --git a/src/networks/network_header.template b/src/networks/network_header.template index f3742f61bd..4f76a14e14 100644 --- a/src/networks/network_header.template +++ b/src/networks/network_header.template @@ -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 ChemBandsHeader_{ @CHEM_BANDS@ }; #endif diff --git a/src/networks/photoionization/actual_rhs.H b/src/networks/photoionization/actual_rhs.H index b26509964a..8ffcbf8d01 100644 --- a/src/networks/photoionization/actual_rhs.H +++ b/src/networks/photoionization/actual_rhs.H @@ -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; @@ -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 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void actual_jac(const burn_t &state, MatrixType &jac) @@ -180,7 +184,9 @@ template 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; @@ -204,20 +210,17 @@ template 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 * @@ -227,21 +230,26 @@ template 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 diff --git a/src/problems/DTypeFront/CMakeLists.txt b/src/problems/DTypeFront/CMakeLists.txt index 415d0fd3fa..96a8815d94 100644 --- a/src/problems/DTypeFront/CMakeLists.txt +++ b/src/problems/DTypeFront/CMakeLists.txt @@ -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 @@ -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}) diff --git a/src/problems/OneZonePhotoionization/CMakeLists.txt b/src/problems/OneZonePhotoionization/CMakeLists.txt index b1192a44cc..35762bb8e3 100644 --- a/src/problems/OneZonePhotoionization/CMakeLists.txt +++ b/src/problems/OneZonePhotoionization/CMakeLists.txt @@ -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}) diff --git a/src/problems/StromgrenSphere/CMakeLists.txt b/src/problems/StromgrenSphere/CMakeLists.txt index 4c6e742512..4bd9c6bb15 100644 --- a/src/problems/StromgrenSphere/CMakeLists.txt +++ b/src/problems/StromgrenSphere/CMakeLists.txt @@ -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 @@ -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}) diff --git a/src/problems/StromgrenSphereRSLA/CMakeLists.txt b/src/problems/StromgrenSphereRSLA/CMakeLists.txt index 1ebe7d7964..158cffdc7c 100644 --- a/src/problems/StromgrenSphereRSLA/CMakeLists.txt +++ b/src/problems/StromgrenSphereRSLA/CMakeLists.txt @@ -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) if(AMReX_GPU_BACKEND MATCHES "CUDA") setup_target_for_cuda_compilation(${JOB_NAME}) diff --git a/src/radiation/photochemistry.hpp b/src/radiation/photochemistry.hpp index 42d6580d2a..b27cc28050 100644 --- a/src/radiation/photochemistry.hpp +++ b/src/radiation/photochemistry.hpp @@ -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::scalar0_index + nn) / spmasses[nn]; } +#ifdef SKIP_PHOTOCHEMFLUX + amrex::GpuArray 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; @@ -125,17 +132,19 @@ auto computePhotoChemistry(amrex::MultiFab &mf, const Real dt, const int stage, state(i, j, k, RadSystem::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; + 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;