Add singularity EOS and unit conversion

cmake/InstallTPLs.cmake

@@ -222,6 +222,69 @@ set_property(GLOBAL PROPERTY SPHERAL_FP_TPLS ${SPHERAL_FP_TPLS})
 set_property(GLOBAL PROPERTY SPHERAL_FP_DIRS ${SPHERAL_FP_DIRS})
 
 message("-----------------------------------------------------------------------------")
+if (singularity_eos_DIR)
+  find_package(ports-of-call REQUIRED NO_DEFAULT_PATH PATHS ${ports_of_call_DIR})
+  if (ports-of-call_FOUND)
+    message("Found ports-of-call External Package")
+    list(APPEND SPHERAL_BLT_DEPENDS ports-of-call::ports-of-call)
+    list(APPEND SPHERAL_FP_TPLS ports-of-call::ports-of-call)
+    list(APPEND SPHERAL_FP_DIRS ${ports_of_call_DIR})
+    blt_convert_to_system_includes(TARGET ports-of-call::ports-of-call)
+  endif()
+
+  find_package(Kokkos REQUIRED NO_DEFAULT_PATH PATHS ${kokkos_DIR})
+  if (Kokkos_FOUND)
+    message("Found kokkos External Package")
+    list(APPEND SPHERAL_FP_DIRS ${ports_of_call_DIR})
+    foreach(_tar IN ITEMS Kokkos::kokkos Kokkos::kokkoscore)
+      list(APPEND SPHERAL_BLT_DEPENDS ${_tar})
+      list(APPEND SPHERAL_FP_TPLS ${_tar})
+      blt_convert_to_system_includes(TARGET ${_tar})
+    endforeach()
+  endif()
+
+  find_package(KokkosKernels REQUIRED NO_DEFAULT_PATH PATHS ${kokkos_kernels_DIR})
+  if (KokkosKernels_FOUND)
+    message("Found KokkosKernels External Package")
+    list(APPEND SPHERAL_BLT_DEPENDS Kokkos::kokkoskernels)
+    list(APPEND SPHERAL_FP_TPLS Kokkos::kokkoskernels)
+    list(APPEND SPHERAL_FP_DIRS ${kokkos_kernels_DIR})
+    blt_convert_to_system_includes(TARGET Kokkos::kokkoskernels)
+  endif()
+
+  find_package(spiner REQUIRED NO_DEFAULT_PATH PATHS ${spiner_DIR})
+  if (spiner_FOUND)
+    message("Found spiner External Package")
+    list(APPEND SPHERAL_BLT_DEPENDS spiner::spiner)
+    list(APPEND SPHERAL_FP_TPLS spiner::spiner)
+    list(APPEND SPHERAL_FP_DIRS ${spiner_DIR})
+    blt_convert_to_system_includes(TARGET spiner::spiner)
+  endif()
+
+  find_package(mpark_variant REQUIRED NO_DEFAULT_PATH PATHS ${mpark_variant_DIR})
+  if (mpark_variant_FOUND)
+    message("Found mpark_variant External Package")
+    list(APPEND SPHERAL_BLT_DEPENDS mpark_variant)
+    list(APPEND SPHERAL_FP_TPLS mpark_variant)
+    list(APPEND SPHERAL_FP_DIRS ${mpark_variant_DIR})
+    blt_convert_to_system_includes(TARGET mpark_variant)
+  endif()
+
+  find_package(singularity-eos REQUIRED COMPONENTS Interface NO_DEFAULT_PATH PATHS ${singularity_eos_DIR})
+  if(singularity-eos_FOUND)
+    message("Found Singularity-EOS External Package")
+    list(APPEND SPHERAL_FP_DIRS ${singularity_eos_DIR})
+    list(APPEND SPHERAL_BLT_DEPENDS singularity-eos::singularity-eos_Interface singularity-eos::singularity-eos_Common)
+    list(APPEND SPHERAL_FP_TPLS singularity-eos::singularity-eos_Interface singularity-eos::singularity-eos_Common)
+    blt_convert_to_system_includes(TARGET singularity-eos::singularity-eos_Interface)
+    blt_convert_to_system_includes(TARGET singularity-eos::singularity-eos_Common)
+  endif()
+
+  list(APPEND SPHERAL_EXTERN_LIBS eospac)
+endif()
+
+message("-----------------------------------------------------------------------------")
+
 # In case we start using find_package on Silo, we should save the silo_DIR path
 set(CONFIG_SILO_DIR "${silo_DIR}" CACHE PATH "Configuration Silo directory")
 # TPLs that must be imported

