Merge branch 'main' into bgclayton_simple_macaw

Author: Yurlungur

CHANGELOG.md

@@ -3,15 +3,17 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
-- [[PR556](https://github.com/lanl/singularity-eos/pull/556) Add introspection into types available in the variant
+- [[PR556]](https://github.com/lanl/singularity-eos/pull/556) Add introspection into types available in the variant
 
 ### Fixed (Repair bugs, etc)
+- [[PR561]](https://github.com/lanl/singularity-eos/pull/561) Fix logic for kokkos-kernels in spackage so that it is only required for closure models on GPU
+- [[PR563]](https://github.com/lanl/singularity-eos/pull/563) Fixed DensityFromPressureTemperature for the Carnahan-Starling EOS.
 
 ### Changed (changing behavior/API/variables/...)
 
 ### Infrastructure (changes irrelevant to downstream codes)
-- [[559]](https://github.com/lanl/singularity-eos/pull/559) Document the intent of the virtual keyword in solvers
-- [[558]](https://github.com/lanl/singularity-eos/pull/558) Make EOSPAC CMake options depend on SINGULARITY_USE_EOSPAC option
+- [[PR559]](https://github.com/lanl/singularity-eos/pull/559) Document the intent of the virtual keyword in solvers
+- [[PR558]](https://github.com/lanl/singularity-eos/pull/558) Make EOSPAC CMake options depend on SINGULARITY_USE_EOSPAC option
 
 ### Removed (removing behavior/API/varaibles/...)
 
@@ -61,7 +63,7 @@ Date: 4/7/2025
 - [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS
 
 ### Fixed (Repair bugs, etc)
-- [PR492](https://github.com/lanl/singularity-eos/pull/492) Fix import issue in spack's package.py.
+- [[PR492]](https://github.com/lanl/singularity-eos/pull/492) Fix import issue in spack's package.py.
 - [[PR485]](https://github.com/lanl/singularity-eos/pull/485) Fix segfault in Fortran interface related to EOSPAC initialization
 - [[PR478]](https://github.com/lanl/singularity-eos/pull/478) Fix bug in KPT test. Add more extensive clang build to github CI matrix.
 - [[PR473]](https://github.com/lanl/singularity-eos/pull/473) Resolve memory issue. Thanks for the catch, Richard!

singularity-eos/eos/eos_carnahan_starling.hpp

@@ -71,9 +71,9 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
     return std::max(robust::SMALL(), (sie - _qq) / _Cv);
   }
   PORTABLE_INLINE_FUNCTION void CheckParams() const {
-    PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive");
-    PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive");
-    PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive");
+    PORTABLE_ALWAYS_REQUIRE(_Cv > 0, "Heat capacity must be strictly positive");
+    PORTABLE_ALWAYS_REQUIRE(_gm1 > 0, "Gruneisen parameter must be positive");
+    PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be strictly non-negative");
     _AZbar.CheckParams();
   }
   template <typename Indexer_t = Real *>
@@ -94,24 +94,51 @@ class CarnahanStarling : public EosBase<CarnahanStarling> {
   PORTABLE_INLINE_FUNCTION Real DensityFromPressureTemperature(
       const Real press, const Real temperature, const Real guess = robust::SMALL(),
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    Real real_root;
-    auto poly = [=](Real dens) {
-      return _Cv * temperature * _gm1 * dens * ZedFromDensity(dens);
+    static constexpr Real rho_low = 0.;
+    // Need to identify an appropriate density scaling for the tolerances.
+    // Note: cannot use `_rho0` since `_rho0` is computed from this function itself.
+    const Real xtol = 1.0e-9;
+    const Real ytol = 1.0e-9;
+
+    // At absolute zero, P = 0, so assumed vacuum state (no cold curve exists).
+    if (temperature <= std::numeric_limits<Real>::min()) {
+      return rho_low;
+    }
+
+    // Use ideal gas value when _bb is zero
+    if (_bb <= std::numeric_limits<Real>::min()) {
+      return robust::ratio(press, _Cv * _gm1 * temperature);
+    }
+    const Real rho_high = 1.0 / _bb;
+
+    // Setup lambda function for finding rho.
+    // The equation has been rewritten to avoid division by zero (polynomial equation).
+    auto f = [=](const Real x /* density */) {
+      const Real eta = _bb * x;
+      const Real term1 = x * (1. + eta + eta * eta - eta * eta * eta);
+      const Real term2 =
+          robust::ratio(press * math_utils::pow<3>(1. - eta), _Cv * _gm1 * temperature);
+      return term1 - term2;
     };
-    using RootFinding1D::findRoot;
-    using RootFinding1D::Status;
-    static constexpr Real xtol = 1.0e-12;
-    static constexpr Real ytol = 1.0e-12;
-    static constexpr Real lo_bound = robust::SMALL();
-    const Real hi_bound = robust::ratio(1.0, _bb);
-    auto status = findRoot(poly, press, guess, lo_bound, hi_bound, xtol, ytol, real_root);
-    if (status != Status::SUCCESS) {
+
+    const RootFinding1D::RootCounts root_info;
+    Real rho; // `rho` will be the root of the equation $f=0$.
+
+    /* The regula falsi method should always return a unique real root in the interval (0,
+     * b^{-1}) since f(0) < 0 and f(b^{-1}) > 0. It can be shown that the derivative of f
+     * (as a function of \eta) is always positive since \eta \in (0,1) hence the root is
+     * unique. */
+    auto status =
+        regula_falsi(f, 0.0 /*target*/, guess, rho_low /*left bracket*/,
+                     rho_high /*right bracket*/, xtol, ytol, rho, &root_info, true);
+
+    if (status != RootFinding1D::Status::SUCCESS) {
       // Root finder failed even though the solution was bracketed... this is an error
       EOS_ERROR("*** (Warning) DensityFromPressureTemperature :: Convergence not met in "
-                "Carnahan-Starling EoS (root finder util) ***\n");
-      real_root = -1.0;
+                "Carnahan-Starling EoS (root finder `regula_falsi`) ***\n");
+      rho = -1.0; // guarantees zero density is returned on convergence failure.
     }
-    return std::max(robust::SMALL(), real_root);
+    return std::max(rho_low, rho);
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(

spack-repo/v1/packages/singularity-eos/package.py

@@ -57,7 +57,7 @@ class SingularityEos(CMakePackage, CudaPackage, ROCmPackage):
     # build with kokkos, kokkos-kernels for offloading support
     variant("kokkos", default=False, description="Enable kokkos")
     variant(
-        "kokkos-kernels", default=False, description="Enable kokkos-kernals for linear algebra"
+        "kokkos-kernels", default=False, description="Enable kokkos-kernels for linear algebra"
     )
 
     # for compatibility with downstream projects
@@ -120,10 +120,13 @@ class SingularityEos(CMakePackage, CudaPackage, ROCmPackage):
     depends_on("catch2@3.0.1:", when="@1.9.0: +tests")
     depends_on("py-numpy", when="+python+tests")
 
-    # linear algebra when not using GPUs
-    depends_on("eigen@3.3.8:", when="~kokkos-kernels")
-    requires("+kokkos-kernels", when="+cuda")
-    requires("+kokkos-kernels", when="+rocm")
+    # Require kokkos for device/offloading support
+    requires("+kokkos", when="+cuda")
+    requires("+kokkos", when="+rocm")
+
+    # linear algebra when using closure models. Eigen without kokkos
+    depends_on("eigen@3.3.8:", when="~kokkos-kernels+closure")
+    requires("+kokkos-kernels", when="+kokkos+closure")
 
     depends_on("eospac", when="+eospac")
 

spack-repo/v2/spack_repo/singularity_eos/packages/singularity_eos/package.py

@@ -60,7 +60,7 @@ class SingularityEos(CMakePackage, CudaPackage, ROCmPackage):
     # build with kokkos, kokkos-kernels for offloading support
     variant("kokkos", default=False, description="Enable kokkos")
     variant(
-        "kokkos-kernels", default=False, description="Enable kokkos-kernals for linear algebra"
+        "kokkos-kernels", default=False, description="Enable kokkos-kernels for linear algebra"
     )
 
     # for compatibility with downstream projects
@@ -123,10 +123,13 @@ class SingularityEos(CMakePackage, CudaPackage, ROCmPackage):
     depends_on("catch2@3.0.1:", when="@1.9.0: +tests")
     depends_on("py-numpy", when="+python+tests")
 
-    # linear algebra when not using GPUs
-    depends_on("eigen@3.3.8:", when="~kokkos-kernels")
-    requires("+kokkos-kernels", when="+cuda")
-    requires("+kokkos-kernels", when="+rocm")
+    # Require kokkos for device/offloading support
+    requires("+kokkos", when="+cuda")
+    requires("+kokkos", when="+rocm")
+
+    # linear algebra when using closure models. Eigen without kokkos
+    depends_on("eigen@3.3.8:", when="~kokkos-kernels+closure")
+    requires("+kokkos-kernels", when="+kokkos+closure")
 
     depends_on("eospac", when="+eospac")