Update singularity-eos and singularity-opac

CMakeLists.txt

@@ -123,9 +123,13 @@ if(PHOEBUS_ENABLE_HDF5)
       "-DPHOEBUS_ENABLE_HDF5=OFF".)
   endif()
   set(SINGULARITY_USE_HDF5 ON CACHE BOOL "" FORCE)
+  set(SINGULARITY_USE_SPINER ON CACHE BOOL "" FORCE)
+  set(SINGULARITY_USE_SPINER_WITH_HDF5 ON CACHE BOOL "" FORCE)
   set(PARTHENON_DISABLE_HDF5 OFF CACHE BOOL "" FORCE)
 else()
   set(SINGULARITY_USE_HDF5 OFF CACHE BOOL "" FORCE)
+  set(SINGULARITY_USE_SPINER OFF CACHE BOOL "" FORCE)
+  set(SINGULARITY_USE_SPINER_WITH_HDF5 OFF CACHE BOOL "" FORCE)
   set(PARTHENON_DISABLE_HDF5 ON CACHE BOOL "" FORCE)
 endif()
 
@@ -159,23 +163,16 @@ add_subdirectory(external/parthenon parthenon)
 message("Configuring singularity")
 set(SINGULARITY_USE_KOKKOS ON CACHE BOOL "" FORCE)
 set(SINGULARITY_USE_FORTRAN OFF CACHE BOOL "" FORCE)