cmake/tpl/eospac.cmake

@@ -0,0 +1 @@
+set(eospac_libs libeospac6.a)

scripts/devtools/tpl-manager.py

@@ -263,13 +263,34 @@ def set_concretize(loader):
                               "blas": ["openblas"],
                               "lapack": ["openblas"]})
         if (provider_dict):
-            self.config_env_providers(provider_dict)
+            # Extract compiler from spec (e.g., "gcc" from "spheral+mpi%gcc")
+            compiler_name = None
+            spec_str = self.args.spec
+            if "%" in spec_str:
+                compiler_name = spec_str.split("%")[-1].split()[0]
+            self.config_env_providers(provider_dict, compiler=compiler_name)
         self.args.add_spec = True
 
-    def config_env_providers(self, config_dict):
+    def config_env_providers(self, config_dict, compiler=None):
         env_file = os.path.join(self.env_dir, "spack.yaml")
         def set_providers(loader):
             new_dict = {"all": {"providers": config_dict}}
+            if compiler:
+                # WORKAROUND: Spack 1.1.0 has a bug where compiler_mixing:false
+                # does not prevent llvm from being assigned to some packages when
+                # both gcc and llvm are detected as externals (spack/spack#51995).
+                # This was fixed upstream in spack/spack#52015 (commit bf1345e6,
+                # merged 2026-03-03) and should ship in spack >=1.2.0.
+                # Remove this require block once we update past spack 1.1.0.
+                #
+                # Use conditional require syntax so packages without compiled
+                # code (e.g. externals) are not affected. See spack docs:
+                # "require: %clang will fail for packages without compiled code"
+                new_dict["all"]["require"] = [
+                    f"%[when=%c]c={compiler} "
+                    f"%[when=%cxx]cxx={compiler} "
+                    f"%[when=%fortran]fortran={compiler}"
+                ]
             loader["spack"]["packages"].update(new_dict)
             return loader
         self.modify_env_file(env_file, set_providers)

scripts/spack/spack_repo/spheral/packages/spheral/package.py

@@ -45,6 +45,7 @@ class Spheral(CachedCMakePackage, CudaPackage, ROCmPackage):
     variant('opensubdiv', default=True, description='Enable use of opensubdiv to do refinement.')
     variant('network', default=True, description='Disable to build Spheral from a local buildcache.')
     variant('sundials', default=True, when="+mpi", description='Enable use of SUNDIALS solvers.')
+    variant('singularity-eos', default=True, description='Build Singularity EOS.')
     variant('leos', default=LEOSpresent, when="+mpi", description='Build LEOS package.')
 
     # -------------------------------------------------------------------------
@@ -99,6 +100,15 @@ class Spheral(CachedCMakePackage, CudaPackage, ROCmPackage):
     depends_on('sundials@7.0.0 ~shared cxxstd=17 cppflags="-fPIC"', type='build', when='+sundials')
     depends_on('sundials build_type=Debug', when='+sundials build_type=Debug')
 
