AMReX-Astro · James471 · Mar 2, 2026 · Mar 2, 2026 · Mar 2, 2026 · Mar 6, 2026
diff --git 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
@@ -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
@@ -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
@@ -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 <AMReX.H>
+#include <network.H>
+#include <fundamental_constants.H>
+#include <extern_parameters.H>
+#include <cmath>
+#include <actual_eos_data.H>
+
+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 <typename I>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+bool is_input_valid (I input)
+{
+  static_assert(std::is_same_v<I, eos_input_t>, "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 <typename I, typename T>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void actual_eos (I input, T& state)
+{
+    static_assert(std::is_same_v<I, eos_input_t>, "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<T>::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<T>::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<T>::value) {
+        state.e = eint;
+    }
+
+    if constexpr (has_pressure<T>::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<T>::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<T>::value) {
+        state.dedT = dedT;
+        state.dedr = dedr;
+        if constexpr (has_pressure<T>::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
@@ -0,0 +1,11 @@
+#ifndef actual_eos_data_H
+#define actual_eos_data_H
+
+#include <AMReX.H>
+#include <AMReX_REAL.H>
+#include <network.H>
+
+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
@@ -0,0 +1,4 @@
+#include <actual_eos_data.H>
+
+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
@@ -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
@@ -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 {