-set(SINGULARITY_SUBMODULE_MODE ON CACHE BOOL "" FORCE)
+set(SINGULARITY_FORCE_SUBMODULE_MODE ON CACHE BOOL "" FORCE)
 set(SINGULARITY_USE_KOKKOSKERNELS OFF CACHE BOOL "" FORCE)
 set(SINGULARITY_BUILD_CLOSURE OFF CACHE BOOL "" FORCE)
 
 # Kill cmake's package registry because it can interfere
 set(CMAKE_FIND_PACKAGE_NO_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
 set(CMAKE_FIND_PACKAGE_NO_SYSTEM_PACKAGE_REGISTRY ON CACHE BOOL "" FORCE)
 
-# singularity eos
-message(STATUS "Patching mpark::variant to support GPUs")
-execute_process(COMMAND patch -N -s -V never
-  ${CMAKE_CURRENT_SOURCE_DIR}/../external/singularity-eos/utils/variant/include/mpark/variant.hpp
-  ${CMAKE_CURRENT_SOURCE_DIR}/../external/singularity-eos/utils/cuda_compatibility.patch
-)
-
 # message("Configuring singularity-eos")
-# add_subdirectory(external/singularity-eos singularity-eos)
+add_subdirectory(external/singularity-eos singularity-eos)
 
 # singularity opac
 # message("Configuring singularity-opac")

src/CMakeLists.txt

@@ -154,6 +154,7 @@ configure_file(compile_constants.hpp.in generated/compile_constants.hpp @ONLY)
 target_include_directories(phoebus PRIVATE
   $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
   $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/generated>
+  # JMM: Not needed, but can't hurt
   "${CMAKE_CURRENT_SOURCE_DIR}/../external/singularity-eos"
   "${CMAKE_CURRENT_SOURCE_DIR}/../external/singularity-opac"
   "${CMAKE_CURRENT_SOURCE_DIR}/../external/singularity-eos/utils"
@@ -183,7 +184,7 @@ endif()
 # target_compile_options(phoebus PUBLIC -G -lineinfo)
 
 target_link_libraries(phoebus PRIVATE parthenon
-  # singularity-eos::singularity-eos
+  singularity-eos::singularity-eos
   # singularity-opac::singularity-opac
 )
 if (PHOEBUS_ENABLE_HDF5)

src/microphysics/eos_phoebus/eos_phoebus.cpp

@@ -20,9 +20,6 @@
 // parthenon includes
 #include <parthenon/package.hpp>
 
-// singularity includes
-#include <singularity-eos/eos/eos.hpp>
-
 // phoebus includes
 #include "microphysics/eos_phoebus/eos_phoebus.hpp"
 #include "phoebus_utils/unit_conversions.hpp"

src/microphysics/eos_phoebus/eos_phoebus.hpp

@@ -17,7 +17,13 @@
 #include <memory>
 #include <parthenon/package.hpp>
 
-#include <singularity-eos/eos/eos.hpp>
+#include <singularity-eos/eos/eos_ideal.hpp>
+#include <singularity-eos/eos/eos_variant.hpp>
+#include <singularity-eos/eos/modifiers/eos_unitsystem.hpp>
+
+#ifdef SPINER_USE_HDF
+#include <singularity-eos/eos/eos_stellar_collapse.hpp>
+#endif
 
 using namespace parthenon::package::prelude;
 

src/monopole_gr/monopole_gr_base.hpp

@@ -72,7 +72,7 @@ enum Gradients {
 };
 
 // TODO(JMM): Do we want this?
-using Radius = Spiner::RegularGrid1D;
+using Radius = Spiner::RegularGrid1D<Real>;
 
 template <typename Array_t>
 PORTABLE_INLINE_FUNCTION Real Interpolate(const Real r, const Array_t &A,

src/pgen/progenitor.cpp

@@ -22,6 +22,8 @@
 #include "pgen/pgen.hpp"
 #include "progenitor/progenitordata.hpp"
 
+using DataBox = Spiner::DataBox<Real>;
+
 // Homologously collapsing star.
 namespace homologous {
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
@@ -67,15 +69,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   // Interpolated Initial Data
   auto progenitor_pkg = pmb->packages.Get("progenitor");
 
-  auto adm_density = progenitor_pkg->Param<Spiner::DataBox>("adm_density");
-  auto adm_momentum = progenitor_pkg->Param<Spiner::DataBox>("adm_momentum");
-  auto S_adm = progenitor_pkg->Param<Spiner::DataBox>("S_adm");
-  auto Srr_adm = progenitor_pkg->Param<Spiner::DataBox>("Srr_adm");
+  auto adm_density = progenitor_pkg->Param<DataBox>("adm_density");
+  auto adm_momentum = progenitor_pkg->Param<DataBox>("adm_momentum");
+  auto S_adm = progenitor_pkg->Param<DataBox>("S_adm");
+  auto Srr_adm = progenitor_pkg->Param<DataBox>("Srr_adm");
 
-  auto mass_density_dev = progenitor_pkg->Param<Spiner::DataBox>("mass_density_dev");
-  auto Ye_dev = progenitor_pkg->Param<Spiner::DataBox>("Ye_dev");
-  auto velocity_dev = progenitor_pkg->Param<Spiner::DataBox>("velocity_dev");
-  auto pressure_dev = progenitor_pkg->Param<Spiner::DataBox>("pressure_dev");
+  auto mass_density_dev = progenitor_pkg->Param<DataBox>("mass_density_dev");
+  auto Ye_dev = progenitor_pkg->Param<DataBox>("Ye_dev");
+  auto velocity_dev = progenitor_pkg->Param<DataBox>("velocity_dev");
+  auto pressure_dev = progenitor_pkg->Param<DataBox>("pressure_dev");
 
   // MonopoleGR
   static bool monopole_initialized = false;

src/pgen/torus.cpp

@@ -37,6 +37,7 @@
 typedef Kokkos::Random_XorShift64_Pool<> RNGPool;
 
 namespace torus {
+using DataBox = Spiner::DataBox<Real>;
 
 using pc = parthenon::constants::PhysicalConstants<parthenon::constants::CGS>;
 
@@ -78,8 +79,8 @@ class GasRadTemperatureResidual {
 class EnthalpyResidual {
  public:
   KOKKOS_FUNCTION
-  EnthalpyResidual(const Real hm1, const Spiner::DataBox T, const Real Ye,
-                   const Real h_min_sc, const EOS &eos)
+  EnthalpyResidual(const Real hm1, const DataBox T, const Real Ye, const Real h_min_sc,
+                   const EOS &eos)
       : hm1_(hm1), T_(T), Ye_(Ye), h_min_sc_(h_min_sc), eos_(eos) {}
 
   KOKKOS_INLINE_FUNCTION
@@ -90,7 +91,7 @@ class EnthalpyResidual {
 
  private:
   const Real hm1_;
-  const Spiner::DataBox T_;
+  const DataBox T_;
   const Real Ye_;
   const Real h_min_sc_;
   const EOS &eos_;
@@ -121,10 +122,9 @@ void ComputeBetas(Mesh *pmesh, const Real rho_min_bnorm, Real &beta_min_global,
                   Real &beta_pmax);
 KOKKOS_FUNCTION
 void GetStateFromEnthalpy(const EOS &eos, const EosType eos_type, const Real hm1,
-                          const Spiner::DataBox rho, const Spiner::DataBox temp,
-                          const Real Ye, const Real h_min_sc, const Real kappa,
-                          const Real gam, const Real Cv, const Real rho_rmax,
-                          Real &rho_out, Real &u_out);
+                          const DataBox rho, const DataBox temp, const Real Ye,
+                          const Real h_min_sc, const Real kappa, const Real gam,
+                          const Real Cv, const Real rho_rmax, Real &rho_out, Real &u_out);
 // ----------------------------------------------------------------------
 
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
@@ -239,8 +239,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
 
   RNGPool rng_pool(seed);
 
-  Spiner::DataBox rho_h(nsamps);
-  Spiner::DataBox temp_h(nsamps);
+  DataBox rho_h(nsamps);
+  DataBox temp_h(nsamps);
 
   const Real rho_min = pmb->packages.Get("eos")->Param<Real>("rho_min");
   const Real rho_max = pmb->packages.Get("eos")->Param<Real>("rho_max");
@@ -622,10 +622,10 @@ Real log_enthalpy(const Real r, const Real th, const Real a, const Real rin, con
  **/
 KOKKOS_FUNCTION
 void GetStateFromEnthalpy(const EOS &eos, const EosType eos_type, const Real hm1,
-                          const Spiner::DataBox rho, const Spiner::DataBox temp,
-                          const Real Ye, const Real h_min_sc, const Real kappa,
-                          const Real gam, const Real Cv, const Real rho_rmax,
-                          Real &rho_out, Real &u_out) {
+                          const DataBox rho, const DataBox temp, const Real Ye,
+                          const Real h_min_sc, const Real kappa, const Real gam,
+                          const Real Cv, const Real rho_rmax, Real &rho_out,
+                          Real &u_out) {
   if (eos_type == IdealGas) { // Ideal Gas
     rho_out = std::pow(hm1 * (gam - 1.) / (kappa * gam), 1. / (gam - 1.));
     u_out = kappa * std::pow(rho_out, gam) / (gam - 1.) / rho_rmax;

src/phoebus_utils/adiabats.hpp

@@ -28,6 +28,7 @@
 namespace Adiabats {
 
 using Microphysics::EOS::EOS;
+using DataBox = Spiner::DataBox<Real>;
 
 template <typename EOS>
 inline void GetRhoBounds(const EOS &eos, const Real rho_min, const Real rho_max,
@@ -71,7 +72,7 @@ inline void GetRhoBounds(const EOS &eos, const Real rho_min, const Real rho_max,
 }
 
 // sample log rho
-inline void SampleRho(Spiner::DataBox lrho, const Real lrho_min, const Real lrho_max,
+inline void SampleRho(DataBox lrho, const Real lrho_min, const Real lrho_max,
                       const int n_samps) {
   const Real dlrho = (lrho_max - lrho_min) / n_samps;
 
@@ -84,9 +85,9 @@ inline void SampleRho(Spiner::DataBox lrho, const Real lrho_min, const Real lrho
  * Given Ye and a target entropy, compute density and temperature of constant entropy
  **/
 template <typename EOS>
-inline void ComputeAdiabats(Spiner::DataBox lrho, Spiner::DataBox temp, const EOS &eos,
-                            const Real Ye, const Real S0, const Real T_min,
-                            const Real T_max, const int n_samps) {
+inline void ComputeAdiabats(DataBox lrho, DataBox temp, const EOS &eos, const Real Ye,
+                            const Real S0, const Real T_min, const Real T_max,
+                            const int n_samps) {
 
   const Real guess0 = (T_max - T_min) / 2.0;
   const Real epsilon = std::numeric_limits<Real>::epsilon();
@@ -110,8 +111,8 @@ inline void ComputeAdiabats(Spiner::DataBox lrho, Spiner::DataBox temp, const EO
  * Find the minimum enthalpy along an adiabat as computed above
  **/
 template <typename EOS>
-inline Real MinEnthalpy(Spiner::DataBox lrho, Spiner::DataBox temp, const Real Ye,
-                        const EOS &eos, const int n_samps) {
+inline Real MinEnthalpy(DataBox lrho, DataBox temp, const Real Ye, const EOS &eos,
+                        const int n_samps) {
   Real min_enthalpy = 1e30;
   for (int i = 0; i < n_samps; i++) {
     const Real Rho = std::pow(10.0, lrho(i));

src/phoebus_utils/phoebus_interpolation.hpp

@@ -30,7 +30,7 @@ namespace interpolation {
 
 using namespace parthenon::package::prelude;
 using parthenon::Coordinates_t;
-using Spiner::weights_t;
+using weights_t = Spiner::weights_t<Real>;
 
 /// Base class for providing interpolation methods on uniformly spaced data.
 /// Constructor is provided with spacing, number of support points, and desired

src/progenitor/progenitordata.cpp

@@ -25,18 +25,18 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   const int npoints = mydata[0].size();
 
   // DataBoxes for interpolation functions
-  Spiner::DataBox r(npoints);
-  Spiner::DataBox mass_density(npoints);
-  Spiner::DataBox temp(npoints);
-  Spiner::DataBox Ye(npoints);
-  Spiner::DataBox specific_internal_energy(npoints);
-  Spiner::DataBox velocity(npoints);
-  Spiner::DataBox pressure(npoints);
-
-  Spiner::DataBox adm_density(npoints);
-  Spiner::DataBox adm_momentum(npoints);
-  Spiner::DataBox S_adm(npoints);
-  Spiner::DataBox Srr_adm(npoints);
+  Spiner::DataBox<Real> r(npoints);
+  Spiner::DataBox<Real> mass_density(npoints);
+  Spiner::DataBox<Real> temp(npoints);
+  Spiner::DataBox<Real> Ye(npoints);
+  Spiner::DataBox<Real> specific_internal_energy(npoints);
+  Spiner::DataBox<Real> velocity(npoints);
+  Spiner::DataBox<Real> pressure(npoints);
+
+  Spiner::DataBox<Real> adm_density(npoints);
+  Spiner::DataBox<Real> adm_momentum(npoints);
+  Spiner::DataBox<Real> S_adm(npoints);
+  Spiner::DataBox<Real> Srr_adm(npoints);
 
   // Fill in the variables
   // Requires same grid for primitive and adm quantities in the input table

tst/unit/geometry/CMakeLists.txt

@@ -47,6 +47,7 @@ configure_file(../../../src/compile_constants.hpp.in generated/compile_constants
 target_link_libraries(geometry_unit_tests
   PRIVATE
   parthenon
+  singularity-eos::singularity-eos
   Catch2::Catch2
   Kokkos::kokkos)
 

tst/unit/monopole_gr/CMakeLists.txt

@@ -33,6 +33,7 @@ target_link_libraries(monopole_gr_unit_tests
   parthenon
   Catch2::Catch2
   Kokkos::kokkos
+  singularity-eos::singularity-eos
   )
 
 target_include_directories(monopole_gr_unit_tests PRIVATE

tst/unit/phoebus_utils/CMakeLists.txt

@@ -48,6 +48,7 @@ target_link_libraries(phoebus_utils_unit_tests
     parthenon
     Catch2::Catch2
     Kokkos::kokkos
+    singularity-eos::singularity-eos
     )
 if (PHOEBUS_ENABLE_HDF5)
   target_link_libraries(phoebus_utils_unit_tests PRIVATE hdf5::hdf5 hdf5::hdf5_hl )

tst/unit/phoebus_utils/test_adiabats.cpp

@@ -48,12 +48,12 @@ TEST_CASE("ADIABATS", "[compute_adiabats]") {
     const int nsamps = 180;
     Spiner::DataBox rho_d(Spiner::AllocationTarget::Device, nsamps);
     Spiner::DataBox temp_d(Spiner::AllocationTarget::Device, nsamps);
-    const Real lrho_min = eos.lRhoMin();
-    const Real lrho_max = eos.lRhoMax();
-    const Real rho_min = std::pow(10.0, lrho_min);
-    const Real rho_max = std::pow(10.0, lrho_max);
-    const Real T_min = std::pow(10.0, eos.lTMin());
-    const Real T_max = std::pow(10.0, eos.lTMax());
+    const Real rho_min = eos.rhoMin();
+    const Real rho_max = eos.rhoMax();
+    const Real lrho_min = std::log10(rho_min);
+    const Real lrho_max = std::log10(rho_max);
+    const Real T_min = eos.TMin();
+    const Real T_max = eos.TMax();
 
     const Real Ye = 0.1;
     const Real S0 = 20.0;

tst/unit/radiation/CMakeLists.txt

@@ -21,7 +21,9 @@ target_link_libraries(radiation_unit_tests
   PRIVATE
   parthenon
   Catch2::Catch2
-  Kokkos::kokkos)
+  Kokkos::kokkos
+  singularity-eos::singularity-eos
+  )
 
 target_include_directories(radiation_unit_tests PRIVATE
   "${CMAKE_CURRENT_SOURCE_DIR}/../../../src"