+    # Doing all these explicitly to avoid linker issues
+    depends_on('singularity-eos@1.10.0 +spiner +hdf5 ~fortran cppflags="-fPIC" cflags="-fPIC"', type='build', when="+singularity-eos")
+    depends_on('spiner cppflags="-fPIC" cflags="-fPIC"', when='+singularity-eos')
+    depends_on('kokkos +pic cppflags="-fPIC" cflags="-fPIC"', when='+singularity-eos')
+    depends_on('kokkos-kernels cppflags="-fPIC" cflags="-fPIC"', when='+singularity-eos')
+    depends_on('eospac cppflags="-fPIC" cflags="-fPIC"', when='+singularity-eos')
+    depends_on('ports-of-call cppflags="-fPIC" cflags="-fPIC"', when='+singularity-eos')
+    depends_on('mpark-variant', when='+singularity-eos')
+
     # Forward MPI Variants
     mpi_tpl_list = ["hdf5", "conduit", "axom", "adiak~shared", "chai", "umpire"]
     for ctpl in mpi_tpl_list:
@@ -297,6 +307,19 @@ def initconfig_package_entries(self):
         if spec.satisfies("+sundials"):
             entries.append(cmake_cache_path('sundials_DIR', spec['sundials'].prefix))
 
+        if spec.satisfies("+singularity-eos"):
+            entries.append(cmake_cache_path('ports_of_call_DIR', spec['ports-of-call'].prefix))
+            entries.append(cmake_cache_path('kokkos_DIR', spec['kokkos'].prefix))
+            entries.append(cmake_cache_path('kokkos_kernels_DIR', spec['kokkos-kernels'].prefix))
+            entries.append(cmake_cache_path('spiner_DIR', spec['spiner'].prefix))
+            entries.append(cmake_cache_path('singularity_eos_DIR', spec['singularity-eos'].prefix))
+            entries.append(cmake_cache_path('eospac_DIR', spec['eospac'].prefix))
+            entries.append(cmake_cache_path('mpark_variant_DIR', spec['mpark-variant'].prefix + '/lib/cmake/mpark_variant'))
+
+            # These are needed for singularity to work correctly
+            entries.append(cmake_cache_path('EOSPAC_INCLUDE_DIR', spec['eospac'].prefix.include))
+            entries.append(cmake_cache_path('EOSPAC_LIBRARY', spec['eospac'].prefix.lib.join('libeospac6.a')))
+
         if spec.satisfies("+leos"):
             entries.append(cmake_cache_path('leos_DIR', spec['leos'].prefix))
             entries.append(cmake_cache_option('SPHERAL_ENABLE_LEOS', True))

src/Material/CMakeLists.txt

