From 5577a4e1d511b48f128831e46362c7f529078c01 Mon Sep 17 00:00:00 2001 From: James471 Date: Mon, 2 Mar 2026 23:27:33 +1100 Subject: [PATCH 1/8] Bug fix: Last character of eos parameters file is not read --- util/build_scripts/write_probin.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/util/build_scripts/write_probin.py b/util/build_scripts/write_probin.py index 58c4a140cf..3dcebe5d40 100755 --- a/util/build_scripts/write_probin.py +++ b/util/build_scripts/write_probin.py @@ -62,7 +62,10 @@ def get_next_line(fin): line = fin.readline() pos = line.find("#") - return line[:pos] + if pos > -1: + return line[:pos] + else: + return line def parse_param_file(params_list, param_file): From c259dcb4f7aeec75c963a5de3593167330b7f78d Mon Sep 17 00:00:00 2001 From: James471 Date: Mon, 2 Mar 2026 23:36:23 +1100 Subject: [PATCH 2/8] WIP: Implement radiation integration in VODE --- integration/VODE/vode_dvhin.H | 20 ++++++++- integration/VODE/vode_dvode.H | 20 +++++++++ integration/VODE/vode_dvstep.H | 60 +++++++++++++++++++++++++++ integration/VODE/vode_type.H | 8 +++- integration/_parameters | 24 +++++++++++ integration/integrator_data.H | 14 +++++++ integration/integrator_setup_strang.H | 40 ++++++++++++++++++ integration/integrator_type.H | 6 +++ integration/integrator_type_strang.H | 12 ++++++ interfaces/burn_type.H | 24 ++++++++++- interfaces/eos.H | 8 +++- 11 files changed, 230 insertions(+), 6 deletions(-) diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H index 96c18b5756..180fa267ef 100644 --- a/integration/VODE/vode_dvhin.H +++ b/integration/VODE/vode_dvhin.H @@ -45,7 +45,25 @@ void dvhin (BurnT& state, DvodeT& vstate, amrex::Real& H0, int& NITER, int& IER) amrex::Real HUB = PT1 * TDIST; for (int i = 1; i <= int_neqs; ++i) { - const amrex::Real atol = i <= NumSpec ? vstate.atol_spec : vstate.atol_enuc; + // TODO(#3)(james471): This can be optimized to reduce branching + amrex::Real atol{}; +#if defined (PHOTOCHEMISTRY) + if (i <= NumSpec) { + atol = vstate.atol_spec; + } else if (i == NumSpec + 1) { + atol = vstate.atol_enuc; + } else if ((i - NumSpec - 2) % MicrophysicsNumRadVarsPerGroup == 0) { + atol = vstate.atol_rad_num; + } else { + atol = vstate.atol_rad_flux; + } +#else + if (i <= NumSpec) { + atol = vstate.atol_spec; + } else { + atol = vstate.atol_enuc; + } +#endif const amrex::Real DELYI = PT1 * std::abs(vstate.yh(i,1)) + atol; const amrex::Real AFI = std::abs(vstate.yh(i,2)); if (AFI * HUB > DELYI) { diff --git a/integration/VODE/vode_dvode.H b/integration/VODE/vode_dvode.H index 361b7cba50..3dc8a37e84 100644 --- a/integration/VODE/vode_dvode.H +++ b/integration/VODE/vode_dvode.H @@ -81,6 +81,16 @@ int dvode (BurnT& state, DvodeT& vstate) } vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc; vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1); +#ifdef PHOTOCHEMISTRY + for (int i = 1; i <= NumChemActiveRadGroups; i+=MicrophysicsNumRadVarsPerGroup) { + int rad_n_index = NumSpec + 2 + (i-1) * MicrophysicsNumRadVarsPerGroup; + vstate.ewt(rad_n_index) = vstate.rtol_rad_num * std::abs(vstate.yh(rad_n_index,1)) + vstate.atol_rad_num; + for (int j = 2; j <= MicrophysicsNumRadVarsPerGroup; ++j){ + int rad_flux_index = rad_n_index + j - 1; + vstate.ewt(rad_flux_index) = vstate.rtol_rad_flux * std::abs(vstate.yh(rad_flux_index,1)) + vstate.atol_rad_flux; + } + } +#endif // Call DVHIN to set initial step size H0 to be attempted. H0 = 0.0_rt; @@ -154,6 +164,16 @@ int dvode (BurnT& state, DvodeT& vstate) } vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc; vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1); +#ifdef PHOTOCHEMISTRY + for (int i = 1; i <= NumChemActiveRadGroups; i+=MicrophysicsNumRadVarsPerGroup) { + int rad_n_index = NumSpec + 2 + (i-1) * MicrophysicsNumRadVarsPerGroup; + vstate.ewt(rad_n_index) = vstate.rtol_rad_num * std::abs(vstate.yh(rad_n_index,1)) + vstate.atol_rad_num; + for (int j = 2; j <= MicrophysicsNumRadVarsPerGroup; ++j){ + int rad_flux_index = rad_n_index + j - 1; + vstate.ewt(rad_flux_index) = vstate.rtol_rad_flux * std::abs(vstate.yh(rad_flux_index,1)) + vstate.atol_rad_flux; + } + } +#endif } else { diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index 18faa5d420..964098955f 100644 --- a/integration/VODE/vode_dvstep.H +++ b/integration/VODE/vode_dvstep.H @@ -148,6 +148,7 @@ int dvstep (BurnT& state, DvodeT& vstate) } #endif #ifdef SDC +// TODO(#3)(james471): Check what needs to be done for radiation variables amrex::Real rho_old = state.rho_orig + TOLD * state.ydot_a[SRHO]; for (int i = 1; i <= int_neqs; ++i) { y_save(i) = vstate.y(i); @@ -234,6 +235,7 @@ int dvstep (BurnT& state, DvodeT& vstate) } #endif + // TODO(#2)(james471): Add radiation checks here for (int i = 1; i <= NumSpec; ++i) { #ifdef STRANG @@ -308,6 +310,64 @@ int dvstep (BurnT& state, DvodeT& vstate) } + for (int i = net_ienuc + 1; i <= int_neqs; i+=2) { + +#ifdef STRANG + // Constrain radiation quantities such that they don't change by more than a certain + // factor in a given step (as long as their initial value was large enough + // for us to worry about). + + if (std::abs(y_save(i)) > integrator_rp::X_reject_buffer * vstate.atol_rad_num && + std::abs(vstate.y(i)) > integrator_rp::X_reject_buffer * vstate.atol_rad_num && + (std::abs(vstate.y(i)) > vode_increase_change_factor * std::abs(y_save(i)) || + std::abs(vstate.y(i)) < vode_decrease_change_factor * std::abs(y_save(i)))) { +#ifdef MICROPHYSICS_DEBUG +#ifndef AMREX_USE_GPU + std::cout << "rejecting step based on species " << i << " from " << y_save(i) << " to " << vstate.y(i) << std::endl; +#endif +#endif + valid_update = false; + break; + } + + if (std::abs(y_save(i+1)) > integrator_rp::X_reject_buffer * vstate.atol_rad_flux && + std::abs(vstate.y(i+1)) > integrator_rp::X_reject_buffer * vstate.atol_rad_flux && + (std::abs(vstate.y(i+1)) > vode_increase_change_factor * std::abs(y_save(i+1)) || + std::abs(vstate.y(i+1)) < vode_decrease_change_factor * std::abs(y_save(i+1)))) { +#ifdef MICROPHYSICS_DEBUG +#ifndef AMREX_USE_GPU + std::cout << "rejecting step based on species " << i+1 << " from " << y_save(i+1) << " to " << vstate.y(i+1) << std::endl; +#endif +#endif + valid_update = false; + break; + } + + // Constrain abundances such that they are not negative (within a tolerance) + // or greater than one (within a tolerance). + if (vstate.y(i) < -radiation_failure_tolerance) { + valid_update = false; + break; + } + + // Don't enforce the condition below if + // vstate.y contains number densities + if (! integrator_rp::use_number_densities) { + if (vstate.y(i) > 1.0_rt + radiation_failure_tolerance) { + valid_update = false; + break; + } + } + +#endif + +#ifdef SDC + // Abort the program. This is unimplemented + amrex::Abort("SDC does not yet support radiation updates"); +#endif + + } + // The corrector has converged (NFLAG = 0). The local error test is // made and control passes to statement 500 if it fails. diff --git a/integration/VODE/vode_type.H b/integration/VODE/vode_type.H index 3c246bfb50..8f389cfa36 100644 --- a/integration/VODE/vode_type.H +++ b/integration/VODE/vode_type.H @@ -16,8 +16,8 @@ constexpr amrex::Real CCMXJ = 0.2e0_rt; constexpr amrex::Real HMIN = 0.0_rt; -// For each species whose abundance is above a certain threshold, we do -// not allow its mass fraction to change by more than a certain amount +// For each species or radiation quantity whose value is above a certain threshold, +// we do not allow its mass fraction to change by more than a certain amount // in any integration step. constexpr amrex::Real vode_increase_change_factor = 4.0_rt; constexpr amrex::Real vode_decrease_change_factor = 0.25_rt; @@ -135,6 +135,10 @@ struct dvode_t // Tolerances amrex::Real rtol_spec, atol_spec; amrex::Real rtol_enuc, atol_enuc; +#ifdef PHOTOCHEMISTRY + amrex::Real rtol_rad_num, atol_rad_num; + amrex::Real rtol_rad_flux, atol_rad_flux; +#endif // Local time and integration end time amrex::Real t, tout; diff --git a/integration/_parameters b/integration/_parameters index a11b5f3ea8..9f2ccb04d6 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -34,6 +34,18 @@ atol_spec real 1.e-8 # absolute tolerance for energy atol_enuc real 1.e-6 +# relative tolerance for radiation density +rtol_rad_num real 1.e-12 + +# relative teolerance for radiation flux +rtol_rad_flux real 1.e-6 + +# absolute tolerance for radiation density +atol_rad_num real 1.e-8 + +# absolute tolerance for radiation flux +atol_rad_flux real 1.e-6 + # Whether to renormalize the mass fractions at each step in the evolution # so that they sum to unity. renormalize_abundances bool 0 @@ -89,6 +101,18 @@ retry_atol_spec real -1 # absolute tolerance for energy on retry retry_atol_enuc real -1 +# relative tolerance for radiation density on retry +retry_rtol_rad_num real -1 + +# relative tolerance for radiation flux on retry +retry_rtol_rad_flux real -1 + +# absolute tolerance for radiation density on retry +retry_atol_rad_num real -1 + +# absolute tolerance for radiation flux on retry +retry_atol_rad_flux real -1 + # in the clean_state process, do we clip the species such that they # are in [0, 1]? do_species_clip bool 1 diff --git a/integration/integrator_data.H b/integration/integrator_data.H index 2d4f7795cd..28950af3ec 100644 --- a/integration/integrator_data.H +++ b/integration/integrator_data.H @@ -5,12 +5,20 @@ // Define the size of the ODE system that VODE will integrate +#ifndef PHOTOCHEMISTRY constexpr int INT_NEQS = NumSpec + 1; +#else +constexpr int INT_NEQS = NumRadEqs + NumSpec + 1; +#endif // We will use this parameter to determine if a given species // abundance is unreasonably small or large (each X must satisfy // -failure_tolerance <= X <= 1.0 + failure_tolerance). constexpr amrex::Real species_failure_tolerance = 1.e-2_rt; +// We will use this parameter to determine if +// a given radiation quantity (number density or flux) is +// unreasonably small (each quantity must satisfy X >= -radiation_failure_tolerance). +constexpr amrex::Real radiation_failure_tolerance = 1.e-2_rt; #ifdef NSE constexpr int MIN_NSE_BAILOUT_STEPS = 10; @@ -38,6 +46,9 @@ constexpr int integrator_neqs () #ifdef STRANG if constexpr (has_xn::value) { int_neqs += NumSpec; +#ifdef PHOTOCHEMISTRY + int_neqs += NumRadEqs; +#endif } if constexpr (has_energy::value) { @@ -45,6 +56,9 @@ constexpr int integrator_neqs () } #else int_neqs = NumSpec + 1; +#ifdef PHOTOCHEMISTRY + int_neqs += NumRadEqs; +#endif #endif return int_neqs; diff --git a/integration/integrator_setup_strang.H b/integration/integrator_setup_strang.H index 66df898327..bd97c7ae23 100644 --- a/integration/integrator_setup_strang.H +++ b/integration/integrator_setup_strang.H @@ -20,6 +20,7 @@ struct state_backup_t { amrex::Real e_in{}; #ifndef AMREX_USE_GPU amrex::Real xn_in[NumSpec]{}; + amrex::Real rn_in[NumRadEqs]{}; #endif }; @@ -38,6 +39,13 @@ IntegratorT integrator_setup (BurnT& state, amrex::Real dt, bool is_retry) int_state.rtol_spec = integrator_rp::rtol_spec; // mass fractions int_state.rtol_enuc = integrator_rp::rtol_enuc; // energy generated +#ifdef PHOTOCHEMISTRY + int_state.atol_rad_num = integrator_rp::atol_rad_num; // radiation number density + int_state.atol_rad_flux = integrator_rp::atol_rad_flux; // radiation flux + + int_state.rtol_rad_num = integrator_rp::rtol_rad_num; // radiation number density + int_state.rtol_rad_flux = integrator_rp::rtol_rad_flux; // radiation flux +#endif } else { int_state.atol_spec = integrator_rp::retry_atol_spec > 0 ? integrator_rp::retry_atol_spec : integrator_rp::atol_spec; // mass fractions @@ -48,6 +56,16 @@ IntegratorT integrator_setup (BurnT& state, amrex::Real dt, bool is_retry) integrator_rp::retry_rtol_spec : integrator_rp::rtol_spec; // mass fractions int_state.rtol_enuc = integrator_rp::retry_rtol_enuc > 0 ? integrator_rp::retry_rtol_enuc : integrator_rp::rtol_enuc; // energy generated +#ifdef PHOTOCHEMISTRY + int_state.atol_rad_num = integrator_rp::retry_atol_rad_num > 0 ? + integrator_rp::retry_atol_rad_num : integrator_rp::atol_rad_num; // radiation number density + int_state.atol_rad_flux = integrator_rp::retry_atol_rad_flux > 0 ? + integrator_rp::retry_atol_rad_flux : integrator_rp::atol_rad_flux; // radiation flux + int_state.rtol_rad_num = integrator_rp::retry_rtol_rad_num > 0 ? + integrator_rp::retry_rtol_rad_num : integrator_rp::rtol_rad_num; // radiation number density + int_state.rtol_rad_flux = integrator_rp::retry_rtol_rad_flux > 0 ? + integrator_rp::retry_rtol_rad_flux : integrator_rp::rtol_rad_flux; // radiation flux +#endif } // set the Jacobian type @@ -107,6 +125,9 @@ state_backup_t integrator_backup (const BurnT& state) { for (int n = 0; n < NumSpec; ++n) { state_save.xn_in[n] = state.xn[n]; } + for (int n = 0; n < NumRadEqs; ++n) { + state_save.rn_in[n] = state.rn[n]; + } #endif state_save.T_in = state.T; state_save.e_in = state.e; @@ -180,6 +201,11 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state, } } } + for (int n = net_ienuc+1; n <= net_ienuc+NumRadEqs; n += MicrophysicsNumRadVarsPerGroup) { + if (int_state.y(n) < -radiation_failure_tolerance) { + state.success = false; + } + } #ifndef AMREX_USE_GPU if (integrator_rp::burner_verbose) { @@ -210,6 +236,13 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state, for (const double X : state_save.xn_in) { std::cout << std::scientific << std::setprecision(16) << X << " "; } +#ifdef PHOTOCHEMISTRY + std::cout << std::endl; + std::cout << "rn start = "; + for (const double X : state_save.rn_in) { + std::cout << std::scientific << std::setprecision(16) << X << " "; + } +#endif std::cout << std::endl; std::cout << "dens current = " << std::setprecision(16) << state.rho << std::endl; std::cout << "temp current = " << std::setprecision(16) << state.T << std::endl; @@ -217,6 +250,13 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state, for (const double X : state.xn) { std::cout << std::scientific << std::setprecision(16) << X << " "; } +#ifdef PHOTOCHEMISTRY + std::cout << std::endl; + std::cout << "rn current = "; + for (const double X : state.rn) { + std::cout << std::scientific << std::setprecision(16) << X << " "; + } +#endif std::cout << std::endl; std::cout << "energy generated = " << state.e << std::endl; #elif defined(ALLOW_GPU_PRINTF) diff --git a/integration/integrator_type.H b/integration/integrator_type.H index 8c6876adea..fe3d4081b5 100644 --- a/integration/integrator_type.H +++ b/integration/integrator_type.H @@ -34,6 +34,7 @@ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void update_thermodynamics (BurnT& state, const I& int_state) { + // TODO(#3)(james471): Add radiation check for aux thermo #ifdef AUX_THERMO // we are handling the thermodynamics via the aux quantities, which // are stored in the burn_t here, so we need to update those based on @@ -58,6 +59,11 @@ void update_thermodynamics (BurnT& state, const I& int_state) state.xn[n-1] = int_state.y(n); } state.e = int_state.y(net_ienuc); +#ifdef PHOTOCHEMISTRY + for (int n = 1; n <= NumRadEqs; ++n) { + state.rn[n-1] = int_state.y(net_ienuc + n); + } +#endif if (integrator_rp::scale_system) { state.e *= state.e_scale; diff --git a/integration/integrator_type_strang.H b/integration/integrator_type_strang.H index 7eb3c6d60f..8fb8f15d35 100644 --- a/integration/integrator_type_strang.H +++ b/integration/integrator_type_strang.H @@ -18,6 +18,12 @@ void integrator_to_burn (const I& int_state, BurnT& state) state.e = int_state.y(net_ienuc); +#ifdef PHOTOCHEMISTRY + for (int n = 1; n <= NumRadEqs; ++n) { + state.rn[n-1] = int_state.y(net_ienuc + n); + } +#endif + if (integrator_rp::scale_system) { state.e *= state.e_scale; } @@ -38,6 +44,12 @@ void burn_to_integrator (const BurnT& state, I& int_state) int_state.y(net_ienuc) = state.e; +#ifdef PHOTOCHEMISTRY + for (int n = 1; n <= NumRadEqs; ++n) { + int_state.y(net_ienuc + n) = state.rn[n-1]; + } +#endif + if (integrator_rp::scale_system) { int_state.y(net_ienuc) /= state.e_scale; } diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H index 316053530d..8665548147 100644 --- a/interfaces/burn_type.H +++ b/interfaces/burn_type.H @@ -24,7 +24,11 @@ using namespace network_rp; // For Strang evolution, this will be the number of equations we are // evolving. For simplified-SDC, we will need neqs when getting the // reaction sources from the network's RHS. +#if defined (PHOTOCHEMISTRY) +const int neqs = 1 + NumSpec + NumRadEqs; +#else const int neqs = 1 + NumSpec; +#endif // Index of the energy variable in the work arrays. @@ -43,7 +47,11 @@ using YdotNetArray1D = amrex::Array1D; // carry momentum as an unevolved variable const int SFS = 0; +#if defined (PHOTOCHEMISTRY) +const int SEINT = NumSpec + NumRadEqs; +#else const int SEINT = NumSpec; +#endif // the following are not evolved const int SFX = SEINT + 1; const int SRHO = SFX + NumAux; // this is SFS + NumSpec if NumAux = 0; @@ -76,6 +84,9 @@ struct burn_t #if NAUX_NET > 0 amrex::Real aux[NumAux]{}; #endif +#ifdef PHOTOCHEMISTRY + amrex::Real rn[NumRadEqs]{}; +#endif // scaling used to make e we integrate dimensionless amrex::Real e_scale{}; @@ -158,7 +169,7 @@ struct burn_t }; - +// TODO(james471): Add radiation variables to the printout inline std::ostream& operator<< (std::ostream& o, burn_t const& burn_state) { @@ -233,6 +244,11 @@ void eos_to_burn (const T& eos_state, BurnT& burn_state) burn_state.aux[n] = eos_state.aux[n]; } #endif +#ifdef PHOTOCHEMISTRY + for (int n = 0; n < NumRadEqs; ++n) { + burn_state.rn[n] = eos_state.rn[n]; + } +#endif burn_state.cv = eos_state.cv; @@ -271,6 +287,12 @@ void burn_to_eos (const BurnT& burn_state, T& eos_state) } #endif +#ifdef PHOTOCHEMISTRY + for (int n = 0; n < NumRadEqs; ++n) { + eos_state.rn[n] = burn_state.rn[n]; + } +#endif + eos_state.cv = burn_state.cv; eos_state.y_e = burn_state.y_e; diff --git a/interfaces/eos.H b/interfaces/eos.H index 3abccc2c77..d08120644f 100644 --- a/interfaces/eos.H +++ b/interfaces/eos.H @@ -1,6 +1,7 @@ #ifndef EOS_H #define EOS_H +#include #include #include #include @@ -72,8 +73,9 @@ void eos_reset (T& state, bool& has_been_reset) { state.T = amrex::min(EOSData::maxtemp, amrex::max(EOSData::mintemp, state.T)); state.rho = amrex::min(EOSData::maxdens, amrex::max(EOSData::mindens, state.rho)); - + // std::cout << "eos_reset()--> Calling actual_eos" << std::endl; actual_eos(eos_input_rt, state); + // std::cout << "eos_reset()--> Finished actual_eos" << std::endl; has_been_reset = true; } @@ -424,9 +426,11 @@ void eos (const I input, T& state) eos_override(state); // Call the EOS. - + // std::cout << "state.xn[1] in eos.h: " << state.xn[1] << std::endl; if (!has_been_reset) { + // std::cout << "eos()--> Calling actual_eos" << std::endl; actual_eos(input, state); + // std::cout << "eos()--> Finished actual_eos" << std::endl; } } From 8d4fd1c771dc5248270d2a63b24241b5559de095 Mon Sep 17 00:00:00 2001 From: James471 Date: Mon, 2 Mar 2026 23:36:54 +1100 Subject: [PATCH 3/8] WIP: Add photoionization EOS --- CMakeLists.txt | 102 +++++++++- EOS/photoionization/Make.package | 3 + EOS/photoionization/_parameters | 16 ++ EOS/photoionization/actual_eos.H | 236 ++++++++++++++++++++++++ EOS/photoionization/actual_eos_data.H | 11 ++ EOS/photoionization/actual_eos_data.cpp | 4 + 6 files changed, 371 insertions(+), 1 deletion(-) create mode 100644 EOS/photoionization/Make.package create mode 100644 EOS/photoionization/_parameters create mode 100644 EOS/photoionization/actual_eos.H create mode 100644 EOS/photoionization/actual_eos_data.H create mode 100644 EOS/photoionization/actual_eos_data.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f4eb5cc45b..3a58aa10e6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20) project(Microphysics VERSION 1.0.0 - DESCRIPTION "building primordial_chemistry or metal_chemistry networks in Microphysics with CMake" + DESCRIPTION "building primordial_chemistry or metal_chemistry or photoionization networks in Microphysics with CMake" LANGUAGES CXX C) #---------------------------------------------------------------------------------------------------------------------- @@ -93,6 +93,106 @@ function(setup_target_for_microphysics_compilation network_name output_dir) execute_process(COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${PYTHONPATH}:${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null" python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/get_naux.py" --net "${networkdir}" --microphysics_path "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/" WORKING_DIRECTORY ${output_dir}/) execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/write_network.py" --header_template "${networkheadertemplatefile}" --header_output "${networkpropfile}" -s "${networkfile}" WORKING_DIRECTORY ${output_dir}/) + elseif (${network_name} STREQUAL "H_ionization_TInd") + #need these to write extern_parameters.H + set(paramfile "${CMAKE_CURRENT_LIST_DIR}/_parameters") + # TODO(#1)(james471): See how to do eosparamfile correctly + set(EOSparamfile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization/_parameters") + set(networkpcparamfile "${CMAKE_SOURCE_DIR}/src/networks/H_ionization_TInd/_parameters") + + #similarly, we want network_properties.H + set(networkpropfile "${output_dir}/network_properties.H") + # set(networkfile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/primordial_chem/pynucastro.net") + set(networkdir "${CMAKE_SOURCE_DIR}/src/networks/H_ionization_TInd/") + set(networkheadertemplatefile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/network_header.template") + + #DO NOT change the order of the directories below! + set (photoionization_dirs ${CMAKE_BINARY_DIR} ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/gcem/include + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/VODE ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/utils + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization + ${CMAKE_SOURCE_DIR}/src/networks/H_ionization_TInd ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/constants + PARENT_SCOPE) + + #we need to have extern_parameters.cpp be available at configure time + #the script write_probin.py writes this .cpp file so we call it here + #note, execute_process only works on 'cmake' and not 'make' + #so, if any of the _parameter files are changed, one needs to re-run 'cmake' + #to generate updated header files + + if(BUILD_UNIT_TEST_PC) + message(STATUS "In BUILD_UNIT_TEST_PC") + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${paramfile} ${EOSparamfile} + ${networkpcparamfile} ${networkparamfile} ${VODEparamfile} ${integrationparamfile} ${unittestparamfile}" --use_namespace WORKING_DIRECTORY ${output_dir}/) + else() + message(STATUS "Not in BUILD_UNIT_TEST_PC") + #do not need paramfile and unittestparamfile + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${EOSparamfile} ${networkpcparamfile} + ${networkparamfile} ${VODEparamfile} ${integrationparamfile} " --use_namespace WORKING_DIRECTORY ${output_dir}/) + endif() + + set(photoionization_sources ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces/eos_data.cpp + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces/network_initialization.cpp + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization/actual_eos_data.cpp + ${CMAKE_SOURCE_DIR}/src/networks/H_ionization_TInd/actual_network_data.cpp + ${output_dir}/extern_parameters.cpp PARENT_SCOPE) + + + #below for NAUX + execute_process(COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${PYTHONPATH}:${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null" python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/get_naux.py" --net "${networkdir}" --microphysics_path "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/" WORKING_DIRECTORY ${output_dir}/) + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/write_network.py" --header_template "${networkheadertemplatefile}" --header_output "${networkpropfile}" -s "${networkfile}" WORKING_DIRECTORY ${output_dir}/) + + elseif (${network_name} STREQUAL "photoionization") + #need these to write extern_parameters.H + set(paramfile "${CMAKE_CURRENT_LIST_DIR}/_parameters") + # TODO(#1)(james471): See how to do eosparamfile correctly + set(EOSparamfile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization/_parameters") + set(networkpcparamfile "${CMAKE_SOURCE_DIR}/src/networks/photoionization/_parameters") + + #similarly, we want network_properties.H + set(networkpropfile "${output_dir}/network_properties.H") + # set(networkfile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/primordial_chem/pynucastro.net") + set(networkdir "${CMAKE_SOURCE_DIR}/src/networks/photoionization/") + set(networkheadertemplatefile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/network_header.template") + + #DO NOT change the order of the directories below! + set (photoionization_dirs ${CMAKE_BINARY_DIR} ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/gcem/include + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/VODE ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/utils + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization + ${CMAKE_SOURCE_DIR}/src/networks/photoionization ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/constants + PARENT_SCOPE) + + #we need to have extern_parameters.cpp be available at configure time + #the script write_probin.py writes this .cpp file so we call it here + #note, execute_process only works on 'cmake' and not 'make' + #so, if any of the _parameter files are changed, one needs to re-run 'cmake' + #to generate updated header files + + if(BUILD_UNIT_TEST_PC) + message(STATUS "In BUILD_UNIT_TEST_PC") + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${paramfile} ${EOSparamfile} + ${networkpcparamfile} ${networkparamfile} ${VODEparamfile} ${integrationparamfile} ${unittestparamfile}" --use_namespace WORKING_DIRECTORY ${output_dir}/) + else() + message(STATUS "Not in BUILD_UNIT_TEST_PC") + #do not need paramfile and unittestparamfile + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${EOSparamfile} ${networkpcparamfile} + ${networkparamfile} ${VODEparamfile} ${integrationparamfile} " --use_namespace WORKING_DIRECTORY ${output_dir}/) + endif() + + set(photoionization_sources ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces/eos_data.cpp + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces/network_initialization.cpp + ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization/actual_eos_data.cpp + ${CMAKE_SOURCE_DIR}/src/networks/photoionization/actual_network_data.cpp + ${output_dir}/extern_parameters.cpp PARENT_SCOPE) + + + #below for NAUX + execute_process(COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${PYTHONPATH}:${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null" python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/get_naux.py" --net "${networkdir}" --microphysics_path "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/" WORKING_DIRECTORY ${output_dir}/) + execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/write_network.py" --header_template "${networkheadertemplatefile}" --header_output "${networkpropfile}" -s "${networkfile}" WORKING_DIRECTORY ${output_dir}/) + elseif (${network_name} STREQUAL "metal_chem") #need these to write extern_parameters.H set(paramfile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/unit_test/burn_cell_metal_chem/_parameters") diff --git a/EOS/photoionization/Make.package b/EOS/photoionization/Make.package new file mode 100644 index 0000000000..c27842afcb --- /dev/null +++ b/EOS/photoionization/Make.package @@ -0,0 +1,3 @@ +CEXE_headers += actual_eos.H +CEXE_headers += actual_eos_data.H +CEXE_sources += actual_eos_data.cpp diff --git a/EOS/photoionization/_parameters b/EOS/photoionization/_parameters new file mode 100644 index 0000000000..c2eaca28fd --- /dev/null +++ b/EOS/photoionization/_parameters @@ -0,0 +1,16 @@ +@namespace: eos + +eos_gamma_default real 1.4 + +# define the specie names, and their masses and gammas +species_1_name string "Electron" +species_1_gamma real 5./3. +species_1_mass real 9.10938188e-28 + +species_2_name string "Hydrogen" +species_2_gamma real 5./3. +species_2_mass real 1.67353251819e-24 + +species_3_name string "Hydrogen_Ion" +species_3_gamma real 5./3. +species_3_mass real 1.67353251819e-24 diff --git a/EOS/photoionization/actual_eos.H b/EOS/photoionization/actual_eos.H new file mode 100644 index 0000000000..1254b83859 --- /dev/null +++ b/EOS/photoionization/actual_eos.H @@ -0,0 +1,236 @@ +#ifndef ACTUAL_EOS_H +#define ACTUAL_EOS_H + +// This is a mult-gamma EOS for photochemistry. +// Each species can have its own gamma, but +// otherwise, they all act as an ideal gas. + +#include "eos_type.H" +#include +#include +#include +#include +#include +#include +#include +#include + +const std::string eos_name = "multigamma"; + +inline +void actual_eos_init () +{ + + // Set the gammas & masses for the species + + for (int n = 0; n < NumSpec; ++n) { + gammas[n] = eos_rp::eos_gamma_default; + spmasses[n] = 1.67353251819e-24; + } + + int idx; + // Set the gammas & masses for the species + #define GET_SPECIES_PARAMS(num) do { \ + idx = network_spec_index(eos_rp::species_##num##_name); \ + if (idx >= 0) { \ + gammas[idx] = eos_rp::species_##num##_gamma; \ + spmasses[idx] = eos_rp::species_##num##_mass; \ + } \ + } while (0) + + GET_SPECIES_PARAMS(1); + GET_SPECIES_PARAMS(2); + GET_SPECIES_PARAMS(3); + + #undef GET_SPECIES_PARAMS +} + + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +bool is_input_valid (I input) +{ + static_assert(std::is_same_v, "input must be either eos_input_rt or eos_input_re"); + + bool valid = false; + + if (input == eos_input_rt || + input == eos_input_re) { + valid = true; + } + + return valid; +} + + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void actual_eos (I input, T& state) +{ + static_assert(std::is_same_v, "input must be either eos_input_rt or eos_input_re"); + const amrex::Real gasconstant = C::n_A * C::k_B; + const amrex::Real protonmass = C::m_p; + // Special gamma factors + amrex::Real sum_Abarinv = 0.0_rt; + amrex::Real sum_gammasinv = 0.0_rt; + amrex::Real rhotot = 0.0_rt; + + for (int n = 0; n < NumSpec; ++n) { + rhotot += state.xn[n]*spmasses[n]; + // std::cout << "actual_eos()--> xn[" << n << "]: " << state.xn[n] << ", spmasses[" << n << "]: " << spmasses[n] << std::endl; + } + // std::cout << "Density in EOS: " << rhotot << std::endl; + + for (int n = 0; n < NumSpec; ++n) { + sum_Abarinv += state.xn[n]; + sum_gammasinv += (state.xn[n]*protonmass/rhotot) * (1.0/(gammas[n]-1.0)); + } + // std::cout << "actual_eos()--> rhotot: " << rhotot << ", state.xn[0]: " << state.xn[0] << ", state.xn[1]: " << state.xn[1] << ", state.xn[2]: " << state.xn[2] << std::endl; + sum_Abarinv *= protonmass/rhotot; + state.mu = 1.0 / sum_Abarinv; + sum_gammasinv /= sum_Abarinv; + + //------------------------------------------------------------------------- + // For all EOS input modes EXCEPT eos_input_rt, first compute dens + // and temp as needed from the inputs. + //------------------------------------------------------------------------- + + amrex::Real temp = NAN; + amrex::Real dens = NAN; + amrex::Real eint = NAN; + // std::cout << "Eos was called" << std::endl; + switch (input) { + + case eos_input_rt: + + // dens, temp and xmass are inputs + // We don't need to do anything here + temp = state.T; + dens = state.rho; + eint = sum_gammasinv * sum_Abarinv * gasconstant * state.T; + // std::cout << "actual_eos() --> eos_input_rt: temp: " << temp << ", dens: " << dens << ", eint: " << eint << ", sum_gammasinv: " << sum_gammasinv << ", sum_Abarinv: " << sum_Abarinv << std::endl; + break; + + case eos_input_re: + + // dens, energy, and xmass are inputs + // Solve for the temperature + + if constexpr (has_energy::value) { + dens = state.rho; + + // stop the integration if the internal energy < 0 + AMREX_ASSERT(state.e > 0.); + temp = state.e /( sum_gammasinv * gasconstant * sum_Abarinv); + eint = state.e; + // std::cout << "actual_eos() --> eos_input_re: temp: " << temp << ", dens: " << dens << ", eint: " << eint << ", sum_gammasinv: " << sum_gammasinv << ", sum_Abarinv: " << sum_Abarinv << std::endl; + } + + break; + + case eos_input_rh: + // dens, enthalpy, and xmass are inputs +#ifndef AMREX_USE_GPU + amrex::Error("eos_input_rh is not supported"); +#endif + + break; + + case eos_input_tp: + // temp, pressure, and xmass are inputs +#ifndef AMREX_USE_GPU + amrex::Error("eos_input_tp is not supported"); +#endif + + break; + + case eos_input_rp: + // dens, pressure, and xmass are inputs + if constexpr (has_pressure::value) { + dens = state.rho; + + // stop the integration if the pressure < 0 + AMREX_ASSERT(state.p > 0.); + eint = state.p * sum_gammasinv / dens; + temp = eint / (sum_gammasinv * gasconstant * sum_Abarinv); + } + break; + + case eos_input_ps: + // pressure, entropy, and xmass are inputs +#ifndef AMREX_USE_GPU + amrex::Error("eos_input_ps is not supported"); +#endif + + break; + + case eos_input_ph: + // pressure, enthalpy, and xmass are inputs +#ifndef AMREX_USE_GPU + amrex::Error("eos_input_ph is not supported"); +#endif + + break; + + case eos_input_th: + // temp, enthalpy, and xmass are inputs +#ifndef AMREX_USE_GPU + amrex::Error("eos_input_th is not supported"); +#endif + + break; + + default: + +#ifndef AMREX_USE_GPU + amrex::Error("EOS: invalid input."); +#endif + + break; + + } + + //------------------------------------------------------------------------- + // Now we have the density and temperature (and mass fractions / + // mu), regardless of the inputs. + //------------------------------------------------------------------------- + + state.T = temp; + state.rho = dens; + + if constexpr (has_energy::value) { + state.e = eint; + } + + if constexpr (has_pressure::value) { + amrex::Real pressure = state.rho * eint / sum_gammasinv; + state.p = pressure; + + state.dpdT = pressure / temp; + state.dpdr = pressure / dens; + state.cs = std::sqrt((1.0 + 1.0/sum_gammasinv) * state.p /state.rho); + if constexpr (has_G::value) { + state.G = 0.5 * (1.0 + (1.0 + 1.0/sum_gammasinv)); + } + } + + amrex::Real dedT = sum_gammasinv * sum_Abarinv * gasconstant; + amrex::Real dedr = 0.0_rt; + if constexpr (has_energy::value) { + state.dedT = dedT; + state.dedr = dedr; + if constexpr (has_pressure::value) { + state.dpde = state.dpdT / state.dedT; + } + } + +} + + + +inline +void actual_eos_finalize () +{ +} + +#endif diff --git a/EOS/photoionization/actual_eos_data.H b/EOS/photoionization/actual_eos_data.H new file mode 100644 index 0000000000..5a06b0cd44 --- /dev/null +++ b/EOS/photoionization/actual_eos_data.H @@ -0,0 +1,11 @@ +#ifndef actual_eos_data_H +#define actual_eos_data_H + +#include +#include +#include + +extern AMREX_GPU_MANAGED amrex::Real gammas[NumSpec]; +extern AMREX_GPU_MANAGED amrex::Real spmasses[NumSpec]; + +#endif diff --git a/EOS/photoionization/actual_eos_data.cpp b/EOS/photoionization/actual_eos_data.cpp new file mode 100644 index 0000000000..e2559d5ec0 --- /dev/null +++ b/EOS/photoionization/actual_eos_data.cpp @@ -0,0 +1,4 @@ +#include + +AMREX_GPU_MANAGED amrex::Real gammas[NumSpec]; +AMREX_GPU_MANAGED amrex::Real spmasses[NumSpec]; From 2f8765ad1774450ca48913191c6894b59112f21c Mon Sep 17 00:00:00 2001 From: James471 Date: Fri, 6 Mar 2026 16:37:35 +1100 Subject: [PATCH 4/8] WIP: Add reduced speed of light to burn state and some missing radiation variables WIP: Add reduced speed of light to burn state and some missing radiation variables WIP: Add reduced speed of light to burn state and some missing radiation variables --- integration/VODE/vode_dvhin.H | 2 +- integration/VODE/vode_dvstep.H | 15 ++------------- integration/VODE/vode_type.H | 3 +++ integration/_parameters | 2 +- integration/integrator_setup_strang.H | 7 +++++++ integration/integrator_type.H | 2 +- integration/integrator_type_strang.H | 7 +++++++ integration/utils/nonaka_plot.H | 5 +++++ integration/utils/numerical_jacobian.H | 12 ++++++++++++ interfaces/burn_type.H | 3 +++ interfaces/eos.H | 7 +++++++ interfaces/eos_type.H | 20 ++++++++++++++++++++ networks/rhs.H | 2 ++ 13 files changed, 71 insertions(+), 16 deletions(-) diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H index 180fa267ef..8e893b8f47 100644 --- a/integration/VODE/vode_dvhin.H +++ b/integration/VODE/vode_dvhin.H @@ -45,7 +45,7 @@ void dvhin (BurnT& state, DvodeT& vstate, amrex::Real& H0, int& NITER, int& IER) amrex::Real HUB = PT1 * TDIST; for (int i = 1; i <= int_neqs; ++i) { - // TODO(#3)(james471): This can be optimized to reduce branching + // TODO(#3)(james471): This can probably be optimized to reduce branching amrex::Real atol{}; #if defined (PHOTOCHEMISTRY) if (i <= NumSpec) { diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index 964098955f..6a02c48c2f 100644 --- a/integration/VODE/vode_dvstep.H +++ b/integration/VODE/vode_dvstep.H @@ -1,6 +1,7 @@ #ifndef VODE_DVSTEP_H #define VODE_DVSTEP_H +#include "network.H" #include #include #include @@ -148,7 +149,6 @@ int dvstep (BurnT& state, DvodeT& vstate) } #endif #ifdef SDC -// TODO(#3)(james471): Check what needs to be done for radiation variables amrex::Real rho_old = state.rho_orig + TOLD * state.ydot_a[SRHO]; for (int i = 1; i <= int_neqs; ++i) { y_save(i) = vstate.y(i); @@ -235,7 +235,6 @@ int dvstep (BurnT& state, DvodeT& vstate) } #endif - // TODO(#2)(james471): Add radiation checks here for (int i = 1; i <= NumSpec; ++i) { #ifdef STRANG @@ -310,7 +309,7 @@ int dvstep (BurnT& state, DvodeT& vstate) } - for (int i = net_ienuc + 1; i <= int_neqs; i+=2) { + for (int i = net_ienuc + 1; i <= int_neqs; i+=MicrophysicsNumRadVarsPerGroup) { #ifdef STRANG // Constrain radiation quantities such that they don't change by more than a certain @@ -349,16 +348,6 @@ int dvstep (BurnT& state, DvodeT& vstate) valid_update = false; break; } - - // Don't enforce the condition below if - // vstate.y contains number densities - if (! integrator_rp::use_number_densities) { - if (vstate.y(i) > 1.0_rt + radiation_failure_tolerance) { - valid_update = false; - break; - } - } - #endif #ifdef SDC diff --git a/integration/VODE/vode_type.H b/integration/VODE/vode_type.H index 8f389cfa36..8152573492 100644 --- a/integration/VODE/vode_type.H +++ b/integration/VODE/vode_type.H @@ -143,6 +143,9 @@ struct dvode_t // Local time and integration end time amrex::Real t, tout; + // Reduced speed of light + amrex::Real c_hat; + // Integration array amrex::Array1D y; diff --git a/integration/_parameters b/integration/_parameters index 9f2ccb04d6..ac332727fc 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -41,7 +41,7 @@ rtol_rad_num real 1.e-12 rtol_rad_flux real 1.e-6 # absolute tolerance for radiation density -atol_rad_num real 1.e-8 +atol_rad_num real 1.e-12 # absolute tolerance for radiation flux atol_rad_flux real 1.e-6 diff --git a/integration/integrator_setup_strang.H b/integration/integrator_setup_strang.H index bd97c7ae23..b24fec90eb 100644 --- a/integration/integrator_setup_strang.H +++ b/integration/integrator_setup_strang.H @@ -20,7 +20,10 @@ struct state_backup_t { amrex::Real e_in{}; #ifndef AMREX_USE_GPU amrex::Real xn_in[NumSpec]{}; +#ifdef PHOTOCHEMISTRY amrex::Real rn_in[NumRadEqs]{}; + amrex::Real c_hat; +#endif #endif }; @@ -125,9 +128,12 @@ state_backup_t integrator_backup (const BurnT& state) { for (int n = 0; n < NumSpec; ++n) { state_save.xn_in[n] = state.xn[n]; } +#ifdef PHOTOCHEMISTRY for (int n = 0; n < NumRadEqs; ++n) { state_save.rn_in[n] = state.rn[n]; } + state_save.c_hat = state.c_hat; +#endif #endif state_save.T_in = state.T; state_save.e_in = state.e; @@ -238,6 +244,7 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state, } #ifdef PHOTOCHEMISTRY std::cout << std::endl; + std::cout << "(Reduced) Speed of light = " << std::setprecision(16) << state.c_hat << std::endl; std::cout << "rn start = "; for (const double X : state_save.rn_in) { std::cout << std::scientific << std::setprecision(16) << X << " "; diff --git a/integration/integrator_type.H b/integration/integrator_type.H index fe3d4081b5..de0c83700c 100644 --- a/integration/integrator_type.H +++ b/integration/integrator_type.H @@ -34,7 +34,6 @@ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void update_thermodynamics (BurnT& state, const I& int_state) { - // TODO(#3)(james471): Add radiation check for aux thermo #ifdef AUX_THERMO // we are handling the thermodynamics via the aux quantities, which // are stored in the burn_t here, so we need to update those based on @@ -63,6 +62,7 @@ void update_thermodynamics (BurnT& state, const I& int_state) for (int n = 1; n <= NumRadEqs; ++n) { state.rn[n-1] = int_state.y(net_ienuc + n); } + state.c_hat = int_state.c_hat; #endif if (integrator_rp::scale_system) { diff --git a/integration/integrator_type_strang.H b/integration/integrator_type_strang.H index 8fb8f15d35..63396c7f87 100644 --- a/integration/integrator_type_strang.H +++ b/integration/integrator_type_strang.H @@ -22,6 +22,7 @@ void integrator_to_burn (const I& int_state, BurnT& state) for (int n = 1; n <= NumRadEqs; ++n) { state.rn[n-1] = int_state.y(net_ienuc + n); } + state.c_hat = int_state.c_hat; #endif if (integrator_rp::scale_system) { @@ -48,6 +49,7 @@ void burn_to_integrator (const BurnT& state, I& int_state) for (int n = 1; n <= NumRadEqs; ++n) { int_state.y(net_ienuc + n) = state.rn[n-1]; } + int_state.c_hat = state.c_hat; #endif if (integrator_rp::scale_system) { @@ -77,6 +79,11 @@ void clean_state (const amrex::Real time, BurnT& state, I& int_state) int_state.y(n) = amrex::max(int_state.y(n), integrator_rp::SMALL_X_SAFE); } } +#ifdef PHOTOCHEMISTRY + for (int n = net_ienuc+1; n <= net_ienuc+NumRadEqs; n += MicrophysicsNumRadVarsPerGroup) { + int_state.y(n) = amrex::max(int_state.y(n), integrator_rp::SMALL_X_SAFE); + } +#endif // Renormalize the abundances as necessary. diff --git a/integration/utils/nonaka_plot.H b/integration/utils/nonaka_plot.H index dd70eeb23f..0fa5d42d1e 100644 --- a/integration/utils/nonaka_plot.H +++ b/integration/utils/nonaka_plot.H @@ -96,6 +96,11 @@ void nonaka_rhs(const amrex::Real time, const BurnT& state, const YdotNetArray1D for (int n = 0; n < NumSpec; ++n) { nf << std::setw(FIELD_WIDTH) << state.xn[n] << " "; } +#ifdef PHOTOCHEMISTRY + for (int n = 0; n < NumRadEqs; ++n) { + nf << std::setw(FIELD_WIDTH) << state.rn[n] << " "; + } +#endif // Convert molar fraction rhs to mass fraction rhs dX/dt diff --git a/integration/utils/numerical_jacobian.H b/integration/utils/numerical_jacobian.H index d396c69bcd..598566d26c 100644 --- a/integration/utils/numerical_jacobian.H +++ b/integration/utils/numerical_jacobian.H @@ -29,6 +29,8 @@ /// difference in terms of X and T and then convert the Jacobian /// elements to be in terms of X and e +// TODO(#2)(james471): This needs fixes for radiation + struct jac_info_t { amrex::Real h; }; @@ -76,6 +78,16 @@ void numerical_jac(BurnT& state, const jac_info_t& jac_info, JacNetArray2D& jac) w = 1.0_rt / (integrator_rp::rtol_spec * std::abs(state.xn[i-1]) + integrator_rp::atol_spec); fac += (ydotm(i) * w) * (ydotm(i) * w); } +#ifdef PHOTOCHEMISTRY + for (int i = net_ienuc + 1; i <= NumRadEqs; i += MicrophysicsNumRadVarsPerGroup) { + w = 1.0_rt / (integrator_rp::rtol_rad_num * std::abs(state.rn[i-1]) + integrator_rp::atol_rad_num); + fac += (ydotm(i) * w) * (ydotm(i) * w); + for (int j = i + 1; j <= i + MicrophysicsNumRadVarsPerGroup; j++) { + w = 1.0_rt / (integrator_rp::rtol_rad_flux * std::abs(state.rn[j-1]) + integrator_rp::atol_rad_flux); + fac += (ydotm(j) * w) * (ydotm(j) * w); + } + } +#endif // even though we difference wrt T below, for |f|, we want to use // energy to get the units correct diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H index 8665548147..0cf5079ac0 100644 --- a/interfaces/burn_type.H +++ b/interfaces/burn_type.H @@ -86,6 +86,7 @@ struct burn_t #endif #ifdef PHOTOCHEMISTRY amrex::Real rn[NumRadEqs]{}; + amrex::Real c_hat; #endif // scaling used to make e we integrate dimensionless @@ -248,6 +249,7 @@ void eos_to_burn (const T& eos_state, BurnT& burn_state) for (int n = 0; n < NumRadEqs; ++n) { burn_state.rn[n] = eos_state.rn[n]; } + burn_state.c_hat = eos_state.c_hat; #endif burn_state.cv = eos_state.cv; @@ -291,6 +293,7 @@ void burn_to_eos (const BurnT& burn_state, T& eos_state) for (int n = 0; n < NumRadEqs; ++n) { eos_state.rn[n] = burn_state.rn[n]; } + eos_state.c_hat = burn_state.c_hat; #endif eos_state.cv = burn_state.cv; diff --git a/interfaces/eos.H b/interfaces/eos.H index d08120644f..0b5e75eacb 100644 --- a/interfaces/eos.H +++ b/interfaces/eos.H @@ -218,6 +218,13 @@ void print_state (T& state) { std::cout << state.xn[n] << " "; } std::cout << std::endl; +#ifdef PHOTOCHEMISTRY + std::cout << "RAD = "; + for (int n = 0; n < NumRadEqs; n ++) { + std::cout << state.rn[n] << " "; + } + std::cout << std::endl; +#endif std::cout << "Y_E = " << state.y_e; } diff --git a/interfaces/eos_type.H b/interfaces/eos_type.H index 7db2db24a6..ebf844264e 100644 --- a/interfaces/eos_type.H +++ b/interfaces/eos_type.H @@ -14,6 +14,10 @@ struct eos_t:eos_base_t { amrex::Real h{}; amrex::Real s{}; amrex::Real xn[NumSpec]{}; +#ifdef PHOTOCHEMISTRY + amrex::Real rn[NumRadEqs]{}; + amrex::Real c_hat{}; +#endif #if NAUX_NET > 0 amrex::Real aux[NumAux]{}; #endif @@ -58,6 +62,10 @@ struct eos_re_t:eos_base_t { amrex::Real T{}; amrex::Real e{}; amrex::Real xn[NumSpec]{}; +#ifdef PHOTOCHEMISTRY + amrex::Real rn[NumRadEqs]{}; + amrex::Real c_hat{}; +#endif #if NAUX_NET > 0 amrex::Real aux[NumAux]{}; #endif @@ -88,6 +96,10 @@ struct eos_rep_t:eos_base_t { amrex::Real xn[NumSpec]{}; #if NAUX_NET > 0 amrex::Real aux[NumAux]{}; +#endif +#ifdef PHOTOCHEMISTRY + amrex::Real rn[NumRadEqs]{}; + amrex::Real c_hat{}; #endif amrex::Real dedT{}; amrex::Real dedr{}; @@ -302,6 +314,14 @@ std::ostream& print_state (std::ostream& o, T const& eos_state) o << eos_state.aux[n] << " "; } o << std::endl; +#endif +#ifdef PHOTOCHEMISTRY + o << "rn = "; + for (int n = 0; n < NumRadEqs; ++n) { + o << eos_state.rn[n] << " "; + } + o << std::endl; + o << "c_hat = " << eos_state.c_hat << std::endl; #endif } if constexpr (has_pressure::value) { diff --git a/networks/rhs.H b/networks/rhs.H index 5d57aa5ac5..c2780a49fb 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -20,6 +20,8 @@ #include #include +// TODO(#2)(james471): Check if this requires changes for radiation variables + #ifdef NEW_NETWORK_IMPLEMENTATION // Forward declarations From 4783e8ebd6ef64ff1d7ced99a8c036218beb0029 Mon Sep 17 00:00:00 2001 From: James471 Date: Fri, 10 Apr 2026 20:29:26 +1000 Subject: [PATCH 5/8] WIP: Minor fixes WIP: Minor fixes WIP: Minor fixes --- integration/VODE/vode_dvstep.H | 2 ++ integration/_parameters | 22 +++++++++++----------- integration/integrator_data.H | 2 ++ integration/integrator_setup_strang.H | 2 ++ 4 files changed, 17 insertions(+), 11 deletions(-) diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index 6a02c48c2f..7f00ae0847 100644 --- a/integration/VODE/vode_dvstep.H +++ b/integration/VODE/vode_dvstep.H @@ -309,6 +309,7 @@ int dvstep (BurnT& state, DvodeT& vstate) } +#ifdef PHOTOCHEMISTRY for (int i = net_ienuc + 1; i <= int_neqs; i+=MicrophysicsNumRadVarsPerGroup) { #ifdef STRANG @@ -349,6 +350,7 @@ int dvstep (BurnT& state, DvodeT& vstate) break; } #endif +#endif #ifdef SDC // Abort the program. This is unimplemented diff --git a/integration/_parameters b/integration/_parameters index ac332727fc..3b79841575 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -23,28 +23,28 @@ burner_verbose bool 0 # species and energy equations. # relative tolerance for species -rtol_spec real 1.e-12 +rtol_spec real 1.e-10 # relative tolerance for energy -rtol_enuc real 1.e-6 +rtol_enuc real 1.e-10 # absolute tolerance for species -atol_spec real 1.e-8 +atol_spec real 1.e-10 # absolute tolerance for energy -atol_enuc real 1.e-6 +atol_enuc real 1.e-25 -# relative tolerance for radiation density -rtol_rad_num real 1.e-12 +# relative tolerance for radiation number density +rtol_rad_num real 1.e-10 -# relative teolerance for radiation flux -rtol_rad_flux real 1.e-6 +# relative tolerance for radiation flux +rtol_rad_flux real 1.e-10 -# absolute tolerance for radiation density -atol_rad_num real 1.e-12 +# absolute tolerance for radiation number density +atol_rad_num real 1.e-10 # absolute tolerance for radiation flux -atol_rad_flux real 1.e-6 +atol_rad_flux real 1.e-10 # Whether to renormalize the mass fractions at each step in the evolution # so that they sum to unity. diff --git a/integration/integrator_data.H b/integration/integrator_data.H index 28950af3ec..75d070af6c 100644 --- a/integration/integrator_data.H +++ b/integration/integrator_data.H @@ -15,10 +15,12 @@ constexpr int INT_NEQS = NumRadEqs + NumSpec + 1; // abundance is unreasonably small or large (each X must satisfy // -failure_tolerance <= X <= 1.0 + failure_tolerance). constexpr amrex::Real species_failure_tolerance = 1.e-2_rt; +#ifdef PHOTOCHEMISTRY // We will use this parameter to determine if // a given radiation quantity (number density or flux) is // unreasonably small (each quantity must satisfy X >= -radiation_failure_tolerance). constexpr amrex::Real radiation_failure_tolerance = 1.e-2_rt; +#endif #ifdef NSE constexpr int MIN_NSE_BAILOUT_STEPS = 10; diff --git a/integration/integrator_setup_strang.H b/integration/integrator_setup_strang.H index b24fec90eb..1ef864c241 100644 --- a/integration/integrator_setup_strang.H +++ b/integration/integrator_setup_strang.H @@ -207,11 +207,13 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state, } } } +#ifdef PHOTOCHEMISTRY for (int n = net_ienuc+1; n <= net_ienuc+NumRadEqs; n += MicrophysicsNumRadVarsPerGroup) { if (int_state.y(n) < -radiation_failure_tolerance) { state.success = false; } } +#endif #ifndef AMREX_USE_GPU if (integrator_rp::burner_verbose) { From a01aae94caf00c0bec0979af881542c7c795da90 Mon Sep 17 00:00:00 2001 From: James471 Date: Wed, 15 Apr 2026 01:28:06 +1000 Subject: [PATCH 6/8] Cleanup --- interfaces/eos.H | 5 ----- 1 file changed, 5 deletions(-) diff --git a/interfaces/eos.H b/interfaces/eos.H index 0b5e75eacb..4de946910e 100644 --- a/interfaces/eos.H +++ b/interfaces/eos.H @@ -73,9 +73,7 @@ void eos_reset (T& state, bool& has_been_reset) { state.T = amrex::min(EOSData::maxtemp, amrex::max(EOSData::mintemp, state.T)); state.rho = amrex::min(EOSData::maxdens, amrex::max(EOSData::mindens, state.rho)); - // std::cout << "eos_reset()--> Calling actual_eos" << std::endl; actual_eos(eos_input_rt, state); - // std::cout << "eos_reset()--> Finished actual_eos" << std::endl; has_been_reset = true; } @@ -433,11 +431,8 @@ void eos (const I input, T& state) eos_override(state); // Call the EOS. - // std::cout << "state.xn[1] in eos.h: " << state.xn[1] << std::endl; if (!has_been_reset) { - // std::cout << "eos()--> Calling actual_eos" << std::endl; actual_eos(input, state); - // std::cout << "eos()--> Finished actual_eos" << std::endl; } } From cc5d13720803a36c170f202f743481c67c265fc6 Mon Sep 17 00:00:00 2001 From: James471 Date: Wed, 15 Apr 2026 01:29:58 +1000 Subject: [PATCH 7/8] Simplify photoionization eos implementation, change species masses Simplify photoionization eos implementation, change species masses Simplify photoionization eos implementation, change species masses --- CMakeLists.txt | 53 -------------------------- EOS/photoionization/_parameters | 6 +-- EOS/photoionization/actual_eos.H | 40 +++++++------------ integration/VODE/vode_dvhin.H | 1 - integration/VODE/vode_dvstep.H | 2 +- integration/utils/numerical_jacobian.H | 2 +- interfaces/burn_type.H | 7 +++- networks/rhs.H | 2 +- 8 files changed, 26 insertions(+), 87 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3a58aa10e6..a1be9c7455 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -88,57 +88,6 @@ function(setup_target_for_microphysics_compilation network_name output_dir) ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/primordial_chem/actual_network_data.cpp ${output_dir}/extern_parameters.cpp PARENT_SCOPE) - - #below for NAUX - execute_process(COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${PYTHONPATH}:${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null" python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/get_naux.py" --net "${networkdir}" --microphysics_path "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/" WORKING_DIRECTORY ${output_dir}/) - execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/write_network.py" --header_template "${networkheadertemplatefile}" --header_output "${networkpropfile}" -s "${networkfile}" WORKING_DIRECTORY ${output_dir}/) - - elseif (${network_name} STREQUAL "H_ionization_TInd") - #need these to write extern_parameters.H - set(paramfile "${CMAKE_CURRENT_LIST_DIR}/_parameters") - # TODO(#1)(james471): See how to do eosparamfile correctly - set(EOSparamfile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization/_parameters") - set(networkpcparamfile "${CMAKE_SOURCE_DIR}/src/networks/H_ionization_TInd/_parameters") - - #similarly, we want network_properties.H - set(networkpropfile "${output_dir}/network_properties.H") - # set(networkfile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/primordial_chem/pynucastro.net") - set(networkdir "${CMAKE_SOURCE_DIR}/src/networks/H_ionization_TInd/") - set(networkheadertemplatefile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/network_header.template") - - #DO NOT change the order of the directories below! - set (photoionization_dirs ${CMAKE_BINARY_DIR} ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/gcem/include - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/VODE ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration/utils - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/integration ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization - ${CMAKE_SOURCE_DIR}/src/networks/H_ionization_TInd ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/constants - PARENT_SCOPE) - - #we need to have extern_parameters.cpp be available at configure time - #the script write_probin.py writes this .cpp file so we call it here - #note, execute_process only works on 'cmake' and not 'make' - #so, if any of the _parameter files are changed, one needs to re-run 'cmake' - #to generate updated header files - - if(BUILD_UNIT_TEST_PC) - message(STATUS "In BUILD_UNIT_TEST_PC") - execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${paramfile} ${EOSparamfile} - ${networkpcparamfile} ${networkparamfile} ${VODEparamfile} ${integrationparamfile} ${unittestparamfile}" --use_namespace WORKING_DIRECTORY ${output_dir}/) - else() - message(STATUS "Not in BUILD_UNIT_TEST_PC") - #do not need paramfile and unittestparamfile - execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/util/build_scripts/write_probin.py" --pa "${EOSparamfile} ${networkpcparamfile} - ${networkparamfile} ${VODEparamfile} ${integrationparamfile} " --use_namespace WORKING_DIRECTORY ${output_dir}/) - endif() - - set(photoionization_sources ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces/eos_data.cpp - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/interfaces/network_initialization.cpp - ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization/actual_eos_data.cpp - ${CMAKE_SOURCE_DIR}/src/networks/H_ionization_TInd/actual_network_data.cpp - ${output_dir}/extern_parameters.cpp PARENT_SCOPE) - - #below for NAUX execute_process(COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${PYTHONPATH}:${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null" python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/get_naux.py" --net "${networkdir}" --microphysics_path "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/" WORKING_DIRECTORY ${output_dir}/) execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/write_network.py" --header_template "${networkheadertemplatefile}" --header_output "${networkpropfile}" -s "${networkfile}" WORKING_DIRECTORY ${output_dir}/) @@ -146,7 +95,6 @@ function(setup_target_for_microphysics_compilation network_name output_dir) elseif (${network_name} STREQUAL "photoionization") #need these to write extern_parameters.H set(paramfile "${CMAKE_CURRENT_LIST_DIR}/_parameters") - # TODO(#1)(james471): See how to do eosparamfile correctly set(EOSparamfile "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/photoionization/_parameters") set(networkpcparamfile "${CMAKE_SOURCE_DIR}/src/networks/photoionization/_parameters") @@ -188,7 +136,6 @@ function(setup_target_for_microphysics_compilation network_name output_dir) ${CMAKE_SOURCE_DIR}/src/networks/photoionization/actual_network_data.cpp ${output_dir}/extern_parameters.cpp PARENT_SCOPE) - #below for NAUX execute_process(COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${PYTHONPATH}:${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null" python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/get_naux.py" --net "${networkdir}" --microphysics_path "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/" WORKING_DIRECTORY ${output_dir}/) execute_process(COMMAND python3 "${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null/write_network.py" --header_template "${networkheadertemplatefile}" --header_output "${networkpropfile}" -s "${networkfile}" WORKING_DIRECTORY ${output_dir}/) diff --git a/EOS/photoionization/_parameters b/EOS/photoionization/_parameters index c2eaca28fd..1da7528347 100644 --- a/EOS/photoionization/_parameters +++ b/EOS/photoionization/_parameters @@ -5,12 +5,12 @@ eos_gamma_default real 1.4 # define the specie names, and their masses and gammas species_1_name string "Electron" species_1_gamma real 5./3. -species_1_mass real 9.10938188e-28 +species_1_mass real 9.10938291e-28 species_2_name string "Hydrogen" species_2_gamma real 5./3. -species_2_mass real 1.67353251819e-24 +species_2_mass real 1.673532715291e-24 species_3_name string "Hydrogen_Ion" species_3_gamma real 5./3. -species_3_mass real 1.67353251819e-24 +species_3_mass real 1.672621777e-24 diff --git a/EOS/photoionization/actual_eos.H b/EOS/photoionization/actual_eos.H index 1254b83859..7680168fd3 100644 --- a/EOS/photoionization/actual_eos.H +++ b/EOS/photoionization/actual_eos.H @@ -7,13 +7,11 @@ #include "eos_type.H" #include -#include #include #include #include #include #include -#include const std::string eos_name = "multigamma"; @@ -68,27 +66,20 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void actual_eos (I input, T& state) { static_assert(std::is_same_v, "input must be either eos_input_rt or eos_input_re"); - const amrex::Real gasconstant = C::n_A * C::k_B; - const amrex::Real protonmass = C::m_p; - // Special gamma factors + constexpr amrex::Real protonmass = C::m_p; + amrex::Real sum_Abarinv = 0.0_rt; - amrex::Real sum_gammasinv = 0.0_rt; + amrex::Real sum_ni_fi_over_2 = 0.0_rt; amrex::Real rhotot = 0.0_rt; for (int n = 0; n < NumSpec; ++n) { - rhotot += state.xn[n]*spmasses[n]; - // std::cout << "actual_eos()--> xn[" << n << "]: " << state.xn[n] << ", spmasses[" << n << "]: " << spmasses[n] << std::endl; - } - // std::cout << "Density in EOS: " << rhotot << std::endl; - - for (int n = 0; n < NumSpec; ++n) { + rhotot += state.xn[n] * spmasses[n]; + sum_ni_fi_over_2 += state.xn[n] / (gammas[n] - 1.0); sum_Abarinv += state.xn[n]; - sum_gammasinv += (state.xn[n]*protonmass/rhotot) * (1.0/(gammas[n]-1.0)); } - // std::cout << "actual_eos()--> rhotot: " << rhotot << ", state.xn[0]: " << state.xn[0] << ", state.xn[1]: " << state.xn[1] << ", state.xn[2]: " << state.xn[2] << std::endl; + sum_ni_fi_over_2 /= sum_Abarinv; sum_Abarinv *= protonmass/rhotot; state.mu = 1.0 / sum_Abarinv; - sum_gammasinv /= sum_Abarinv; //------------------------------------------------------------------------- // For all EOS input modes EXCEPT eos_input_rt, first compute dens @@ -98,7 +89,6 @@ void actual_eos (I input, T& state) amrex::Real temp = NAN; amrex::Real dens = NAN; amrex::Real eint = NAN; - // std::cout << "Eos was called" << std::endl; switch (input) { case eos_input_rt: @@ -107,8 +97,7 @@ void actual_eos (I input, T& state) // We don't need to do anything here temp = state.T; dens = state.rho; - eint = sum_gammasinv * sum_Abarinv * gasconstant * state.T; - // std::cout << "actual_eos() --> eos_input_rt: temp: " << temp << ", dens: " << dens << ", eint: " << eint << ", sum_gammasinv: " << sum_gammasinv << ", sum_Abarinv: " << sum_Abarinv << std::endl; + eint = sum_ni_fi_over_2 * C::k_B * temp / rhotot; break; case eos_input_re: @@ -121,9 +110,8 @@ void actual_eos (I input, T& state) // stop the integration if the internal energy < 0 AMREX_ASSERT(state.e > 0.); - temp = state.e /( sum_gammasinv * gasconstant * sum_Abarinv); + temp = state.e * rhotot / (sum_ni_fi_over_2 * C::k_B); eint = state.e; - // std::cout << "actual_eos() --> eos_input_re: temp: " << temp << ", dens: " << dens << ", eint: " << eint << ", sum_gammasinv: " << sum_gammasinv << ", sum_Abarinv: " << sum_Abarinv << std::endl; } break; @@ -151,8 +139,8 @@ void actual_eos (I input, T& state) // stop the integration if the pressure < 0 AMREX_ASSERT(state.p > 0.); - eint = state.p * sum_gammasinv / dens; - temp = eint / (sum_gammasinv * gasconstant * sum_Abarinv); + eint = state.p * sum_ni_fi_over_2 / dens; + temp = eint * rhotot / (sum_ni_fi_over_2 * C::k_B); } break; @@ -203,18 +191,18 @@ void actual_eos (I input, T& state) } if constexpr (has_pressure::value) { - amrex::Real pressure = state.rho * eint / sum_gammasinv; + amrex::Real pressure = state.rho * eint / sum_ni_fi_over_2; state.p = pressure; state.dpdT = pressure / temp; state.dpdr = pressure / dens; - state.cs = std::sqrt((1.0 + 1.0/sum_gammasinv) * state.p /state.rho); + state.cs = std::sqrt((1.0 + 1.0/sum_ni_fi_over_2) * state.p /state.rho); if constexpr (has_G::value) { - state.G = 0.5 * (1.0 + (1.0 + 1.0/sum_gammasinv)); + state.G = 0.5 * (1.0 + (1.0 + 1.0/sum_ni_fi_over_2)); } } - amrex::Real dedT = sum_gammasinv * sum_Abarinv * gasconstant; + amrex::Real dedT = sum_ni_fi_over_2 * C::k_B / rhotot; amrex::Real dedr = 0.0_rt; if constexpr (has_energy::value) { state.dedT = dedT; diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H index 8e893b8f47..cf18f1f708 100644 --- a/integration/VODE/vode_dvhin.H +++ b/integration/VODE/vode_dvhin.H @@ -45,7 +45,6 @@ void dvhin (BurnT& state, DvodeT& vstate, amrex::Real& H0, int& NITER, int& IER) amrex::Real HUB = PT1 * TDIST; for (int i = 1; i <= int_neqs; ++i) { - // TODO(#3)(james471): This can probably be optimized to reduce branching amrex::Real atol{}; #if defined (PHOTOCHEMISTRY) if (i <= NumSpec) { diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H index 7f00ae0847..db664f5dc2 100644 --- a/integration/VODE/vode_dvstep.H +++ b/integration/VODE/vode_dvstep.H @@ -350,7 +350,6 @@ int dvstep (BurnT& state, DvodeT& vstate) break; } #endif -#endif #ifdef SDC // Abort the program. This is unimplemented @@ -358,6 +357,7 @@ int dvstep (BurnT& state, DvodeT& vstate) #endif } +#endif // The corrector has converged (NFLAG = 0). The local error test is // made and control passes to statement 500 if it fails. diff --git a/integration/utils/numerical_jacobian.H b/integration/utils/numerical_jacobian.H index 598566d26c..1d52531861 100644 --- a/integration/utils/numerical_jacobian.H +++ b/integration/utils/numerical_jacobian.H @@ -29,7 +29,7 @@ /// difference in terms of X and T and then convert the Jacobian /// elements to be in terms of X and e -// TODO(#2)(james471): This needs fixes for radiation +// TODO(james471): This needs fixes for radiation struct jac_info_t { amrex::Real h; diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H index 0cf5079ac0..1359b63fd5 100644 --- a/interfaces/burn_type.H +++ b/interfaces/burn_type.H @@ -170,7 +170,6 @@ struct burn_t }; -// TODO(james471): Add radiation variables to the printout inline std::ostream& operator<< (std::ostream& o, burn_t const& burn_state) { @@ -181,6 +180,12 @@ std::ostream& operator<< (std::ostream& o, burn_t const& burn_state) for (double X : burn_state.xn) { o << X << " "; } +#ifdef PHOTOCHEMISTRY + o << "rn = "; + for (double Rn : burn_state.rn) { + o << Rn << " "; + } +#endif o << std::endl; o << "(i, j, k) = " << burn_state.i << " " << burn_state.j << " " << burn_state.k << std::endl; #if NAUX_NET > 0 diff --git a/networks/rhs.H b/networks/rhs.H index c2780a49fb..da145f3d1b 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -20,7 +20,7 @@ #include #include -// TODO(#2)(james471): Check if this requires changes for radiation variables +// TODO(james471): Check if this requires changes for radiation variables #ifdef NEW_NETWORK_IMPLEMENTATION From e177e7a729c8e4ee3521bc7da295a4e9b6c4b550 Mon Sep 17 00:00:00 2001 From: James471 Date: Tue, 12 May 2026 19:06:26 +1000 Subject: [PATCH 8/8] Refactor radiation variable names --- integration/VODE/vode_dvode.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/integration/VODE/vode_dvode.H b/integration/VODE/vode_dvode.H index 3dc8a37e84..81db1d75fc 100644 --- a/integration/VODE/vode_dvode.H +++ b/integration/VODE/vode_dvode.H @@ -82,7 +82,7 @@ int dvode (BurnT& state, DvodeT& vstate) vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc; vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1); #ifdef PHOTOCHEMISTRY - for (int i = 1; i <= NumChemActiveRadGroups; i+=MicrophysicsNumRadVarsPerGroup) { + for (int i = 1; i <= NumChemBands; i+=MicrophysicsNumRadVarsPerGroup) { int rad_n_index = NumSpec + 2 + (i-1) * MicrophysicsNumRadVarsPerGroup; vstate.ewt(rad_n_index) = vstate.rtol_rad_num * std::abs(vstate.yh(rad_n_index,1)) + vstate.atol_rad_num; for (int j = 2; j <= MicrophysicsNumRadVarsPerGroup; ++j){ @@ -165,7 +165,7 @@ int dvode (BurnT& state, DvodeT& vstate) vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc; vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1); #ifdef PHOTOCHEMISTRY - for (int i = 1; i <= NumChemActiveRadGroups; i+=MicrophysicsNumRadVarsPerGroup) { + for (int i = 1; i <= NumChemBands; i+=MicrophysicsNumRadVarsPerGroup) { int rad_n_index = NumSpec + 2 + (i-1) * MicrophysicsNumRadVarsPerGroup; vstate.ewt(rad_n_index) = vstate.rtol_rad_num * std::abs(vstate.yh(rad_n_index,1)) + vstate.atol_rad_num; for (int j = 2; j <= MicrophysicsNumRadVarsPerGroup; ++j){