diff --git a/CMakeLists.txt b/CMakeLists.txt index f4eb5cc45b..a1be9c7455 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) #---------------------------------------------------------------------------------------------------------------------- @@ -88,6 +88,53 @@ 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 "photoionization") + #need these to write extern_parameters.H + set(paramfile "${CMAKE_CURRENT_LIST_DIR}/_parameters") + 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}/) 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..1da7528347 --- /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.10938291e-28 + +species_2_name string "Hydrogen" +species_2_gamma real 5./3. +species_2_mass real 1.673532715291e-24 + +species_3_name string "Hydrogen_Ion" +species_3_gamma real 5./3. +species_3_mass real 1.672621777e-24 diff --git a/EOS/photoionization/actual_eos.H b/EOS/photoionization/actual_eos.H new file mode 100644 index 0000000000..7680168fd3 --- /dev/null +++ b/EOS/photoionization/actual_eos.H @@ -0,0 +1,224 @@ +#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 + +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"); + constexpr amrex::Real protonmass = C::m_p; + + amrex::Real sum_Abarinv = 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]; + sum_ni_fi_over_2 += state.xn[n] / (gammas[n] - 1.0); + sum_Abarinv += state.xn[n]; + } + sum_ni_fi_over_2 /= sum_Abarinv; + sum_Abarinv *= protonmass/rhotot; + state.mu = 1.0 / 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; + 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_ni_fi_over_2 * C::k_B * temp / rhotot; + 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 * rhotot / (sum_ni_fi_over_2 * C::k_B); + eint = state.e; + } + + 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_ni_fi_over_2 / dens; + temp = eint * rhotot / (sum_ni_fi_over_2 * C::k_B); + } + 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_ni_fi_over_2; + state.p = pressure; + + state.dpdT = pressure / temp; + state.dpdr = pressure / dens; + 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_ni_fi_over_2)); + } + } + + 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; + 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]; diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H index 96c18b5756..cf18f1f708 100644 --- a/integration/VODE/vode_dvhin.H +++ b/integration/VODE/vode_dvhin.H @@ -45,7 +45,24 @@ 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; + 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..81db1d75fc 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 <= 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){ + 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 <= 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){ + 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..db664f5dc2 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 @@ -308,6 +309,56 @@ int dvstep (BurnT& state, DvodeT& vstate) } +#ifdef PHOTOCHEMISTRY + 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 + // 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; + } +#endif + +#ifdef SDC + // Abort the program. This is unimplemented + amrex::Abort("SDC does not yet support radiation updates"); +#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/VODE/vode_type.H b/integration/VODE/vode_type.H index 3c246bfb50..8152573492 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,10 +135,17 @@ 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; + // Reduced speed of light + amrex::Real c_hat; + // Integration array amrex::Array1D y; diff --git a/integration/_parameters b/integration/_parameters index a11b5f3ea8..3b79841575 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -23,16 +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 number density +rtol_rad_num real 1.e-10 + +# relative tolerance for radiation flux +rtol_rad_flux real 1.e-10 + +# absolute tolerance for radiation number density +atol_rad_num real 1.e-10 + +# absolute tolerance for radiation flux +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. @@ -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..75d070af6c 100644 --- a/integration/integrator_data.H +++ b/integration/integrator_data.H @@ -5,12 +5,22 @@ // 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; +#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; @@ -38,6 +48,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 +58,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..1ef864c241 100644 --- a/integration/integrator_setup_strang.H +++ b/integration/integrator_setup_strang.H @@ -20,6 +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 }; @@ -38,6 +42,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 +59,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 +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; @@ -180,6 +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) { @@ -210,6 +244,14 @@ 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 << "(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 << " "; + } +#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 +259,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..de0c83700c 100644 --- a/integration/integrator_type.H +++ b/integration/integrator_type.H @@ -58,6 +58,12 @@ 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); + } + state.c_hat = int_state.c_hat; +#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..63396c7f87 100644 --- a/integration/integrator_type_strang.H +++ b/integration/integrator_type_strang.H @@ -18,6 +18,13 @@ 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); + } + state.c_hat = int_state.c_hat; +#endif + if (integrator_rp::scale_system) { state.e *= state.e_scale; } @@ -38,6 +45,13 @@ 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]; + } + int_state.c_hat = state.c_hat; +#endif + if (integrator_rp::scale_system) { int_state.y(net_ienuc) /= state.e_scale; } @@ -65,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..1d52531861 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(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 316053530d..1359b63fd5 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,10 @@ struct burn_t #if NAUX_NET > 0 amrex::Real aux[NumAux]{}; #endif +#ifdef PHOTOCHEMISTRY + amrex::Real rn[NumRadEqs]{}; + amrex::Real c_hat; +#endif // scaling used to make e we integrate dimensionless amrex::Real e_scale{}; @@ -158,7 +170,6 @@ struct burn_t }; - inline std::ostream& operator<< (std::ostream& o, burn_t const& burn_state) { @@ -169,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 @@ -233,6 +250,12 @@ 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]; + } + burn_state.c_hat = eos_state.c_hat; +#endif burn_state.cv = eos_state.cv; @@ -271,6 +294,13 @@ 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]; + } + eos_state.c_hat = burn_state.c_hat; +#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..4de946910e 100644 --- a/interfaces/eos.H +++ b/interfaces/eos.H @@ -1,6 +1,7 @@ #ifndef EOS_H #define EOS_H +#include #include #include #include @@ -72,7 +73,6 @@ 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)); - actual_eos(eos_input_rt, state); has_been_reset = true; @@ -216,6 +216,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; } @@ -424,7 +431,6 @@ void eos (const I input, T& state) eos_override(state); // Call the EOS. - if (!has_been_reset) { actual_eos(input, state); } 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..da145f3d1b 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -20,6 +20,8 @@ #include #include +// TODO(james471): Check if this requires changes for radiation variables + #ifdef NEW_NETWORK_IMPLEMENTATION // Forward declarations 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):