@@ -9,6 +9,7 @@ set(Material_inst
 
 set(Material_sources
     PhysicalConstants.cc
+    UnitConversion.cc
     )
 
 set(Material_headers
@@ -23,6 +24,7 @@ set(Material_headers
     PhysicalConstantsInline.hh
     PolytropicEquationOfState.hh
     PolytropicEquationOfStateInline.hh
+    UnitConversion.hh
     )
 
 spheral_install_python_files(

src/Material/MaterialEquationsOfState.py

@@ -73,7 +73,6 @@ def __init__(self,
                                                    maximumPressure)
         return
 
-
 #-------------------------------------------------------------------------------
 # Stiffened Gas EOS
 #-------------------------------------------------------------------------------
@@ -91,16 +90,34 @@ def __init__(self,
                                       P0,
                                       Cv,
                                       self._units,
+                                      minimumPressure,
+                                      maximumPressure)
+        return
+"""
+
+NullFactoryString = """
+#-------------------------------------------------------------------------------
+# Null EOS
+#-------------------------------------------------------------------------------
+class NullEquationOfState%(dim)id(IsothermalEquationOfState%(dim)id):
+    def __init__(self,
+                 constants = "PhysicalConstants&",
+                 minimumPressure = 0.0,
+                 maximumPressure =  0.0):
+        IsothermalEquationOfState%(dim)id.__init__(self,
+                                                   1.0e-80,
+                                                   1.0,
+                                                   constants,
                                                    minimumPressure,
                                                    maximumPressure)
         return
-
 """
 
 #-------------------------------------------------------------------------------
 # Create the different instantiations.
 #-------------------------------------------------------------------------------
 for dim in dims:
+    exec(NullFactoryString % {"dim" : dim})
     for units in ("MKS", "CGS", "Cosmological", "Solar"):
         exec(EOSFactoryString % {"dim"   : dim,
                                  "units" : units})

src/Material/PhysicalConstants.cc

@@ -34,7 +34,8 @@ PhysicalConstants(const double unitLm,
   UnitMassDensity(unitMkg/(unitLm*unitLm*unitLm)),
   Sigma(StefanBoltzmannMKS*unitTeK*unitTeK*unitTeK*unitTeK/unitMkg*unitTsec*unitTsec*unitTsec),
   BlackBody(4*StefanBoltzmannMKS*unitTeK*unitTeK*unitTeK*unitTeK/cMKS*unitTsec*unitTsec*unitLm/unitMkg),
-  Planck(PlanckMKS*unitTsec/(unitMkg*unitLm*unitLm)) {
+  Planck(PlanckMKS*unitTsec/(unitMkg*unitLm*unitLm)),
+  Mu0(Mu0MKS*unitCcou*unitCcou/(unitMkg*unitLm)) {
 }
 
 //------------------------------------------------------------------------------
@@ -52,5 +53,6 @@ const double PhysicalConstants::RgasMKS =   8.314462618;//8.3144621       // J/m
 const double PhysicalConstants::NAvogadro = 6.02214076e23;//6.02214129   // mol^-1
 const double PhysicalConstants::StefanBoltzmannMKS = 5.670374419e-8;//5.67   // W/m^2/K^4
 const double PhysicalConstants::PlanckMKS = 6.62607015e-34; // J*s
+const double PhysicalConstants::Mu0MKS = 1.25663706127e-6; // N*A^-2
 
 }

src/Material/PhysicalConstants.hh

@@ -43,6 +43,7 @@ public:
   double blackBodyConstant() const;
   double planckConstant() const;
   double unitEnergyJ() const;
+  double vacuumPermeability() const;
 
 private:
   //--------------------------- Private Interface ---------------------------//
@@ -63,6 +64,7 @@ private:
   const double Sigma;
   const double BlackBody;
   const double Planck;
+  const double Mu0;
 
   // The reference MKS data we base our values on.
   static const double mpMKS;
@@ -75,6 +77,7 @@ private:
   static const double NAvogadro;
   static const double StefanBoltzmannMKS;
   static const double PlanckMKS;
+  static const double Mu0MKS;
 
 };
 

src/Material/PhysicalConstantsInline.hh

@@ -171,4 +171,13 @@ PhysicalConstants::planckConstant() const {
     return Planck;
 }
 
+//------------------------------------------------------------------------------
+// Vacuum permeability (magnetic constant)
+//------------------------------------------------------------------------------
+inline
+double
+PhysicalConstants::vacuumPermeability() const {
+    return Mu0;
+}
+
 }

src/Material/UnitConversion.cc

@@ -0,0 +1,53 @@
+//---------------------------------Spheral++----------------------------------//
+// UnitConversion
+//
+// Convert between two unit systems
+//----------------------------------------------------------------------------//
+
+#include "UnitConversion.hh"
+
+namespace Spheral {
+
+UnitConversion::
+UnitConversion(const PhysicalConstants& unitsFrom,
+               const PhysicalConstants& unitsTo):
+  mUnitsFrom(unitsFrom),
+  mUnitsTo(unitsTo) {
+  // From these units (comments to right of code for common-sense check)
+  const auto l0 = unitsFrom.unitLengthMeters(); // 0.01 (cm)
+  const auto m0 = unitsFrom.unitMassKg();       // 0.001 (g)
+  const auto t0 = unitsFrom.unitTimeSec();
+  const auto k0 = unitsFrom.unitTemperatureKelvin();
+  const auto c0 = unitsFrom.unitChargeCoulomb();
+
+  // To these units
+  const auto l = unitsTo.unitLengthMeters(); // 1000 (km)
+  const auto m = unitsTo.unitMassKg();       // 1.0 (kg)
+  const auto t = unitsTo.unitTimeSec();
+  const auto k = unitsTo.unitTemperatureKelvin();
+  const auto c = unitsTo.unitChargeCoulomb();
+
+  // Conversion factors
+  const double lC = l0 / l; // 0.01/1000=0.00001 cm/km
+  const double mC = m0 / m; // 0.001/1.0=0.001 g/kg
+  const double tC = t0 / t;
+  const double kC = k0 / k;
+  const double cC = c0 / c;
+
+  // Store these
+  mLength = lC; // From cm to km: multiply by 0.00001
+  mMass = mC;   // From g to kg:  multiply by 0.001
+  mTime = tC;
+  mTemperature = kC;
+  mCharge = cC;
+
+  // Computed factors
+  mMassDensity = mC / (lC * lC * lC); // 1e-3/(1e-5)^3=10^12 g/cm^3 / (kg/km^3)
+  mSpecificEnergy = (lC * lC) / (tC * tC);
+  mEnergyDensity = mC / (lC * tC * tC);
+  mSpecificHeat = (lC * lC) / (tC * tC * kC);
+  mSpeed = lC / tC;
+  mVolume = 1.0 / (lC * lC * lC);
+}
+
+}

src/Material/UnitConversion.hh

@@ -0,0 +1,52 @@
+//---------------------------------Spheral++----------------------------------//
+// UnitConversion
+//
+// Convert between two unit systems
+//----------------------------------------------------------------------------//
+
+#ifndef __Spheral_UnitConversion_hh__
+#define __Spheral_UnitConversion_hh__
+
+#include "PhysicalConstants.hh"
+
+namespace Spheral {
+
+class UnitConversion {
+public:
+
+  // Constructor
+  UnitConversion(const PhysicalConstants& unitsFrom,
+                 const PhysicalConstants& unitsTo);
+
+  // Fundamental quantities
+  double length() const { return mLength; }
+  double mass() const { return mMass; }
+  double time() const { return mTime; }
+  double temperature() const { return mTemperature; }
+  double charge() const { return mCharge; }
+
+  // Derived quantities
+  double massDensity() const { return mMassDensity; }
+  double specificEnergy() const { return mSpecificEnergy; }
+  double energyDensity() const { return mEnergyDensity; }
+  double pressure() const { return mEnergyDensity; }
+  double bulkModulus() const { return mEnergyDensity; }
+  double specificHeat() const { return mSpecificHeat; }
+  double entropy() const { return mSpecificHeat; }
+  double speed() const { return mSpeed; }
+  double volume() const { return mVolume; }
+
+private:
+  // Input
+  PhysicalConstants mUnitsFrom;
+  PhysicalConstants mUnitsTo;
+
+  // Precomputed unit data
+  double mLength, mMass, mTime, mTemperature, mCharge;
+  double mMassDensity, mSpecificEnergy, mEnergyDensity, mSpecificHeat, mSpeed, mVolume;
+};
+
+}
+
+#endif
+

src/PYB11/Material/Material_PYB11.py

@@ -13,6 +13,7 @@
 # Includes
 #-------------------------------------------------------------------------------
 PYB11includes += ['"Material/PhysicalConstants.hh"',
+                  '"Material/UnitConversion.hh"',
                   '"Material/EquationOfState.hh"',
                   '"Material/GammaLawGas.hh"',
                   '"Material/PolytropicEquationOfState.hh"',
@@ -35,6 +36,7 @@
 # Instantiate our types
 #-------------------------------------------------------------------------------
 from PhysicalConstants import *
+from UnitConversion import *
 from EquationOfState import *
 from GammaLawGas import *
 from PolytropicEquationOfState import *

src/PYB11/Material/PhysicalConstants.py

@@ -51,3 +51,5 @@ def pyinit(self,
                                    doc="h: the Planck constant")
     unitEnergyJ = PYB11property("double", "unitEnergyJ",
                                 doc="unit of energy in SI")
+    vacuumPermeability = PYB11property("double", "vacuumPermeability",
+                                       doc="mu_0: the vacuum permeability (magnetic constant)")