diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 34123969f3..8a77dd7792 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -10,6 +10,7 @@ variables:
   SINGULARITY_EOS_GCC_VERSION: "9.3.0"
   SINGULARITY_EOS_OPENMPI_VERSION: "4.1.3"
   SINGULARITY_EOS_CUDA_VERSION: "11.6.0"
+  SINGULARITY_EOS_BUILD_TYPE: "RelWithDebInfo"
   SINGULARITY_EOS_SPACK_SPEC: "singularity-eos@main+mpi+hdf5+python+tests%gcc@${SINGULARITY_EOS_GCC_VERSION} ^openmpi@${SINGULARITY_EOS_OPENMPI_VERSION}"
   COLOR_CYAN: '\033[1;36m'
   COLOR_PLAIN: '\033[0m'
@@ -45,10 +46,13 @@ before_script:
 ##########################
 
 .buildtype_debug: &buildtype_debug
-  BUILD_TYPE: "Debug"
+  SINGULARITY_EOS_BUILD_TYPE: "Debug"
 
 .buildtype_release: &buildtype_release
-  BUILD_TYPE: "Release"
+  SINGULARITY_EOS_BUILD_TYPE: "Release"
+  
+.buildtype_relDebugInfo: &buildtype_relDebugInfo
+  SINGULARITY_EOS_BUILD_TYPE: "RelWithDebInfo"
 
 .gpu: &gpu
   SINGULARITY_USE_CUDA: "true"
@@ -133,6 +137,7 @@ before_script:
     - cd build
     - |
       cmake -DCMAKE_INSTALL_PREFIX=${CI_PROJECT_DIR}/install ${SINGULARITY_CMAKE_CXX_COMPILER_FLAG} \
+            -DCMAKE_BUILD_TYPE=${SINGULARITY_EOS_BUILD_TYPE} \
             -DSINGULARITY_USE_EOSPAC=ON \
             -DSINGULARITY_USE_SPINER=ON \
             -DSINGULARITY_USE_SPINER_WITH_HDF5=ON \
@@ -190,6 +195,7 @@ test_gnu_power9_gpu:
   extends: .test
   variables:    
     <<: *power9
+    <<: *buildtype_release
     <<: *gpu
 
 test_x86_volta_gpu:
@@ -197,6 +203,7 @@ test_x86_volta_gpu:
   extends: .test
   variables:    
     <<: *x86volta
+    <<: *buildtype_release
     <<: *gpu
 
 install_gnu_skylake_fort:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1580521083..ba1c4ac246 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -43,6 +43,7 @@
 
 ### Infrastructure (changes irrelevant to downstream codes)
 - [[PR190]](https://github.com/lanl/singularity-eos/pull/190) update CI on re-git
+- [[PR245]](https://github.com/lanl/singularity-eos/pull/245) Separating get_sg_eos to other files. Build/compilation improvements, warning fixes/suppression.
 
 ### Removed (removing behavior/API/varaibles/...)
 
diff --git a/cmake/submodule_configs.cmake b/cmake/submodule_configs.cmake
index 303df4414b..64d4502b71 100644
--- a/cmake/submodule_configs.cmake
+++ b/cmake/submodule_configs.cmake
@@ -23,7 +23,6 @@ macro(singularity_cmake_config pkg)
     if(SINGULARITY_USE_CUDA)
       set(Kokkos_ENABLE_CUDA ON CACHE BOOL "" FORCE)
       set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "" FORCE)
-      set(Kokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE ON CACHE BOOL "" FORCE)
     endif()
   elseif(${pkg} STREQUAL "KokkosKernels")
     # Disable TPLs
diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt
index 285545c10e..f243d057f7 100644
--- a/singularity-eos/CMakeLists.txt
+++ b/singularity-eos/CMakeLists.txt
@@ -13,7 +13,6 @@
 #------------------------------------------------------------------------------#
 
 set(EOS_HEADERS
-    closure/mixed_cell_models.hpp
     base/fast-math/logs.hpp
     base/robust_utils.hpp
     base/root-finding-1d/root_finding.hpp
@@ -49,12 +48,28 @@ set(EOS_SRCS
 )
 
 if (SINGULARITY_BUILD_CLOSURE)
-  list(APPEND EOS_SRCS eos/singularity_eos.cpp)
-  list(APPEND EOS_HEADERS eos/singularity_eos.hpp closure/mixed_cell_models.hpp)
+  list(APPEND EOS_HEADERS
+       closure/mixed_cell_models.hpp)
+  if (SINGULARITY_USE_FORTRAN)
+    list(APPEND EOS_SRCS 
+         eos/get_sg_eos.cpp)
+    if (SINGULARITY_USE_KOKKOS)
+      list(APPEND EOS_SRCS 
+           eos/get_sg_eos_p_t.cpp
+           eos/get_sg_eos_rho_t.cpp
+           eos/get_sg_eos_rho_p.cpp
+           eos/get_sg_eos_rho_e.cpp)
+    endif()
+    list(APPEND EOS_HEADERS
+         eos/get_sg_eos.hpp
+         eos/get_sg_eos_functors.hpp)
+  endif()
 endif()
 
 if (SINGULARITY_USE_FORTRAN)
   list(APPEND EOS_SRCS eos/singularity_eos.f90)
+  list(APPEND EOS_SRCS eos/singularity_eos.cpp)
+  list(APPEND EOS_HEADERS eos/singularity_eos.hpp)
   # would rather handle this more robustly, being sloppy for now
   list(APPEND EOS_MODS singularity_eos.mod singularity_eos_types.mod)
 endif()
diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp
index 863fc355a6..2d1aa9f7de 100644
--- a/singularity-eos/base/variadic_utils.hpp
+++ b/singularity-eos/base/variadic_utils.hpp
@@ -68,10 +68,6 @@ struct flatten<type_list<Ts...>> {
 };
 
 // filter nested variadic templates
-template <template <typename> class FIRST, template <typename> class... REST>
-constexpr auto remove_first(adapt_list<FIRST, REST...>) {
-  return adapt_list<REST...>{};
-}
 
 // nested specialization filter type
 template <template <typename> class Template, typename T>
@@ -109,18 +105,13 @@ struct filter<ADAPTER, Pred, Variadic, T, Ts...> {
       typename filter<ADAPTER, Pred, Variadic, Ts...>::type>::type;
 };
 
-template <template <typename> class U, typename... Ts>
-constexpr auto filter_nested(type_list<Ts...>) {
-  using f_list = typename filter<U, is_not_duplicate_nested, type_list, Ts...>::type;
-  return f_list{};
-}
-
 template <template <typename> class FIRST, template <typename> class... Us,
           typename... Ts>
 constexpr auto filter_nested_variadic(adapt_list<FIRST, Us...> m, type_list<Ts...> l) {
-  using t1 = type_list<decltype(filter_nested<FIRST>(l))>;
+  using t1 =
+      type_list<typename filter<FIRST, is_not_duplicate_nested, type_list, Ts...>::type>;
   constexpr typename flatten<t1>::type l1{};
-  return filter_nested_variadic(remove_first(m), l1);
+  return filter_nested_variadic(adapt_list<Us...>{}, l1);
 }
 
 template <typename... Ts>
@@ -130,14 +121,14 @@ constexpr auto filter_nested_variadic(adapt_list<>, type_list<Ts...> l) {
 
 // apply class template to typelist
 template <template <typename> class T, typename... Us>
-constexpr auto transform_list(type_list<Us...>) {
-  return type_list<T<Us>...>{};
-}
+struct transform_list_struct {
+  using type = type_list<T<Us>...>;
+};
 
 template <template <typename> class... Ts, typename... Us>
 constexpr auto transform_variadic_list(type_list<Us...> list,
                                        adapt_list<Ts...> mod_list) {
-  using t1 = type_list<decltype(transform_list<Ts>(list))...>;
+  using t1 = type_list<typename transform_list_struct<Ts, Us...>::type...>;
   constexpr typename flatten<t1>::type l1{};
   constexpr auto l2 = filter_nested_variadic(mod_list, l1);
   return l2;
diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp
index c9893d1873..76aeb2d937 100644
--- a/singularity-eos/eos/eos.hpp
+++ b/singularity-eos/eos/eos.hpp
@@ -58,7 +58,7 @@ namespace singularity {
 template <typename... Ts>
 using tl = singularity::detail::type_list<Ts...>;
 
-template <template <typename> typename... Ts>
+template <template <typename> class... Ts>
 using al = singularity::detail::adapt_list<Ts...>;
 
 // transform variadic list: applies modifiers to eos's
@@ -113,8 +113,8 @@ static constexpr const auto combined_list_1 = singularity::detail::concat(
 static constexpr const auto ramped_all =
     transform_variadic_list(combined_list_1, al<BilinearRampEOS>{});
 // final combined list
-static constexpr const auto combined_list = singularity::detail::concat(
-    full_eos_list, shifted, scaled, scaled_of_shifted, unit_or_rel, ramped_all);
+static constexpr const auto combined_list =
+    singularity::detail::concat(combined_list_1, ramped_all);
 // a function that returns a Variant from a typelist
 template <typename... Ts>
 struct tl_to_Variant_struct {
diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
index 0eada63c7e..3b4174e1ed 100644
--- a/singularity-eos/eos/eos_base.hpp
+++ b/singularity-eos/eos/eos_base.hpp
@@ -80,7 +80,10 @@ char *StrCat(char *destination, const char *source) {
   using EosBase<EOSDERIVED>::FillEos;                                                    \
   using EosBase<EOSDERIVED>::EntropyFromDensityTemperature;                              \
   using EosBase<EOSDERIVED>::EntropyFromDensityInternalEnergy;                           \
-  using EosBase<EOSDERIVED>::EntropyIsNotEnabled;
+  using EosBase<EOSDERIVED>::EntropyIsNotEnabled;                                        \
+  using EosBase<EOSDERIVED>::IsModified;                                                 \
+  using EosBase<EOSDERIVED>::UnmodifyOnce;                                               \
+  using EosBase<EOSDERIVED>::GetUnmodifiedObject;
 
 class Factor {
   Real value_ = 1.0;
@@ -606,12 +609,13 @@ class EosBase {
   }
 
   // Tooling for modifiers
-  PORTABLE_FORCEINLINE_FUNCTION
-  bool IsModified() const { return false; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  auto UnmodifyOnce() { return *static_cast<CRTP *>(this); }
-  PORTABLE_FORCEINLINE_FUNCTION
-  auto GetUnmodifiedObject() { return *static_cast<CRTP *>(this); }
+  inline constexpr bool IsModified() const { return false; }
+
+  inline constexpr decltype(auto) UnmodifyOnce() { return *static_cast<CRTP *>(this); }
+
+  inline constexpr decltype(auto) GetUnmodifiedObject() {
+    return *static_cast<CRTP *>(this);
+  }
 };
 } // namespace eos_base
 } // namespace singularity
diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
index 7de2bf282e..c8c59de598 100644
--- a/singularity-eos/eos/eos_eospac.hpp
+++ b/singularity-eos/eos/eos_eospac.hpp
@@ -41,41 +41,64 @@ using namespace eos_base;
 
 // Only really works in serial
 // Not really supported on device
+#ifndef PORTABILITY_STRATEGY_NONE
+#if defined(__CUDACC__)
+#define SG_PIF_NOWARN #pragma nv_exec_check_disable
+#endif // __CUDACC__
+#else
+#define SG_PIF_NOWARN
+#endif // PORTABILITY_STRATEGY_NONE
+
 class EOSPAC : public EosBase<EOSPAC> {
  public:
   inline EOSPAC() = default;
   inline EOSPAC(int matid, bool invert_at_setup = false);
   inline EOSPAC GetOnDevice() { return *this; }
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
       const Real rho, const Real sie, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy(
       const Real rho, const Real sie, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy(
       const Real rho, const Real sie, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(
       const Real rho, const Real sie, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy(
       const Real rho, const Real sie, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature(
       const Real rho, const Real temperature, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(
       const Real rho, const Real sie, Real *lambda = nullptr) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                         Real &cv, Real &bmod, const unsigned long output,
                                         Real *lambda = nullptr) const;
-  PORTABLE_INLINE_FUNCTION
-  void DensityEnergyFromPressureTemperature(const Real press, const Real temp,
-                                            Real *lambda, Real &rho, Real &sie) const;
+  SG_PIF_NOWARN
+  PORTABLE_INLINE_FUNCTION void
+  DensityEnergyFromPressureTemperature(const Real press, const Real temp, Real *lambda,
+                                       Real &rho, Real &sie) const;
+  SG_PIF_NOWARN
   PORTABLE_INLINE_FUNCTION void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie,
                                                        Real &press, Real &cv, Real &bmod,
                                                        Real &dpde, Real &dvdt,
@@ -1094,9 +1117,10 @@ class EOSPAC : public EosBase<EOSPAC> {
         });
   }
 
-  PORTABLE_INLINE_FUNCTION
-  void PTofRE(Real &rho, Real &sie, Real *lambda, Real &press, Real &temp, Real &dpdr,
-              Real &dpde, Real &dtdr, Real &dtde) const {
+  SG_PIF_NOWARN
+  PORTABLE_INLINE_FUNCTION void PTofRE(Real &rho, Real &sie, Real *lambda, Real &press,
+                                       Real &temp, Real &dpdr, Real &dpde, Real &dtdr,
+                                       Real &dtde) const {
 
     press = PressureFromDensityInternalEnergy(rho, sie, lambda);
     temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda);
@@ -1138,7 +1162,9 @@ class EOSPAC : public EosBase<EOSPAC> {
   }
   static std::string EosType() { return std::string("EOSPAC"); }
   static std::string EosPyType() { return EosType(); }
-  inline void PrintParams() const { printf("EOSPAC parameters:\nmatid = %s\n", matid_); }
+  PORTABLE_INLINE_FUNCTION void PrintParams() const {
+    printf("EOSPAC parameters:\nmatid = %i\n", matid_);
+  }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; }
 
@@ -1232,6 +1258,7 @@ inline EOSPAC::EOSPAC(const int matid, bool invert_at_setup) : matid_(matid) {
       robust::ratio(dpde_ref_ * cv_ref_, rho_ref_ * rho_ref_ * pressureFromSesame(DPDR));
 }
 
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::TemperatureFromDensityInternalEnergy(
     const Real rho, const Real sie, Real *lambda) const {
   using namespace EospacWrapper;
@@ -1242,6 +1269,7 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::TemperatureFromDensityInternalEnergy(
   return Real(temperatureFromSesame(T[0]));
 }
 
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityTemperature(const Real rho,
                                                                      const Real temp,
                                                                      Real *lambda) const {
@@ -1253,12 +1281,14 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityTemperature(const Real
   return Real(pressureFromSesame(P[0]));
 }
 
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::EntropyFromDensityTemperature(
     const Real rho, const Real temperature, Real *lambda) const {
   EntropyIsNotEnabled("EOSPAC");
   return 1.0;
 }
 
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION void EOSPAC::FillEos(Real &rho, Real &temp, Real &sie,
                                               Real &press, Real &cv, Real &bmod,
                                               const unsigned long output,
@@ -1268,7 +1298,7 @@ PORTABLE_INLINE_FUNCTION void EOSPAC::FillEos(Real &rho, Real &temp, Real &sie,
   EOS_REAL E[1] = {sie}, P[1] = {pressureToSesame(press)};
   EOS_REAL dx[1], dy[1];
   EOS_INTEGER nxypairs = 1;
-  Real CV, BMOD_T, BMOD, SIE, PRESS, DPDE, DPDT, DPDR, DEDT, DEDR;
+  Real /*CV,*/ BMOD_T, BMOD, SIE, PRESS, DPDE, DPDT, DPDR, DEDT, DEDR;
   const unsigned long input = ~output;
   if (output == thermalqs::none) {
     UNDEFINED_ERROR;
@@ -1342,6 +1372,7 @@ PORTABLE_INLINE_FUNCTION void EOSPAC::FillEos(Real &rho, Real &temp, Real &sie,
   }
 }
 
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::InternalEnergyFromDensityTemperature(
     const Real rho, const Real temp, Real *lambda) const {
   using namespace EospacWrapper;
@@ -1350,6 +1381,7 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::InternalEnergyFromDensityTemperature(
   FillEos(RHO, TEMP, sie, press, cv, bmod, output, lambda);
   return sie;
 }
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::BulkModulusFromDensityTemperature(
     const Real rho, const Real temp, Real *lambda) const {
   using namespace EospacWrapper;
@@ -1358,6 +1390,7 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::BulkModulusFromDensityTemperature(
   FillEos(RHO, TEMP, sie, press, cv, bmod, output, lambda);
   return bmod;
 }
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::SpecificHeatFromDensityTemperature(
     const Real rho, const Real temp, Real *lambda) const {
   using namespace EospacWrapper;
@@ -1366,24 +1399,28 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::SpecificHeatFromDensityTemperature(
   FillEos(RHO, TEMP, sie, press, cv, bmod, output, lambda);
   return cv;
 }
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::PressureFromDensityInternalEnergy(
     const Real rho, const Real sie, Real *lambda) const {
   using namespace EospacWrapper;
   Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda);
   return PressureFromDensityTemperature(rho, temp, lambda);
 }
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::EntropyFromDensityInternalEnergy(
     const Real rho, const Real sie, Real *lambda) const {
   using namespace EospacWrapper;
   const Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda);
   return EntropyFromDensityTemperature(rho, temp, lambda);
 }
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::SpecificHeatFromDensityInternalEnergy(
     const Real rho, const Real sie, Real *lambda) const {
   using namespace EospacWrapper;
   Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda);
   return SpecificHeatFromDensityTemperature(rho, temp, lambda);
 }
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::BulkModulusFromDensityInternalEnergy(
     const Real rho, const Real sie, Real *lambda) const {
   using namespace EospacWrapper;
@@ -1391,6 +1428,7 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::BulkModulusFromDensityInternalEnergy(
   return BulkModulusFromDensityTemperature(rho, temp, lambda);
 }
 
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::GruneisenParamFromDensityTemperature(
     const Real rho, const Real temperature, Real *lambda) const {
   using namespace EospacWrapper;
@@ -1407,22 +1445,22 @@ PORTABLE_INLINE_FUNCTION Real EOSPAC::GruneisenParamFromDensityTemperature(
   DPDE = DPDT / DEDT;
   return robust::ratio(pressureFromSesame(sieToSesame(DPDE)), rho);
 }
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION Real EOSPAC::GruneisenParamFromDensityInternalEnergy(
     const Real rho, const Real sie, Real *lambda) const {
   Real temperature = TemperatureFromDensityInternalEnergy(rho, sie, lambda);
   return GruneisenParamFromDensityTemperature(rho, temperature, lambda);
 }
 
-PORTABLE_INLINE_FUNCTION
-void EOSPAC::DensityEnergyFromPressureTemperature(const Real press, const Real temp,
-                                                  Real *lambda, Real &rho,
-                                                  Real &sie) const {
+SG_PIF_NOWARN
+PORTABLE_INLINE_FUNCTION void
+EOSPAC::DensityEnergyFromPressureTemperature(const Real press, const Real temp,
+                                             Real *lambda, Real &rho, Real &sie) const {
   using namespace EospacWrapper;
   EOS_REAL P[1] = {pressureToSesame(press)};
   EOS_REAL T[1] = {temperatureToSesame(temp)};
   EOS_REAL dx[1], dy[1], R[1], E[1];
   EOS_INTEGER nxypairs = 1;
-  EOS_INTEGER errorCode;
   EOS_INTEGER table;
 
   table = RofPT_table_;
@@ -1434,6 +1472,7 @@ void EOSPAC::DensityEnergyFromPressureTemperature(const Real press, const Real t
   sie = sieFromSesame(E[0]);
 }
 
+SG_PIF_NOWARN
 PORTABLE_INLINE_FUNCTION void
 EOSPAC::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
                                Real &bmod, Real &dpde, Real &dvdt, Real *lambda) const {
diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp
index 43b34b7181..2dddbc3bdc 100644
--- a/singularity-eos/eos/eos_spiner.hpp
+++ b/singularity-eos/eos/eos_spiner.hpp
@@ -201,25 +201,21 @@ class SpinerEOSDependsRhoT : public EosBase<SpinerEOSDependsRhoT> {
   static PORTABLE_FORCEINLINE_FUNCTION Real fromLog_(const Real lx, const Real offset) {
     return FastMath::pow10(lx) - offset;
   }
-  PORTABLE_INLINE_FUNCTION Real __attribute__((always_inline))
-  lRho_(const Real rho) const noexcept {
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real lRho_(const Real rho) const noexcept {
     Real out = toLog_(rho, lRhoOffset_);
     return out;
     // return out < lRhoMin_ ? lRhoMin_ : out;
   }
-  PORTABLE_INLINE_FUNCTION Real __attribute__((always_inline))
-  lT_(const Real T) const noexcept {
-    return toLog_(T, lTOffset_);
-  }
-  PORTABLE_INLINE_FUNCTION Real __attribute__((always_inline))
-  rho_(const Real lRho) const noexcept {
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real lT_(const Real T) const noexcept { return toLog_(T, lTOffset_); }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real rho_(const Real lRho) const noexcept {
     Real rho = fromLog_(lRho, lRhoOffset_);
     return rho < 0 ? 0 : rho;
   }
-  PORTABLE_INLINE_FUNCTION Real __attribute__((always_inline))
-  T_(const Real lT) const noexcept {
-    return fromLog_(lT, lTOffset_);
-  }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real T_(const Real lT) const noexcept { return fromLog_(lT, lTOffset_); }
 
   PORTABLE_INLINE_FUNCTION
   Real lTFromlRhoSie_(const Real lRho, const Real sie, TableStatus &whereAmI,
@@ -1204,7 +1200,7 @@ Real SpinerEOSDependsRhoT::lRhoFromPlT_(const Real P, const Real lT,
                     lRhoMinSearch_, lRhoMax_, ROOT_THRESH, ROOT_THRESH, lRho, pcounts);
   }
   if (status != RootFinding1D::Status::SUCCESS) {
-#if EPINER_EOS_VERBOSE
+#if SPINER_EOS_VERBOSE
     std::stringstream errorMessage;
     errorMessage << "inverting P table for logRho failed\n"
                  << "matid     = " << matid_ << "\n"
@@ -1219,8 +1215,10 @@ Real SpinerEOSDependsRhoT::lRhoFromPlT_(const Real P, const Real lT,
     lambda[Lambda::lRho] = lRho;
     lambda[Lambda::lT] = lT;
   }
-  status_ = status;
-  whereAmI_ = whereAmI;
+  if (memoryStatus_ != DataStatus::OnDevice) {
+    status_ = status;
+    whereAmI_ = whereAmI;
+  }
   return lRho;
 }
 
@@ -1281,8 +1279,8 @@ Real SpinerEOSDependsRhoT::lTFromlRhoSie_(const Real lRho, const Real sie,
   }
   if (memoryStatus_ != DataStatus::OnDevice) {
     status_ = status;
+    whereAmI_ = whereAmI;
   }
-  whereAmI_ = whereAmI;
   return lT;
 }
 
@@ -1339,8 +1337,8 @@ Real SpinerEOSDependsRhoT::lTFromlRhoP_(const Real lRho, const Real press,
   }
   if (memoryStatus_ != DataStatus::OnDevice) {
     status_ = status;
+    whereAmI_ = whereAmI;
   }
-  whereAmI_ = whereAmI;
   return lT;
 }
 
@@ -1414,8 +1412,9 @@ Real SpinerEOSDependsRhoT::bModFromRholRhoTlT_(const Real rho, const Real lRho,
   return bMod > robust::EPS() ? bMod : robust::EPS();
 }
 
-PORTABLE_INLINE_FUNCTION TableStatus
-SpinerEOSDependsRhoT::getLocDependsRhoSie_(const Real lRho, const Real sie) const {
+PORTABLE_INLINE_FUNCTION
+TableStatus SpinerEOSDependsRhoT::getLocDependsRhoSie_(const Real lRho,
+                                                       const Real sie) const {
   TableStatus whereAmI;
   const Real sielTMax = sielTMax_.interpToReal(lRho);
   const Real sieCold = sieCold_.interpToReal(lRho);
@@ -1427,7 +1426,9 @@ SpinerEOSDependsRhoT::getLocDependsRhoSie_(const Real lRho, const Real sie) cons
   } else {
     whereAmI = TableStatus::OnTable;
   }
-  whereAmI_ = whereAmI;
+  if (memoryStatus_ != DataStatus::OnDevice) {
+    whereAmI_ = whereAmI;
+  }
   return whereAmI;
 }
 
@@ -1440,7 +1441,9 @@ SpinerEOSDependsRhoT::getLocDependsRhoT_(const Real lRho, const Real lT) const {
     whereAmI = TableStatus::OffTop;
   else
     whereAmI = TableStatus::OnTable;
-  whereAmI_ = whereAmI;
+  if (memoryStatus_ != DataStatus::OnDevice) {
+    whereAmI_ = whereAmI;
+  }
   return whereAmI;
 }
 
diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp
index bb6b035039..dbab82d1ce 100644
--- a/singularity-eos/eos/eos_variant.hpp
+++ b/singularity-eos/eos/eos_variant.hpp
@@ -960,17 +960,16 @@ class Variant {
   }
 
   // Tooling for modifiers
-  PORTABLE_FORCEINLINE_FUNCTION
-  bool IsModified() const {
+  inline constexpr bool IsModified() const {
     return mpark::visit([](const auto &eos) { return eos.IsModified(); }, eos_);
   }
-  PORTABLE_FORCEINLINE_FUNCTION
-  Variant UnmodifyOnce() {
+
+  inline constexpr Variant UnmodifyOnce() {
     return mpark::visit(
         [](auto &eos) { return eos_variant<EOSs...>(eos.UnmodifyOnce()); }, eos_);
   }
-  PORTABLE_FORCEINLINE_FUNCTION
-  Variant GetUnmodifiedObject() {
+
+  inline constexpr Variant GetUnmodifiedObject() {
     return mpark::visit(
         [](auto &eos) { return eos_variant<EOSs...>(eos.GetUnmodifiedObject()); }, eos_);
   }
diff --git a/singularity-eos/eos/get_sg_eos.cpp b/singularity-eos/eos/get_sg_eos.cpp
new file mode 100644
index 0000000000..7e697d56fe
--- /dev/null
+++ b/singularity-eos/eos/get_sg_eos.cpp
@@ -0,0 +1,275 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#include <map>
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/closure/mixed_cell_models.hpp>
+#include <singularity-eos/eos/eos.hpp>
+#include <singularity-eos/eos/get_sg_eos.hpp>
+#include <singularity-eos/eos/get_sg_eos_functors.hpp>
+#include <singularity-eos/eos/singularity_eos.hpp>
+
+using namespace singularity;
+
+// mapping from EAP integer to
+static const std::map<const int, const unsigned long> EAPInputToBD = {
+    {-3, thermalqs::temperature | thermalqs::density},
+    {-2, thermalqs::density | thermalqs::pressure},
+    {-1, thermalqs::pressure | thermalqs::temperature},
+    {0, thermalqs::specific_internal_energy | thermalqs::density},
+    {1, thermalqs::specific_internal_energy | thermalqs::density},
+};
+
+// EAP centric arguments and function signature
+int get_sg_eos( // sizing information
+    int nmat, int ncell, int cell_dim,
+    // Input parameters
+    int input_int,
+    // eos index offsets
+    int *eos_offsets,
+    // equation of state array
+    EOS *eos,
+    // index offsets
+    int *offsets,
+    // per cell quantities
+    double *press, double *pmax, double *vol, double *spvol, double *sie, double *temp,
+    double *bmod, double *dpde, double *cv,
+    // per material quantities
+    double *frac_mass, double *frac_vol, double *frac_ie,
+    // optional per material quantities
+    double *frac_bmod, double *frac_dpde, double *frac_cv) {
+  // printBacktrace();
+  // kernel return value will be the number of failures
+  int ret{0};
+  // handle optionals
+  bool do_frac_bmod{true};
+  if (frac_bmod == NULL || frac_bmod == nullptr) do_frac_bmod = false;
+  bool do_frac_dpde{true};
+  if (frac_dpde == NULL || frac_dpde == nullptr) do_frac_dpde = false;
+  bool do_frac_cv{true};
+  if (frac_dpde == NULL || frac_dpde == nullptr) do_frac_cv = false;
+  // get inputs
+  enum class input_condition {
+    RHO_T_INPUT = -3,
+    RHO_P_INPUT = -2,
+    P_T_INPUT = -1,
+    NORM_RHO_E_INPUT = 0,
+    RHO_E_INPUT = 1
+  };
+  const auto input{EAPInputToBD.at(input_int)};
+  const bool p_is_inp{static_cast<bool>(input & thermalqs::pressure)};
+  const bool r_is_inp{static_cast<bool>(input & thermalqs::density)};
+  const bool t_is_inp{static_cast<bool>(input & thermalqs::temperature)};
+  const bool s_is_inp{static_cast<bool>(input & thermalqs::specific_internal_energy)};
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+  // convert pointers to host side views
+  Kokkos::View<int *, Llft, HS, Unmgd> eos_offsets_hv(eos_offsets, nmat);
+  Kokkos::View<int *, Llft, HS, Unmgd> offsets_hv(offsets, ncell);
+  host_v press_hv(press, cell_dim);
+  host_v pmax_hv(pmax, cell_dim);
+  host_v vol_hv(vol, cell_dim);
+  host_v spvol_hv(spvol, cell_dim);
+  host_v sie_hv(sie, cell_dim);
+  host_v temp_hv(temp, cell_dim);
+  host_v bmod_hv(bmod, cell_dim);
+  host_v dpde_hv(dpde, cell_dim);
+  host_v cv_hv(cv, cell_dim);
+  host_frac_v frac_mass_hv(frac_mass, cell_dim, nmat);
+  host_frac_v frac_vol_hv(frac_vol, cell_dim, nmat);
+  host_frac_v frac_ie_hv(frac_ie, cell_dim, nmat);
+  host_frac_v frac_bmod_hv, frac_dpde_hv, frac_cv_hv;
+  if (do_frac_bmod) frac_bmod_hv = host_frac_v(frac_bmod, cell_dim, nmat);
+  if (do_frac_dpde) frac_dpde_hv = host_frac_v(frac_dpde, cell_dim, nmat);
+  if (do_frac_cv) frac_cv_hv = host_frac_v(frac_cv, cell_dim, nmat);
+
+  // get device views if necessary
+  indirection_v offsets_v{create_mirror_view_and_copy(DMS(), offsets_hv)};
+  indirection_v eos_offsets_v{create_mirror_view_and_copy(DMS(), eos_offsets_hv)};
+  dev_v press_v{create_mirror_view_and_copy(DMS(), press_hv)};
+  dev_v pmax_v{create_mirror_view_and_copy(DMS(), pmax_hv)};
+  dev_v spvol_v{create_mirror_view_and_copy(DMS(), spvol_hv)};
+  dev_v vol_v{create_mirror_view_and_copy(DMS(), vol_hv)};
+  dev_v sie_v{create_mirror_view_and_copy(DMS(), sie_hv)};
+  dev_v temp_v{create_mirror_view_and_copy(DMS(), temp_hv)};
+  dev_v bmod_v{create_mirror_view_and_copy(DMS(), bmod_hv)};
+  dev_v dpde_v{create_mirror_view_and_copy(DMS(), dpde_hv)};
+  dev_v cv_v{create_mirror_view_and_copy(DMS(), cv_hv)};
+  dev_frac_v frac_mass_v{create_mirror_view_and_copy(DMS(), frac_mass_hv)};
+  dev_frac_v frac_vol_v{create_mirror_view_and_copy(DMS(), frac_vol_hv)};
+  dev_frac_v frac_ie_v{create_mirror_view_and_copy(DMS(), frac_ie_hv)};
+  dev_frac_v frac_bmod_v, frac_dpde_v, frac_cv_v;
+  if (do_frac_bmod) frac_bmod_v = create_mirror_view_and_copy(DMS(), frac_bmod_hv);
+  if (do_frac_dpde) frac_dpde_v = create_mirror_view_and_copy(DMS(), frac_dpde_hv);
+  if (do_frac_cv) frac_cv_v = create_mirror_view_and_copy(DMS(), frac_cv_hv);
+  // array of eos's
+  const auto eos_nmat{*std::max_element(eos_offsets, eos_offsets + nmat)};
+  Kokkos::View<EOS *, Llft, HS, Unmgd> eos_hv(eos, eos_nmat);
+  Kokkos::View<EOS *, Llft> eos_v{create_mirror_view_and_copy(DMS(), eos_hv)};
+  // TODO: make this a scatter view
+  constexpr auto at_int_full{0 | Kokkos::Atomic};
+#ifdef KOKKOS_ENABLE_SERIAL
+  constexpr auto at_int{std::is_same<DES, Kokkos::Serial>::value ? 0 : at_int_full};
+#else  // KOKKOS_ENABLE_SERIAL
+  // assume atomics are required for correctness if serial not even enabled
+  constexpr auto at_int{at_int_full};
+#endif // KOKKOS_ENABLE_SERIAL
+  // set up scratch arrays
+  constexpr auto KGlobal = Kokkos::Experimental::UniqueTokenScope::Global;
+  Kokkos::Experimental::UniqueToken<DES, KGlobal> tokens{};
+  using VAWI = Kokkos::ViewAllocateWithoutInitializing;
+
+  const bool small_loop{tokens.size() > ncell};
+  const decltype(tokens)::size_type scratch_size{std::min(tokens.size(), ncell)};
+  ScratchV<int> pte_mats(VAWI("PTE::scratch mats"), scratch_size, nmat);
+  ScratchV<int> pte_idxs(VAWI("PTE::scratch idxs"), scratch_size, nmat);
+  ScratchV<double> mass_pte(VAWI("PTE::scratch mass"), scratch_size, nmat);
+  ScratchV<double> sie_pte(VAWI("PTE::scratch sie"), scratch_size, nmat);
+  ScratchV<double> vfrac_pte(VAWI("PTE::scratch vfrac"), scratch_size, nmat);
+  ScratchV<double> temp_pte(VAWI("PTE::scratch temp"), scratch_size, nmat);
+  ScratchV<double> press_pte(VAWI("PTE::scratch press"), scratch_size, nmat);
+  ScratchV<double> rho_pte(VAWI("PTE::scratch rho"), scratch_size, nmat);
+  // declare solver scratch, allocate in the case statement
+  int pte_solver_scratch_size{};
+  ScratchV<double> solver_scratch;
+  // declare init and final functors
+  auto input_int_enum = static_cast<input_condition>(input_int);
+  init_functor i_func;
+  final_functor f_func(temp_v, press_v, sie_v, bmod_v, cv_v, dpde_v, pte_mats, press_pte,
+                       vfrac_pte, temp_pte, sie_pte, frac_mass_v, frac_ie_v, frac_vol_v,
+                       vol_v, eos_v, pte_idxs, rho_pte, frac_bmod_v, frac_cv_v,
+                       frac_dpde_v, nmat, do_frac_bmod, do_frac_cv, do_frac_dpde);
+  // only initialize init functor when needed
+  if (input_int_enum != input_condition::P_T_INPUT) {
+    i_func = init_functor(frac_mass_v, pte_idxs, eos_offsets_v, frac_vol_v, pte_mats,
+                          vfrac_pte, sie_pte, temp_pte, press_pte, rho_pte, spvol_v,
+                          temp_v, press_v, sie_v, nmat);
+  }
+
+  // create helper lambdas to reduce code duplication
+  Kokkos::View<int, MemoryTraits<at_int>> res("PTE::num fails");
+  Kokkos::View<int, MemoryTraits<at_int>> n_solves("PTE::num solves");
+  const std::string perf_nums =
+      "[" + std::to_string(nmat) + "," + std::to_string(ncell) + "]";
+  switch (input_int_enum) {
+  case input_condition::RHO_T_INPUT: {
+    // T-rho input
+    // set frac_vol = 1/nmat
+    // set rho_i = nmat / spvol * frac_mass_i
+    // iterate PTE solver to obtain internal energies
+    // that results in the input T
+    pte_solver_scratch_size = PTESolverFixedTRequiredScratch(nmat);
+    solver_scratch = ScratchV<double>(VAWI("PTE::scratch solver"), scratch_size,
+                                      pte_solver_scratch_size);
+    const std::string rt_name = "PTE::solve (rho,T) input" + perf_nums;
+    singularity::get_sg_eos_rho_t(rt_name.c_str(), ncell, offsets_v, eos_v, press_v,
+                                  pmax_v, sie_v, frac_mass_v, pte_idxs, pte_mats,
+                                  press_pte, vfrac_pte, rho_pte, sie_pte, temp_pte,
+                                  solver_scratch, tokens, small_loop, i_func, f_func);
+    break;
+  }
+  case input_condition::RHO_P_INPUT: {
+    // rho-P input
+    // set frac_vol = 1/nmat
+    // set rho_i = nmat / spvol * frac_mass_i
+    // iterate PTE solver to obtain internal energies
+    // that results in the input P
+    pte_solver_scratch_size = PTESolverRhoTRequiredScratch(nmat);
+    solver_scratch = ScratchV<double>(VAWI("PTE::scratch solver"), scratch_size,
+                                      pte_solver_scratch_size);
+    const std::string rp_name = "PTE::solve (rho,P) input" + perf_nums;
+    singularity::get_sg_eos_rho_p(rp_name.c_str(), ncell, offsets_v, eos_v, press_v,
+                                  pmax_v, sie_v, frac_mass_v, pte_idxs, pte_mats,
+                                  press_pte, vfrac_pte, rho_pte, sie_pte, temp_pte,
+                                  solver_scratch, tokens, small_loop, i_func, f_func);
+    break;
+  }
+  case input_condition::P_T_INPUT: {
+    // P-T input
+    const int pte_solver_scratch_size = nmat * MAX_NUM_LAMBDAS;
+    solver_scratch = ScratchV<double>(VAWI("PTE::scratch solver"), scratch_size,
+                                      pte_solver_scratch_size);
+    const std::string pt_name = "PTE::solve (P,T) input" + perf_nums;
+    singularity::get_sg_eos_p_t(pt_name.c_str(), ncell, nmat, offsets_v, eos_offsets_v,
+                                eos_v, press_v, pmax_v, vol_v, spvol_v, sie_v, temp_v,
+                                frac_mass_v, pte_idxs, pte_mats, press_pte, vfrac_pte,
+                                rho_pte, sie_pte, temp_pte, solver_scratch, tokens,
+                                small_loop, f_func);
+    break;
+  }
+  case input_condition::NORM_RHO_E_INPUT:
+    // rho-sie input
+    // no break so fallthrough to case 1
+  case input_condition::RHO_E_INPUT: {
+    // rho-sie input
+    pte_solver_scratch_size = PTESolverRhoTRequiredScratch(nmat);
+    solver_scratch = ScratchV<double>(VAWI("PTE::scratch solver"), scratch_size,
+                                      pte_solver_scratch_size);
+    const std::string re_name = "PTE::solve (rho,e) input" + perf_nums;
+    singularity::get_sg_eos_rho_e(re_name.c_str(), ncell, offsets_v, eos_v, press_v,
+                                  pmax_v, sie_v, pte_idxs, press_pte, vfrac_pte, rho_pte,
+                                  sie_pte, temp_pte, solver_scratch, tokens, small_loop,
+                                  i_func, f_func);
+    break;
+  }
+  }
+
+  Kokkos::fence();
+  // copy results back into local values
+  // there is lots of room for performance optimization
+  // in terms of when to copy and when not necessary
+  // this is the return value (number of solve failures)
+  // deep_copy(ret, res);
+  // copy pressure, this is not needed in all cases
+  if (!p_is_inp) {
+    deep_copy(DES(), press_hv, press_v);
+  }
+  // return max pressure, this may be needed
+  deep_copy(DES(), pmax_hv, pmax_v);
+  // I don't think the volume is necessary
+  deep_copy(DES(), vol_hv, vol_v);
+  // specific volume, copy-back not needed in all cases
+  if (!r_is_inp) {
+    deep_copy(DES(), spvol_hv, spvol_v);
+  }
+  // internal energy, copy-back not needed in all cases
+  if (!s_is_inp) {
+    deep_copy(DES(), sie_hv, sie_v);
+  }
+  // temperature, copy-back not needed in all cases
+  if (!t_is_inp) {
+    deep_copy(DES(), temp_hv, temp_v);
+  }
+  // bulk modulus, alwasy copy-back
+  deep_copy(DES(), bmod_hv, bmod_v);
+  // dpde, always copy-back
+  deep_copy(DES(), dpde_hv, dpde_v);
+  // specific heat, always copy-back
+  deep_copy(DES(), cv_hv, cv_v);
+  // volume fractions, always copy-back (maybe not for pure cells)
+  deep_copy(DES(), frac_vol_hv, frac_vol_v);
+  // component internal energies, always copy-back (maybe not for pure cells)
+  deep_copy(DES(), frac_ie_hv, frac_ie_v);
+  // optionally copy-back the component bmod, dpde, and cv
+  if (do_frac_bmod) {
+    deep_copy(DES(), frac_bmod_hv, frac_bmod_v);
+  }
+  if (do_frac_dpde) {
+    deep_copy(DES(), frac_dpde_hv, frac_dpde_v);
+  }
+  if (do_frac_cv) {
+    deep_copy(DES(), frac_cv_hv, frac_cv_v);
+  }
+#endif // PORTABILITY_STRATEGY_KOKKOS
+  return ret;
+}
diff --git a/singularity-eos/eos/get_sg_eos.hpp b/singularity-eos/eos/get_sg_eos.hpp
new file mode 100644
index 0000000000..9baccd83e1
--- /dev/null
+++ b/singularity-eos/eos/get_sg_eos.hpp
@@ -0,0 +1,116 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#ifndef _SINGULARITY_EOS_EOS_GET_SG_EOS_HPP_
+#define _SINGULARITY_EOS_EOS_GET_SG_EOS_HPP_
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/eos/eos.hpp>
+
+using namespace singularity;
+
+// forward declare init/final functors
+namespace singularity {
+struct init_functor;
+struct final_functor;
+} // namespace singularity
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+using Lrgt = Kokkos::LayoutRight;
+using Llft = Kokkos::LayoutLeft;
+template <typename T>
+using ScratchV = Kokkos::View<T **, Lrgt>;
+#endif // PORTABILITY_STRATEGY_KOKKOS
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+namespace singularity {
+struct EOSAccessor_ {
+  PORTABLE_INLINE_FUNCTION
+  EOSAccessor_(const Kokkos::View<EOS *, Llft> &eos_v, int *mats)
+      : eos_v_(eos_v), mats_(mats) {}
+  PORTABLE_INLINE_FUNCTION
+  EOS &operator[](const int m) const { return eos_v_(mats_[m]); }
+  Kokkos::View<EOS *, Llft> eos_v_;
+  int *mats_;
+};
+} // namespace singularity
+#endif // PORTABILITY_STRATEGY_KOKKOS
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+// convenience aliases
+using Unmgd = Kokkos::MemoryUnmanaged;
+using HS = Kokkos::HostSpace;
+using Kokkos::create_mirror_view_and_copy;
+using host_v = Kokkos::View<double *, Llft, HS, Unmgd>;
+using dev_v = Kokkos::View<double *, Llft>;
+using host_frac_v = Kokkos::View<double **, Llft, HS, Unmgd>;
+using dev_frac_v = Kokkos::View<double **, Llft>;
+using DES = Kokkos::DefaultExecutionSpace;
+using DMS = DES::memory_space;
+using Kokkos::MemoryTraits;
+constexpr const unsigned int ra{0 | Kokkos::RandomAccess};
+// constexpr const unsigned int ra_u{Kokkos::Unmanaged | Kokkos::RandomAccess};
+using indirection_v = Kokkos::View<const int *, Llft, MemoryTraits<ra>>;
+using VAWI = Kokkos::ViewAllocateWithoutInitializing;
+using Kokkos::deep_copy;
+constexpr auto KGlobal = Kokkos::Experimental::UniqueTokenScope::Global;
+static constexpr const double ev2k = 1.160451930280894026e4;
+#endif // PORTABILITY_STRATEGY_KOKKOS
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+namespace singularity {
+// rho t input
+void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v,
+                      Kokkos::View<EOS *, Llft> &eos_v, dev_v &press_v, dev_v &pmax_v,
+                      dev_v &sie_v, dev_frac_v &frac_mass_v, ScratchV<int> &pte_idxs,
+                      ScratchV<int> &pte_mats, ScratchV<double> &press_pte,
+                      ScratchV<double> &vfrac_pte, ScratchV<double> &rho_pte,
+                      ScratchV<double> &sie_pte, ScratchV<double> &temp_pte,
+                      ScratchV<double> &solver_scratch,
+                      Kokkos::Experimental::UniqueToken<DES, KGlobal> &tokens,
+                      bool small_loop, init_functor &i_func, final_functor &f_func);
+// rho P input
+void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v,
+                      Kokkos::View<EOS *, Llft> &eos_v, dev_v &press_v, dev_v &pmax_v,
+                      dev_v &sie_v, dev_frac_v &frac_mass_v, ScratchV<int> &pte_idxs,
+                      ScratchV<int> &pte_mats, ScratchV<double> &press_pte,
+                      ScratchV<double> &vfrac_pte, ScratchV<double> &rho_pte,
+                      ScratchV<double> &sie_pte, ScratchV<double> &temp_pte,
+                      ScratchV<double> &solver_scratch,
+                      Kokkos::Experimental::UniqueToken<DES, KGlobal> &tokens,
+                      bool small_loop, init_functor &i_func, final_functor &f_func);
+// PT input
+void get_sg_eos_p_t(const char *name, int ncell, int nmat, indirection_v &offsets_v,
+                    indirection_v &eos_offsets_v, Kokkos::View<EOS *, Llft> &eos_v,
+                    dev_v &press_v, dev_v &pmax_v, dev_v &vol_v, dev_v &spvol_v,
+                    dev_v &sie_v, dev_v &temp_v, dev_frac_v &frac_mass_v,
+                    ScratchV<int> &pte_idxs, ScratchV<int> &pte_mats,
+                    ScratchV<double> &press_pte, ScratchV<double> &vfrac_pte,
+                    ScratchV<double> &rho_pte, ScratchV<double> &sie_pte,
+                    ScratchV<double> &temp_pte, ScratchV<double> &solver_scratch,
+                    Kokkos::Experimental::UniqueToken<DES, KGlobal> &tokens,
+                    bool small_loop, final_functor &f_func);
+// rho e input
+void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v,
+                      Kokkos::View<EOS *, Llft> &eos_v, dev_v &press_v, dev_v &pmax_v,
+                      dev_v &sie_v, ScratchV<int> &pte_idxs, ScratchV<double> &press_pte,
+                      ScratchV<double> &vfrac_pte, ScratchV<double> &rho_pte,
+                      ScratchV<double> &sie_pte, ScratchV<double> &temp_pte,
+                      ScratchV<double> &solver_scratch,
+                      Kokkos::Experimental::UniqueToken<DES, KGlobal> &tokens,
+                      bool small_loop, init_functor &i_func, final_functor &f_func);
+} // namespace singularity
+#endif // PORTABILITY_STRATEGY_KOKKOS
+
+#endif // _SINGULARITY_EOS_EOS_GET_SG_EOS_HPP_
diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp
new file mode 100644
index 0000000000..f8680eacd3
--- /dev/null
+++ b/singularity-eos/eos/get_sg_eos_functors.hpp
@@ -0,0 +1,227 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#ifndef _SINGULARITY_EOS_EOS_GET_SG_EOS_LAMBDAS_HPP_
+#define _SINGULARITY_EOS_EOS_GET_SG_EOS_LAMBDAS_HPP_
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/eos/get_sg_eos.hpp>
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+
+namespace singularity {
+
+struct init_functor {
+ private:
+  dev_frac_v frac_mass_v;
+  ScratchV<int> pte_idxs;
+  indirection_v eos_offsets_v;
+  dev_frac_v frac_vol_v;
+  ScratchV<int> pte_mats;
+  ScratchV<double> vfrac_pte;
+  ScratchV<double> sie_pte;
+  ScratchV<double> temp_pte;
+  ScratchV<double> press_pte;
+  ScratchV<double> rho_pte;
+  dev_v spvol_v;
+  dev_v temp_v;
+  dev_v press_v;
+  dev_v sie_v;
+  int nmat;
+
+ public:
+  init_functor() = default;
+  init_functor(dev_frac_v &frac_mass_v_, ScratchV<int> &pte_idxs_,
+               indirection_v &eos_offsets_v_, dev_frac_v &frac_vol_v_,
+               ScratchV<int> &pte_mats_, ScratchV<double> &vfrac_pte_,
+               ScratchV<double> &sie_pte_, ScratchV<double> &temp_pte_,
+               ScratchV<double> &press_pte_, ScratchV<double> &rho_pte_, dev_v &spvol_v_,
+               dev_v &temp_v_, dev_v &press_v_, dev_v &sie_v_, int &nmat)
+      : frac_mass_v{frac_mass_v_}, pte_idxs{pte_idxs_}, eos_offsets_v{eos_offsets_v_},
+        frac_vol_v{frac_vol_v_}, pte_mats{pte_mats_}, vfrac_pte{vfrac_pte_},
+        sie_pte{sie_pte_}, temp_pte{temp_pte_}, press_pte{press_pte_}, rho_pte{rho_pte_},
+        spvol_v{spvol_v_}, temp_v{temp_v_}, press_v{press_v_}, sie_v{sie_v_}, nmat{nmat} {
+  }
+
+  PORTABLE_INLINE_FUNCTION
+  void operator()(const int i, const int tid, double &mass_sum, int &npte,
+                  const Real t_mult, const Real s_mult, const Real p_mult) const {
+    /* normalize mass fractions */
+    /* first find the mass sum */
+    /* also set idxs as the decrement of the eos offsets */
+    /* to take into account 1 based indexing in fortran */
+    for (int m = 0; m < nmat; ++m) {
+      mass_sum += frac_mass_v(i, m);
+      pte_idxs(tid, m) = eos_offsets_v(m) - 1;
+      frac_vol_v(i, m) = 0.0;
+    }
+    for (int m = 0; m < nmat; ++m) {
+      frac_mass_v(i, m) /= mass_sum;
+    }
+    /* set inputs */
+    npte = 0;
+    for (int m = 0; m < nmat; ++m) {
+      if (frac_mass_v(i, m) > 1.e-12) {
+        pte_idxs(tid, npte) = eos_offsets_v(m) - 1;
+        pte_mats(tid, npte) = m;
+        npte += 1;
+      }
+      vfrac_pte(tid, m) = 0.0;
+      sie_pte(tid, m) = 0.0;
+      temp_pte(tid, m) = 0.0;
+      press_pte(tid, m) = 0.0;
+    }
+    for (int mp = 0; mp < npte; ++mp) {
+      const int m = pte_mats(tid, mp);
+      rho_pte(tid, mp) = npte / spvol_v(i) * frac_mass_v(i, m);
+      vfrac_pte(tid, mp) = 1.0 / npte;
+      temp_pte(tid, mp) = temp_v(i) * ev2k * t_mult;
+      press_pte(tid, mp) = press_v(i) * p_mult;
+      sie_pte(tid, mp) = sie_v(i) * frac_mass_v(i, m) * s_mult;
+    }
+    return;
+  }
+};
+
+struct final_functor {
+ private:
+  dev_v temp_v;
+  dev_v press_v;
+  dev_v sie_v;
+  dev_v bmod_v;
+  dev_v cv_v;
+  dev_v dpde_v;
+  ScratchV<int> pte_mats;
+  ScratchV<double> press_pte;
+  ScratchV<double> vfrac_pte;
+  ScratchV<double> temp_pte;
+  ScratchV<double> sie_pte;
+  dev_frac_v frac_mass_v;
+  dev_frac_v frac_ie_v;
+  dev_frac_v frac_vol_v;
+  dev_v vol_v;
+  Kokkos::View<EOS *, Llft> eos_v;
+  ScratchV<int> pte_idxs;
+  ScratchV<double> rho_pte;
+  dev_frac_v frac_bmod_v;
+  dev_frac_v frac_cv_v;
+  dev_frac_v frac_dpde_v;
+  int nmat;
+  bool do_frac_bmod;
+  bool do_frac_cv;
+  bool do_frac_dpde;
+
+ public:
+  final_functor(dev_v &temp_v_, dev_v &press_v_, dev_v &sie_v_, dev_v &bmod_v_,
+                dev_v &cv_v_, dev_v &dpde_v_, ScratchV<int> &pte_mats_,
+                ScratchV<double> &press_pte_, ScratchV<double> &vfrac_pte_,
+                ScratchV<double> &temp_pte_, ScratchV<double> &sie_pte_,
+                dev_frac_v &frac_mass_v_, dev_frac_v &frac_ie_v_, dev_frac_v &frac_vol_v_,
+                dev_v &vol_v_, Kokkos::View<EOS *, Llft> &eos_v_,
+                ScratchV<int> &pte_idxs_, ScratchV<double> &rho_pte_,
+                dev_frac_v &frac_bmod_v_, dev_frac_v &frac_cv_v_,
+                dev_frac_v &frac_dpde_v_, int &nmat_, bool do_frac_bmod_,
+                bool do_frac_cv_, bool do_frac_dpde_)
+      : temp_v{temp_v_}, press_v{press_v_}, sie_v{sie_v_}, bmod_v{bmod_v_}, cv_v{cv_v_},
+        dpde_v{dpde_v_}, pte_mats{pte_mats_}, press_pte{press_pte_},
+        vfrac_pte{vfrac_pte_}, temp_pte{temp_pte_}, sie_pte{sie_pte_},
+        frac_mass_v{frac_mass_v_}, frac_ie_v{frac_ie_v_}, frac_vol_v{frac_vol_v_},
+        vol_v{vol_v_}, eos_v{eos_v_}, pte_idxs{pte_idxs_}, rho_pte{rho_pte_},
+        frac_bmod_v{frac_bmod_v_}, frac_cv_v{frac_cv_v_},
+        frac_dpde_v{frac_dpde_v_}, nmat{nmat_}, do_frac_bmod{do_frac_bmod_},
+        do_frac_cv{do_frac_cv_}, do_frac_dpde{do_frac_dpde_} {}
+
+ public:
+  PORTABLE_INLINE_FUNCTION
+  void operator()(const int i, const int tid, const int npte, const Real mass_sum,
+                  const Real t_mult, const Real s_mult, const Real p_mult,
+                  singularity::mix_impl::CacheAccessor &cache) const {
+    /* initialize averaged quantities to 0 */
+    const bool do_t = t_mult == 1.0;
+    const bool do_s = s_mult == 1.0;
+    const bool do_p = p_mult == 1.0;
+    if (do_t) {
+      temp_v(i) = 0.0;
+    }
+    if (do_p) {
+      press_v(i) = 0.0;
+    }
+    if (do_s) {
+      sie_v(i) = 0.0;
+    }
+    bmod_v(i) = 0.0;
+    cv_v(i) = 0.0;
+    dpde_v(i) = 0.0;
+    /* material loop for averaging and assigning per mat quantities */
+    for (int mp = 0; mp < npte; ++mp) {
+      const int m = pte_mats(tid, mp);
+      /* pressure contribution from material m */
+      press_v(i) += press_pte(tid, mp) * vfrac_pte(tid, mp) * p_mult;
+      /* temperature contribution from material m */
+      temp_v(i) += temp_pte(tid, mp) * vfrac_pte(tid, mp) * t_mult;
+      const Real ie_m = sie_pte(tid, mp) * frac_mass_v(i, m) * mass_sum;
+      /* sie contribution from material m */
+      sie_v(i) += ie_m * s_mult;
+      /* assign per material specific internal energy */
+      frac_ie_v(i, m) = ie_m;
+      /* assign volume fraction based on pte calculation */
+      frac_vol_v(i, m) = vfrac_pte(tid, mp) * vol_v(i);
+      /* calculate bulk modulus for material m */
+      const Real bmod_m = eos_v(pte_idxs(tid, mp))
+                              .BulkModulusFromDensityTemperature(
+                                  rho_pte(tid, mp), temp_pte(tid, mp), cache[mp]);
+      /* add bmod contribution from material m */
+      bmod_v(i) += bmod_m * vfrac_pte(tid, mp);
+      /* calculate specific heat for material m */
+      const Real cv_m = ev2k * eos_v(pte_idxs(tid, mp))
+                                   .SpecificHeatFromDensityTemperature(
+                                       rho_pte(tid, mp), temp_pte(tid, mp), cache[mp]);
+      /* add mass weighted contribution specific heat for material m */
+      cv_v(i) += cv_m * frac_mass_v(i, m);
+      /* calculate gruneisen parameter for material m */
+      const Real dpde_m = eos_v(pte_idxs(tid, mp))
+                              .GruneisenParamFromDensityTemperature(
+                                  rho_pte(tid, mp), temp_pte(tid, mp), cache[mp]);
+      /* add gruneisen param contribution from material m */
+      dpde_v(i) += dpde_m * vfrac_pte(tid, mp);
+      /* optionally assign per material quantities to per material arrays */
+      if (do_frac_bmod) {
+        frac_bmod_v(i, m) = bmod_m;
+      }
+      if (do_frac_cv) {
+        frac_cv_v(i, m) = cv_m;
+      }
+      if (do_frac_dpde) {
+        frac_dpde_v(i, m) = dpde_m;
+      }
+    }
+    if (do_t) {
+      temp_v(i) /= ev2k;
+    }
+    if (do_s) {
+      sie_v(i) /= mass_sum;
+    }
+    /* reset mass fractions to original values if not normalized to 1 */
+    for (int m = 0; m < nmat; ++m) {
+      frac_mass_v(i, m) *= mass_sum;
+    }
+    return;
+  }
+};
+
+} // namespace singularity
+
+#endif // PORTABILITY_STRATEGY_KOKKOS
+
+#endif // _SINGULARITY_EOS_EOS_GET_SG_EOS_LAMBDAS_HPP_
diff --git a/singularity-eos/eos/get_sg_eos_p_t.cpp b/singularity-eos/eos/get_sg_eos_p_t.cpp
new file mode 100644
index 0000000000..fb1db64d2d
--- /dev/null
+++ b/singularity-eos/eos/get_sg_eos_p_t.cpp
@@ -0,0 +1,92 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/closure/mixed_cell_models.hpp>
+#include <singularity-eos/eos/eos.hpp>
+#include <singularity-eos/eos/get_sg_eos.hpp>
+#include <singularity-eos/eos/get_sg_eos_functors.hpp>
+
+// any bmod vars, dpde, cv, frac vol, frac ie
+
+namespace singularity {
+void get_sg_eos_p_t(const char *name, int ncell, int nmat, indirection_v &offsets_v,
+                    indirection_v &eos_offsets_v, Kokkos::View<EOS *, Llft> &eos_v,
+                    dev_v &press_v, dev_v &pmax_v, dev_v &vol_v, dev_v &spvol_v,
+                    dev_v &sie_v, dev_v &temp_v, dev_frac_v &frac_mass_v,
+                    ScratchV<int> &pte_idxs, ScratchV<int> &pte_mats,
+                    ScratchV<double> &press_pte, ScratchV<double> &vfrac_pte,
+                    ScratchV<double> &rho_pte, ScratchV<double> &sie_pte,
+                    ScratchV<double> &temp_pte, ScratchV<double> &solver_scratch,
+                    Kokkos::Experimental::UniqueToken<DES, KGlobal> &tokens,
+                    bool small_loop, final_functor &f_func) {
+  portableFor(
+      name, 0, ncell, PORTABLE_LAMBDA(const int &iloop) {
+        // cell offset
+        const int i{offsets_v(iloop) - 1};
+        // get "thread-id" like thing with optimization
+        // for small loops
+        const int32_t token{tokens.acquire()};
+        const int32_t tid{small_loop ? iloop : token};
+        // caching mechanism
+        singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0));
+        double mass_sum{0.0};
+        // normalize mass fractions
+        // first find the mass sum
+        // also set idxs as the decrement of the eos offsets
+        // to take into account 1 based indexing in fortran
+        for (int m = 0; m < nmat; ++m) {
+          mass_sum += frac_mass_v(i, m);
+          pte_idxs(tid, m) = eos_offsets_v(m) - 1;
+          pte_mats(tid, m) = m;
+          temp_pte(tid, m) = temp_v(i) * ev2k;
+          press_pte(tid, m) = press_v(i);
+        }
+        for (int m = 0; m < nmat; ++m) {
+          frac_mass_v(i, m) /= mass_sum;
+        }
+        // do r-e of pt for each mat
+        singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
+        Real vfrac_tot{0.0};
+        Real sie_tot{0.0};
+        for (int m = 0; m < nmat; ++m) {
+          // obtain rho and sie from P-T
+          eos_inx[m].DensityEnergyFromPressureTemperature(
+              press_pte(tid, m), temp_pte(tid, m), cache[m], rho_pte(tid, m),
+              sie_pte(tid, m));
+          // assign volume fractions
+          // this is a physical volume
+          vfrac_pte(tid, m) = frac_mass_v(i, m) / rho_pte(tid, m) * mass_sum;
+          vfrac_tot += vfrac_pte(tid, m);
+          // add internal energy component
+          sie_tot += sie_pte(tid, m) * frac_mass_v(i, m);
+        }
+        // assign volume, etc.
+        // total sie is known
+        sie_v(i) = sie_tot;
+        vol_v(i) = vfrac_tot;
+        spvol_v(i) = vol_v(i) / mass_sum;
+        for (int m = 0; m < nmat; ++m) {
+          vfrac_pte(tid, m) /= vfrac_tot;
+        }
+        // assign remaining outputs
+        f_func(i, tid, nmat, mass_sum, 0.0, 0.0, 0.0, cache);
+        // assign max pressure
+        pmax_v(i) = press_v(i) > pmax_v(i) ? press_v(i) : pmax_v(i);
+        // release the token used for scratch arrays
+        tokens.release(token);
+      });
+  return;
+}
+} // namespace singularity
diff --git a/singularity-eos/eos/get_sg_eos_rho_e.cpp b/singularity-eos/eos/get_sg_eos_rho_e.cpp
new file mode 100644
index 0000000000..43049e31a7
--- /dev/null
+++ b/singularity-eos/eos/get_sg_eos_rho_e.cpp
@@ -0,0 +1,74 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/closure/mixed_cell_models.hpp>
+#include <singularity-eos/eos/eos.hpp>
+#include <singularity-eos/eos/get_sg_eos.hpp>
+#include <singularity-eos/eos/get_sg_eos_functors.hpp>
+
+namespace singularity {
+void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v,
+                      Kokkos::View<EOS *, Llft> &eos_v, dev_v &press_v, dev_v &pmax_v,
+                      dev_v &sie_v, ScratchV<int> &pte_idxs, ScratchV<double> &press_pte,
+                      ScratchV<double> &vfrac_pte, ScratchV<double> &rho_pte,
+                      ScratchV<double> &sie_pte, ScratchV<double> &temp_pte,
+                      ScratchV<double> &solver_scratch,
+                      Kokkos::Experimental::UniqueToken<DES, KGlobal> &tokens,
+                      bool small_loop, init_functor &i_func, final_functor &f_func) {
+  portableFor(
+      name, 0, ncell, PORTABLE_LAMBDA(const int &iloop) {
+        // cell offset
+        const int i{offsets_v(iloop) - 1};
+        // get "thread-id" like thing with optimization
+        // for small loops
+        const int32_t token{tokens.acquire()};
+        const int32_t tid{small_loop ? iloop : token};
+        double mass_sum{0.0};
+        int npte{0};
+        // initialize values for solver / lookup
+        i_func(i, tid, mass_sum, npte, 0.0, 1.0, 0.0);
+        // get cache from offsets into scratch
+        const int neq = npte + 1;
+        singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) +
+                                                   neq * (neq + 4) + 2 * npte);
+        if (npte > 1) {
+          // create solver lambda
+          // eos accessor
+          singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
+          // reset inputs
+          PTESolverRhoT<singularity::EOSAccessor_, Real *, Real **> method(
+              npte, eos_inx, 1.0, sie_v(i), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
+              &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache,
+              &solver_scratch(tid, 0));
+          const bool res_{PTESolver(method)};
+        } else {
+          // pure cell (nmat = 1)
+          temp_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
+                                 .TemperatureFromDensityInternalEnergy(
+                                     rho_pte(tid, 0), sie_pte(tid, 0), cache[0]);
+          press_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
+                                  .PressureFromDensityTemperature(
+                                      rho_pte(tid, 0), temp_pte(tid, 0), cache[0]);
+        }
+        // assign outputs
+        f_func(i, tid, npte, mass_sum, 1.0, 0.0, 1.0, cache);
+        // assign max pressure
+        pmax_v(i) = press_v(i) > pmax_v(i) ? press_v(i) : pmax_v(i);
+        // release the token used for scratch arrays
+        tokens.release(token);
+      });
+  return;
+}
+} // namespace singularity
diff --git a/singularity-eos/eos/get_sg_eos_rho_p.cpp b/singularity-eos/eos/get_sg_eos_rho_p.cpp
new file mode 100644
index 0000000000..dc6da80ecc
--- /dev/null
+++ b/singularity-eos/eos/get_sg_eos_rho_p.cpp
@@ -0,0 +1,95 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/closure/mixed_cell_models.hpp>
+#include <singularity-eos/eos/eos.hpp>
+#include <singularity-eos/eos/get_sg_eos.hpp>
+#include <singularity-eos/eos/get_sg_eos_functors.hpp>
+
+namespace singularity {
+void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v,
+                      Kokkos::View<EOS *, Llft> &eos_v, dev_v &press_v, dev_v &pmax_v,
+                      dev_v &sie_v, dev_frac_v &frac_mass_v, ScratchV<int> &pte_idxs,
+                      ScratchV<int> &pte_mats, ScratchV<double> &press_pte,
+                      ScratchV<double> &vfrac_pte, ScratchV<double> &rho_pte,
+                      ScratchV<double> &sie_pte, ScratchV<double> &temp_pte,
+                      ScratchV<double> &solver_scratch,
+                      Kokkos::Experimental::UniqueToken<DES, KGlobal> &tokens,
+                      bool small_loop, init_functor &i_func, final_functor &f_func) {
+  portableFor(
+      name, 0, ncell, PORTABLE_LAMBDA(const int &iloop) {
+        // cell offset
+        const int i{offsets_v(iloop) - 1};
+        // get "thread-id" like thing with optimization
+        // for small loops
+        const int32_t token{tokens.acquire()};
+        const int32_t tid{small_loop ? iloop : token};
+        double mass_sum{0.0};
+        int npte{0};
+        // initialize values for solver / lookup
+        i_func(i, tid, mass_sum, npte, 0.0, 0.0, 1.0);
+        Real sie_tot_true{0.0};
+        const int neq = npte + 1;
+        singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) +
+                                                   neq * (neq + 4) + 2 * npte);
+        if (npte > 1) {
+          // create solver lambda
+          // eos accessor
+          singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
+          PTESolverFixedP<singularity::EOSAccessor_, Real *, Real *> method(
+              npte, eos_inx, 1.0, press_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
+              &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache[0],
+              &solver_scratch(tid, 0));
+          const bool res_{PTESolver(method)};
+          // calculate total sie
+          for (int mp = 0; mp < npte; ++mp) {
+            const int m = pte_mats(tid, mp);
+            sie_tot_true += sie_pte(tid, mp) * frac_mass_v(i, m);
+          }
+        } else {
+          // pure cell (nmat = 1)
+          // calculate sie from single eos
+          auto p_from_t = [&](const Real &t_i) {
+            return eos_v(pte_idxs(tid, 0))
+                .PressureFromDensityTemperature(rho_pte(tid, 0), t_i, cache[0]);
+          };
+          // calculate sie root bounds
+          Real r_rho{}, r_temp{}, r_sie{}, r_press{}, r_cv{}, r_bmod{};
+          Real r_dpde{}, r_dvdt{};
+          eos_v(pte_idxs(tid, 0))
+              .ValuesAtReferenceState(r_rho, r_temp, r_sie, r_press, r_cv, r_bmod, r_dpde,
+                                      r_dvdt);
+          Real temp_i;
+          RootFinding1D::findRoot(p_from_t, press_v(i), r_temp, 1.e-10 * r_temp,
+                                  1.e10 * r_temp, 1.e-12, 1.e-12, temp_i);
+          sie_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
+                                .InternalEnergyFromDensityTemperature(rho_pte(tid, 0),
+                                                                      temp_i, cache[0]);
+          sie_tot_true = sie_pte(tid, 0);
+          // set temperature
+          temp_pte(tid, 0) = temp_i;
+        }
+        // sie calculate explicitly already
+        sie_v(i) = sie_tot_true;
+        // assign remaining outputs
+        f_func(i, tid, npte, mass_sum, 1.0, 0.0, 0.0, cache);
+        // assign max pressure
+        pmax_v(i) = press_v(i) > pmax_v(i) ? press_v(i) : pmax_v(i);
+        // release the token used for scratch arrays
+        tokens.release(token);
+      });
+  return;
+}
+} // namespace singularity
diff --git a/singularity-eos/eos/get_sg_eos_rho_t.cpp b/singularity-eos/eos/get_sg_eos_rho_t.cpp
new file mode 100644
index 0000000000..f1b91eed1f
--- /dev/null
+++ b/singularity-eos/eos/get_sg_eos_rho_t.cpp
@@ -0,0 +1,86 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/closure/mixed_cell_models.hpp>
+#include <singularity-eos/eos/eos.hpp>
+#include <singularity-eos/eos/get_sg_eos.hpp>
+#include <singularity-eos/eos/get_sg_eos_functors.hpp>
+
+namespace singularity {
+void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v,
+                      Kokkos::View<EOS *, Llft> &eos_v, dev_v &press_v, dev_v &pmax_v,
+                      dev_v &sie_v, dev_frac_v &frac_mass_v, ScratchV<int> &pte_idxs,
+                      ScratchV<int> &pte_mats, ScratchV<double> &press_pte,
+                      ScratchV<double> &vfrac_pte, ScratchV<double> &rho_pte,
+                      ScratchV<double> &sie_pte, ScratchV<double> &temp_pte,
+                      ScratchV<double> &solver_scratch,
+                      Kokkos::Experimental::UniqueToken<DES, KGlobal> &tokens,
+                      bool small_loop, init_functor &i_func, final_functor &f_func) {
+  portableFor(
+      name, 0, ncell, PORTABLE_LAMBDA(const int &iloop) {
+        // cell offset
+        const int i{offsets_v(iloop) - 1};
+        // get "thread-id" like thing with optimization
+        // for small loops
+        const int32_t token{tokens.acquire()};
+        const int32_t tid{small_loop ? iloop : token};
+        double mass_sum{0.0};
+        int npte{0};
+        // initialize values for solver / lookup
+        i_func(i, tid, mass_sum, npte, 1.0, 0.0, 0.0);
+        // calculate pte condition (lookup for 1 mat cell)
+        Real sie_tot_true{0.0};
+        const int neq = npte;
+        singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) +
+                                                   neq * (neq + 4) + 2 * npte);
+        if (npte > 1) {
+          // create solver lambda
+          // eos accessor
+          singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
+          PTESolverFixedT<singularity::EOSAccessor_, Real *, Real **> method(
+              npte, eos_inx, 1.0, temp_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
+              &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache,
+              &solver_scratch(tid, 0));
+          const bool res_{PTESolver(method)};
+          // calculate total internal energy
+          for (int mp = 0; mp < npte; ++mp) {
+            const int m = pte_mats(tid, mp);
+            sie_tot_true += sie_pte(tid, mp) * frac_mass_v(i, m);
+          }
+        } else {
+          // pure cell (nmat = 1)
+          // calculate sie from single eos
+          sie_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
+                                .InternalEnergyFromDensityTemperature(
+                                    rho_pte(tid, 0), temp_pte(tid, 0), cache[0]);
+          // set total sie to material 0 value
+          sie_tot_true = sie_pte(tid, 0);
+          // set pressure
+          press_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
+                                  .PressureFromDensityTemperature(
+                                      rho_pte(tid, 0), temp_pte(tid, 0), cache[0]);
+        }
+        // total sie is known
+        sie_v(i) = sie_tot_true;
+        // assign quantities
+        f_func(i, tid, npte, mass_sum, 0.0, 0.0, 1.0, cache);
+        // assign max pressure
+        pmax_v(i) = press_v(i) > pmax_v(i) ? press_v(i) : pmax_v(i);
+        // release the token used for scratch arrays
+        tokens.release(token);
+      });
+  return;
+}
+} // namespace singularity
diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
index 5e7f6b033a..a01fdf599b 100644
--- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp
+++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
@@ -166,6 +166,7 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
         t_.PressureFromDensityTemperature(rho * rho_unit_, temp * temp_unit_, lambda);
     return inv_press_unit_ * P;
   }
+  PORTABLE_FUNCTION
   Real EntropyFromDensityTemperature(const Real rho, const Real temp,
                                      Real *lambda = nullptr) const {
     const Real S =
@@ -269,12 +270,13 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
     printf("Units = %e %e %e %e\n", rho_unit_, sie_unit_, temp_unit_, press_unit_);
   }
 
-  PORTABLE_FORCEINLINE_FUNCTION
-  bool IsModified() const { return true; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  T UnmodifyOnce() { return t_; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  auto GetUnmodifiedObject() { return t_.GetUnmodifiedObject(); }
+  inline constexpr bool IsModified() const { return true; }
+
+  inline constexpr T UnmodifyOnce() { return t_; }
+
+  inline constexpr decltype(auto) GetUnmodifiedObject() {
+    return t_.GetUnmodifiedObject();
+  }
 
  private:
   T t_;
diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp
index b6cb1d50f8..a0111b8ea5 100644
--- a/singularity-eos/eos/modifiers/ramps_eos.hpp
+++ b/singularity-eos/eos/modifiers/ramps_eos.hpp
@@ -246,10 +246,10 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
     static auto const name = singularity::mfuncname::member_func_name(
         typeid(BilinearRampEOS<T>).name(), __func__);
     static auto const cname = name.c_str();
-
+    auto const copy = *this;
     portableFor(
         cname, 0, num, PORTABLE_LAMBDA(const int i) {
-          const Real p_ramp = get_ramp_pressure(rhos[i]);
+          const Real p_ramp = copy.get_ramp_pressure(rhos[i]);
           pressures[i] = std::max(pressures[i], p_ramp);
         });
   }
@@ -265,10 +265,10 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
     static auto const name = singularity::mfuncname::member_func_name(
         typeid(BilinearRampEOS<T>).name(), __func__);
     static auto const cname = name.c_str();
-
+    auto const copy = *this;
     portableFor(
         cname, 0, num, PORTABLE_LAMBDA(const int i) {
-          const Real p_ramp = get_ramp_pressure(rhos[i]);
+          const Real p_ramp = copy.get_ramp_pressure(rhos[i]);
           pressures[i] = std::max(pressures[i], p_ramp);
         });
   }
@@ -305,12 +305,12 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
     static auto const name = singularity::mfuncname::member_func_name(
         typeid(BilinearRampEOS<T>).name(), __func__);
     static auto const cname = name.c_str();
-
+    auto const copy = *this;
     portableFor(
         cname, 0, num, PORTABLE_LAMBDA(const int i) {
-          const Real p_ramp = get_ramp_pressure(rhos[i]);
+          const Real p_ramp = copy.get_ramp_pressure(rhos[i]);
           if (pressures[i] < p_ramp) {
-            bmods[i] = rhos[i] * get_ramp_dpdrho(rhos[i]);
+            bmods[i] = rhos[i] * copy.get_ramp_dpdrho(rhos[i]);
           }
         });
   }
@@ -329,12 +329,12 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
     static auto const name = singularity::mfuncname::member_func_name(
         typeid(BilinearRampEOS<T>).name(), __func__);
     static auto const cname = name.c_str();
-
+    auto const copy = *this;
     portableFor(
         cname, 0, num, PORTABLE_LAMBDA(const int i) {
-          const Real p_ramp = get_ramp_pressure(rhos[i]);
+          const Real p_ramp = copy.get_ramp_pressure(rhos[i]);
           if (pressures[i] < p_ramp) {
-            bmods[i] = rhos[i] * get_ramp_dpdrho(rhos[i]);
+            bmods[i] = rhos[i] * copy.get_ramp_dpdrho(rhos[i]);
           }
         });
   }
@@ -428,12 +428,13 @@ class BilinearRampEOS : public EosBase<BilinearRampEOS<T>> {
   // Vector functions that overload the scalar versions declared here.
   SG_ADD_BASE_CLASS_USINGS(BilinearRampEOS<T>)
 
-  PORTABLE_FORCEINLINE_FUNCTION
-  bool IsModified() const { return true; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  T UnmodifyOnce() { return t_; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  auto GetUnmodifiedObject() { return t_.GetUnmodifiedObject(); }
+  inline constexpr bool IsModified() const { return true; }
+
+  inline constexpr T UnmodifyOnce() { return t_; }
+
+  inline constexpr decltype(auto) GetUnmodifiedObject() {
+    return t_.GetUnmodifiedObject();
+  }
 
  private:
   T t_;
diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp
index e1fe0b5fa7..174349f5fd 100644
--- a/singularity-eos/eos/modifiers/relativistic_eos.hpp
+++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp
@@ -185,12 +185,13 @@ class RelativisticEOS : public EosBase<RelativisticEOS<T>> {
     t_.ValuesAtReferenceState(rho, temp, sie, press, cv, bmod, dpde, dvdt, lambda);
   }
 
-  PORTABLE_FORCEINLINE_FUNCTION
-  bool IsModified() const { return true; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  T UnmodifyOnce() { return t_; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  auto GetUnmodifiedObject() { return t_.GetUnmodifiedObject(); }
+  inline constexpr bool IsModified() const { return true; }
+
+  inline constexpr T UnmodifyOnce() { return t_; }
+
+  inline constexpr decltype(auto) GetUnmodifiedObject() {
+    return t_.GetUnmodifiedObject();
+  }
 
  private:
   T t_;
diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp
index 730aa1e7d0..c1397c6e85 100644
--- a/singularity-eos/eos/modifiers/scaled_eos.hpp
+++ b/singularity-eos/eos/modifiers/scaled_eos.hpp
@@ -92,6 +92,7 @@ class ScaledEOS : public EosBase<ScaledEOS<T>> {
                                          Real *lambda = nullptr) const {
     return t_.PressureFromDensityInternalEnergy(scale_ * rho, inv_scale_ * sie, lambda);
   }
+  PORTABLE_FUNCTION
   Real EntropyFromDensityInternalEnergy(const Real rho, const Real sie,
                                         Real *lambda = nullptr) const {
     return scale_ *
@@ -334,12 +335,13 @@ class ScaledEOS : public EosBase<ScaledEOS<T>> {
     return t_.MinimumTemperature();
   }
 
-  PORTABLE_FORCEINLINE_FUNCTION
-  bool IsModified() const { return true; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  T UnmodifyOnce() { return t_; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  auto GetUnmodifiedObject() { return t_.GetUnmodifiedObject(); }
+  inline constexpr bool IsModified() const { return true; }
+
+  inline constexpr T UnmodifyOnce() { return t_; }
+
+  inline constexpr decltype(auto) GetUnmodifiedObject() {
+    return t_.GetUnmodifiedObject();
+  }
 
  private:
   T t_;
diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp
index be2d0b4773..45c3ed4d4b 100644
--- a/singularity-eos/eos/modifiers/shifted_eos.hpp
+++ b/singularity-eos/eos/modifiers/shifted_eos.hpp
@@ -158,16 +158,19 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
     static auto const name =
         singularity::mfuncname::member_func_name(typeid(ShiftedEOS<T>).name(), __func__);
     static auto const cname = name.c_str();
+    const auto shift_val = shift_;
     portableFor(
-        cname, 0, num, PORTABLE_LAMBDA(const int i) { shifted[i] = sies[i] - shift_; });
+        cname, 0, num,
+        PORTABLE_LAMBDA(const int i) { shifted[i] = sies[i] - shift_val; });
   }
 
   inline void unshift_sies(Real *sies, const int num) const {
     static auto const name =
         singularity::mfuncname::member_func_name(typeid(ShiftedEOS<T>).name(), __func__);
     static auto const cname = name.c_str();
+    const auto shift_val = shift_;
     portableFor(
-        cname, 0, num, PORTABLE_LAMBDA(const int i) { sies[i] += shift_; });
+        cname, 0, num, PORTABLE_LAMBDA(const int i) { sies[i] += shift_val; });
   }
 
   template <typename LambdaIndexer>
@@ -345,12 +348,13 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
     return t_.MinimumTemperature();
   }
 
-  PORTABLE_FORCEINLINE_FUNCTION
-  bool IsModified() const { return true; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  T UnmodifyOnce() { return t_; }
-  PORTABLE_FORCEINLINE_FUNCTION
-  auto GetUnmodifiedObject() { return t_.GetUnmodifiedObject(); }
+  inline constexpr bool IsModified() const { return true; }
+
+  inline constexpr T UnmodifyOnce() { return t_; }
+
+  inline constexpr decltype(auto) GetUnmodifiedObject() {
+    return t_.GetUnmodifiedObject();
+  }
 
  private:
   T t_;
diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp
index 54867736f0..128e0e0e95 100644
--- a/singularity-eos/eos/singularity_eos.cpp
+++ b/singularity-eos/eos/singularity_eos.cpp
@@ -12,20 +12,14 @@
 // publicly and display publicly, and to permit others to do so.
 //------------------------------------------------------------------------------
 
-#include <algorithm>
 #include <cassert>
-#include <map>
+#include <ports-of-call/portability.hpp>
+#include <singularity-eos/eos/eos.hpp>
 #include <singularity-eos/eos/eos_builder.hpp>
 #include <singularity-eos/eos/singularity_eos.hpp>
 
 using namespace singularity;
 
-template <typename T,
-          typename = typename std::enable_if<std::is_floating_point<T>::value>::type>
-constexpr bool isfinite(const T &a) {
-  return (a == a) && ((a == 0) || (a != 2 * a));
-}
-
 int init_sg_eos(const int nmat, EOS *&eos) {
 #ifdef PORTABILITY_STRATEGY_KOKKOS
   if (!Kokkos::is_initialized()) Kokkos::initialize();
@@ -272,584 +266,6 @@ int init_sg_eospac(const int matindex, EOS *eos, const int id) {
 
 #undef SGAPPLYMOD
 
-#ifdef PORTABILITY_STRATEGY_KOKKOS
-using Lrgt = Kokkos::LayoutRight;
-using Llft = Kokkos::LayoutLeft;
-template <typename T>
-using ScratchV = Kokkos::View<T **, Lrgt>;
-#endif // PORTABILITY_STRATEGY_KOKKOS
-
-// mapping from EAP integer to
-static const std::map<const int, const unsigned long> EAPInputToBD = {
-    {-3, thermalqs::temperature | thermalqs::density},
-    {-2, thermalqs::density | thermalqs::pressure},
-    {-1, thermalqs::pressure | thermalqs::temperature},
-    {0, thermalqs::specific_internal_energy | thermalqs::density},
-    {1, thermalqs::specific_internal_energy | thermalqs::density},
-};
-
-#ifdef PORTABILITY_STRATEGY_KOKKOS
-namespace singularity {
-struct EOSAccessor_ {
-  PORTABLE_INLINE_FUNCTION
-  EOSAccessor_(const Kokkos::View<EOS *, Llft> &eos_v, int *mats)
-      : eos_v_(eos_v), mats_(mats) {}
-  PORTABLE_INLINE_FUNCTION
-  EOS &operator[](const int m) const { return eos_v_(mats_[m]); }
-  Kokkos::View<EOS *, Llft> eos_v_;
-  int *mats_;
-};
-} // namespace singularity
-#endif // PORTABILITY_STRATEGY_KOKKOS
-// EAP centric arguments and function signature
-
-int get_sg_eos( // sizing information
-    int nmat, int ncell, int cell_dim,
-    // Input parameters
-    int input_int,
-    // eos index offsets
-    int *eos_offsets,
-    // equation of state array
-    EOS *eos,
-    // index offsets
-    int *offsets,
-    // per cell quantities
-    double *press, double *pmax, double *vol, double *spvol, double *sie, double *temp,
-    double *bmod, double *dpde, double *cv,
-    // per material quantities
-    double *frac_mass, double *frac_vol, double *frac_ie,
-    // optional per material quantities
-    double *frac_bmod, double *frac_dpde, double *frac_cv) {
-  // printBacktrace();
-  // kernel return value will be the number of failures
-  int ret{0};
-  // handle optionals
-  bool do_frac_bmod{true};
-  if (frac_bmod == NULL || frac_bmod == nullptr) do_frac_bmod = false;
-  bool do_frac_dpde{true};
-  if (frac_dpde == NULL || frac_dpde == nullptr) do_frac_dpde = false;
-  bool do_frac_cv{true};
-  if (frac_dpde == NULL || frac_dpde == nullptr) do_frac_cv = false;
-  // get inputs
-  enum class input_condition {
-    RHO_T_INPUT = -3,
-    RHO_P_INPUT = -2,
-    P_T_INPUT = -1,
-    NORM_RHO_E_INPUT = 0,
-    RHO_E_INPUT = 1
-  };
-  const auto input{EAPInputToBD.at(input_int)};
-  const auto output = thermalqs::all_values - input;
-  const bool p_is_inp{static_cast<bool>(input & thermalqs::pressure)};
-  const bool r_is_inp{static_cast<bool>(input & thermalqs::density)};
-  const bool t_is_inp{static_cast<bool>(input & thermalqs::temperature)};
-  const bool s_is_inp{static_cast<bool>(input & thermalqs::specific_internal_energy)};
-  constexpr const Real min_guess_density{1.e-20};
-#ifdef PORTABILITY_STRATEGY_KOKKOS
-  // convenience aliases
-  using Unmgd = Kokkos::MemoryUnmanaged;
-  using HS = Kokkos::HostSpace;
-  using Kokkos::create_mirror_view_and_copy;
-  using host_v = Kokkos::View<double *, Llft, HS, Unmgd>;
-  using dev_v = Kokkos::View<double *, Llft>;
-  using host_frac_v = Kokkos::View<double **, Llft, HS, Unmgd>;
-  using dev_frac_v = Kokkos::View<double **, Llft>;
-  using DES = Kokkos::DefaultExecutionSpace;
-  using DMS = DES::memory_space;
-  using Kokkos::MemoryTraits;
-  constexpr const unsigned int ra{0 | Kokkos::RandomAccess};
-  // constexpr const unsigned int ra_u{Kokkos::Unmanaged | Kokkos::RandomAccess};
-  using VAWI = Kokkos::ViewAllocateWithoutInitializing;
-  using Kokkos::deep_copy;
-  static constexpr const double ev2k = 1.160451930280894026e4;
-  // convert pointers to host side views
-  Kokkos::View<int *, Llft, HS, Unmgd> eos_offsets_hv(eos_offsets, nmat);
-  Kokkos::View<int *, Llft, HS, Unmgd> offsets_hv(offsets, ncell);
-  host_v press_hv(press, cell_dim);
-  host_v pmax_hv(pmax, cell_dim);
-  host_v vol_hv(vol, cell_dim);
-  host_v spvol_hv(spvol, cell_dim);
-  host_v sie_hv(sie, cell_dim);
-  host_v temp_hv(temp, cell_dim);
-  host_v bmod_hv(bmod, cell_dim);
-  host_v dpde_hv(dpde, cell_dim);
-  host_v cv_hv(cv, cell_dim);
-  host_frac_v frac_mass_hv(frac_mass, cell_dim, nmat);
-  host_frac_v frac_vol_hv(frac_vol, cell_dim, nmat);
-  host_frac_v frac_ie_hv(frac_ie, cell_dim, nmat);
-  host_frac_v frac_bmod_hv, frac_dpde_hv, frac_cv_hv;
-  if (do_frac_bmod) frac_bmod_hv = host_frac_v(frac_bmod, cell_dim, nmat);
-  if (do_frac_dpde) frac_dpde_hv = host_frac_v(frac_dpde, cell_dim, nmat);
-  if (do_frac_cv) frac_cv_hv = host_frac_v(frac_cv, cell_dim, nmat);
-
-  // get device views if necessary
-  Kokkos::View<const int *, Llft, MemoryTraits<ra>> offsets_v{
-      create_mirror_view_and_copy(DMS(), offsets_hv)};
-  Kokkos::View<const int *, Llft, MemoryTraits<ra>> eos_offsets_v{
-      create_mirror_view_and_copy(DMS(), eos_offsets_hv)};
-  dev_v press_v{create_mirror_view_and_copy(DMS(), press_hv)};
-  dev_v pmax_v{create_mirror_view_and_copy(DMS(), pmax_hv)};
-  dev_v spvol_v{create_mirror_view_and_copy(DMS(), spvol_hv)};
-  dev_v vol_v{create_mirror_view_and_copy(DMS(), vol_hv)};
-  dev_v sie_v{create_mirror_view_and_copy(DMS(), sie_hv)};
-  dev_v temp_v{create_mirror_view_and_copy(DMS(), temp_hv)};
-  dev_v bmod_v{create_mirror_view_and_copy(DMS(), bmod_hv)};
-  dev_v dpde_v{create_mirror_view_and_copy(DMS(), dpde_hv)};
-  dev_v cv_v{create_mirror_view_and_copy(DMS(), cv_hv)};
-  dev_frac_v frac_mass_v{create_mirror_view_and_copy(DMS(), frac_mass_hv)};
-  dev_frac_v frac_vol_v{create_mirror_view_and_copy(DMS(), frac_vol_hv)};
-  dev_frac_v frac_ie_v{create_mirror_view_and_copy(DMS(), frac_ie_hv)};
-  dev_frac_v frac_bmod_v, frac_dpde_v, frac_cv_v;
-  if (do_frac_bmod) frac_bmod_v = create_mirror_view_and_copy(DMS(), frac_bmod_hv);
-  if (do_frac_dpde) frac_dpde_v = create_mirror_view_and_copy(DMS(), frac_dpde_hv);
-  if (do_frac_cv) frac_cv_v = create_mirror_view_and_copy(DMS(), frac_cv_hv);
-  // array of eos's
-  const auto eos_nmat{*std::max_element(eos_offsets, eos_offsets + nmat)};
-  Kokkos::View<EOS *, Llft, HS, Unmgd> eos_hv(eos, eos_nmat);
-  Kokkos::View<EOS *, Llft> eos_v{create_mirror_view_and_copy(DMS(), eos_hv)};
-  // TODO: make this a scatter view
-  constexpr auto at_int_full{0 | Kokkos::Atomic};
-#ifdef KOKKOS_ENABLE_SERIAL
-  constexpr auto at_int{std::is_same<DES, Kokkos::Serial>::value ? 0 : at_int_full};
-#else  // KOKKOS_ENABLE_SERIAL
-  // assume atomics are required for correctness if serial not even enabled
-  constexpr auto at_int{at_int_full};
-#endif // KOKKOS_ENABLE_SERIAL
-  // set up scratch arrays
-  constexpr auto KGlobal = Kokkos::Experimental::UniqueTokenScope::Global;
-  Kokkos::Experimental::UniqueToken<DES, KGlobal> tokens{};
-  using VAWI = Kokkos::ViewAllocateWithoutInitializing;
-
-  const bool small_loop{tokens.size() > ncell};
-  const decltype(tokens)::size_type scratch_size{std::min(tokens.size(), ncell)};
-  ScratchV<int> pte_mats(VAWI("PTE::scratch mats"), scratch_size, nmat);
-  ScratchV<int> pte_idxs(VAWI("PTE::scratch idxs"), scratch_size, nmat);
-  ScratchV<double> mass_pte(VAWI("PTE::scratch mass"), scratch_size, nmat);
-  ScratchV<double> sie_pte(VAWI("PTE::scratch sie"), scratch_size, nmat);
-  ScratchV<double> vfrac_pte(VAWI("PTE::scratch vfrac"), scratch_size, nmat);
-  ScratchV<double> temp_pte(VAWI("PTE::scratch temp"), scratch_size, nmat);
-  ScratchV<double> press_pte(VAWI("PTE::scratch press"), scratch_size, nmat);
-  ScratchV<double> rho_pte(VAWI("PTE::scratch rho"), scratch_size, nmat);
-  // declare solver scratch, allocate in the case statement
-  int pte_solver_scratch_size{};
-  ScratchV<double> solver_scratch;
-
-  // create helper lambdas to reduce code duplication
-  const auto init_lambda =
-      PORTABLE_LAMBDA(const int i, const int tid, double &mass_sum, int &npte,
-                      const Real t_mult, const Real s_mult, const Real p_mult) {
-    // normalize mass fractions
-    // first find the mass sum
-    // also set idxs as the decrement of the eos offsets
-    // to take into account 1 based indexing in fortran
-    for (int m = 0; m < nmat; ++m) {
-      mass_sum += frac_mass_v(i, m);
-      pte_idxs(tid, m) = eos_offsets_v(m) - 1;
-      frac_vol_v(i, m) = 0.0;
-    }
-    for (int m = 0; m < nmat; ++m) {
-      frac_mass_v(i, m) /= mass_sum;
-    }
-    // set inputs
-    npte = 0;
-    for (int m = 0; m < nmat; ++m) {
-      if (frac_mass_v(i, m) > 1.e-12) {
-        pte_idxs(tid, npte) = eos_offsets_v(m) - 1;
-        pte_mats(tid, npte) = m;
-        npte += 1;
-      }
-      vfrac_pte(tid, m) = 0.0;
-      sie_pte(tid, m) = 0.0;
-      temp_pte(tid, m) = 0.0;
-      press_pte(tid, m) = 0.0;
-    }
-    for (int mp = 0; mp < npte; ++mp) {
-      const int m = pte_mats(tid, mp);
-      rho_pte(tid, mp) = npte / spvol_v(i) * frac_mass_v(i, m);
-      vfrac_pte(tid, mp) = 1.0 / npte;
-      temp_pte(tid, mp) = temp_v(i) * ev2k * t_mult;
-      press_pte(tid, mp) = press_v(i) * p_mult;
-      sie_pte(tid, mp) = sie_v(i) * frac_mass_v(i, m) * s_mult;
-    }
-  };
-  const auto final_lambda = PORTABLE_LAMBDA(
-      const int i, const int tid, const int npte, const Real mass_sum, const Real t_mult,
-      const Real s_mult, const Real p_mult, singularity::mix_impl::CacheAccessor &cache) {
-    // initialize averaged quantities to 0
-    const bool do_t = t_mult == 1.0;
-    const bool do_s = s_mult == 1.0;
-    const bool do_p = p_mult == 1.0;
-    if (do_t) {
-      temp_v(i) = 0.0;
-    }
-    if (do_p) {
-      press_v(i) = 0.0;
-    }
-    if (do_s) {
-      sie_v(i) = 0.0;
-    }
-    bmod_v(i) = 0.0;
-    cv_v(i) = 0.0;
-    dpde_v(i) = 0.0;
-    // material loop for averaging and assigning per mat quantities
-    for (int mp = 0; mp < npte; ++mp) {
-      const int m = pte_mats(tid, mp);
-      // pressure contribution from material m
-      press_v(i) += press_pte(tid, mp) * vfrac_pte(tid, mp) * p_mult;
-      // temperature contribution from material m
-      temp_v(i) += temp_pte(tid, mp) * vfrac_pte(tid, mp) * t_mult;
-      const Real ie_m = sie_pte(tid, mp) * frac_mass_v(i, m) * mass_sum;
-      // sie contribution from material m
-      sie_v(i) += ie_m * s_mult;
-      // assign per material specific internal energy
-      frac_ie_v(i, m) = ie_m;
-      // assign volume fraction based on pte calculation
-      frac_vol_v(i, m) = vfrac_pte(tid, mp) * vol_v(i);
-      // calculate bulk modulus for material m
-      const Real bmod_m = eos_v(pte_idxs(tid, mp))
-                              .BulkModulusFromDensityTemperature(
-                                  rho_pte(tid, mp), temp_pte(tid, mp), cache[mp]);
-      // add bmod contribution from material m
-      bmod_v(i) += bmod_m * vfrac_pte(tid, mp);
-      // calculate specific heat for material m
-      const Real cv_m = ev2k * eos_v(pte_idxs(tid, mp))
-                                   .SpecificHeatFromDensityTemperature(
-                                       rho_pte(tid, mp), temp_pte(tid, mp), cache[mp]);
-      // add mass weighted contribution specific heat for material m
-      cv_v(i) += cv_m * frac_mass_v(i, m);
-      // calculate gruneisen parameter for material m
-      const Real dpde_m = eos_v(pte_idxs(tid, mp))
-                              .GruneisenParamFromDensityTemperature(
-                                  rho_pte(tid, mp), temp_pte(tid, mp), cache[mp]);
-      // add gruneisen param contribution from material m
-      dpde_v(i) += dpde_m * vfrac_pte(tid, mp);
-      // optionally assign per material quantities to per material arrays
-      if (do_frac_bmod) {
-        frac_bmod_v(i, m) = bmod_m;
-      }
-      if (do_frac_cv) {
-        frac_cv_v(i, m) = cv_m;
-      }
-      if (do_frac_dpde) {
-        frac_dpde_v(i, m) = dpde_m;
-      }
-    }
-    if (do_t) {
-      temp_v(i) /= ev2k;
-    }
-    if (do_s) {
-      sie_v(i) /= mass_sum;
-    }
-    // reset mass fractions to original values if not normalized to 1
-    for (int m = 0; m < nmat; ++m) {
-      frac_mass_v(i, m) *= mass_sum;
-    }
-  };
-  Kokkos::View<int, MemoryTraits<at_int>> res("PTE::num fails");
-  Kokkos::View<int, MemoryTraits<at_int>> n_solves("PTE::num solves");
-  const std::string perf_nums =
-      "[" + std::to_string(nmat) + "," + std::to_string(ncell) + "]";
-  auto input_int_enum = static_cast<input_condition>(input_int);
-  switch (input_int_enum) {
-  case input_condition::RHO_T_INPUT: {
-    // T-rho input
-    // set frac_vol = 1/nmat
-    // set rho_i = nmat / spvol * frac_mass_i
-    // iterate PTE solver to obtain internal energies
-    // that results in the input T
-    pte_solver_scratch_size = PTESolverFixedTRequiredScratch(nmat);
-    solver_scratch = ScratchV<double>(VAWI("PTE::scratch solver"), scratch_size,
-                                      pte_solver_scratch_size);
-    const std::string rt_name = "PTE::solve (rho,T) input" + perf_nums;
-    portableFor(
-        rt_name.c_str(), 0, ncell, PORTABLE_LAMBDA(const int &iloop) {
-          // cell offset
-          const int i{offsets_v(iloop) - 1};
-          // get "thread-id" like thing with optimization
-          // for small loops
-          const int32_t token{tokens.acquire()};
-          const int32_t tid{small_loop ? iloop : token};
-          double mass_sum{0.0};
-          int npte{0};
-          // initialize values for solver / lookup
-          init_lambda(i, tid, mass_sum, npte, 1.0, 0.0, 0.0);
-          // calculate pte condition (lookup for 1 mat cell)
-          Real sie_tot_true{0.0};
-          const int neq = npte;
-          singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) +
-                                                     neq * (neq + 4) + 2 * npte);
-          if (npte > 1) {
-            // create solver lambda
-            // eos accessor
-            singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
-            PTESolverFixedT<singularity::EOSAccessor_, Real *, Real **> method(
-                npte, eos_inx, 1.0, temp_pte(tid, 0), &rho_pte(tid, 0),
-                &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0),
-                &press_pte(tid, 0), cache, &solver_scratch(tid, 0));
-            const bool res_{PTESolver(method)};
-            // calculate total internal energy
-            for (int mp = 0; mp < npte; ++mp) {
-              const int m = pte_mats(tid, mp);
-              sie_tot_true += sie_pte(tid, mp) * frac_mass_v(i, m);
-            }
-          } else {
-            // pure cell (nmat = 1)
-            // calculate sie from single eos
-            sie_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
-                                  .InternalEnergyFromDensityTemperature(
-                                      rho_pte(tid, 0), temp_pte(tid, 0), cache[0]);
-            // set total sie to material 0 value
-            sie_tot_true = sie_pte(tid, 0);
-            // set pressure
-            press_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
-                                    .PressureFromDensityTemperature(
-                                        rho_pte(tid, 0), temp_pte(tid, 0), cache[0]);
-          }
-          // total sie is known
-          sie_v(i) = sie_tot_true;
-          // assign quantities
-          final_lambda(i, tid, npte, mass_sum, 0.0, 0.0, 1.0, cache);
-          // assign max pressure
-          pmax_v(i) = press_v(i) > pmax_v(i) ? press_v(i) : pmax_v(i);
-          // release the token used for scratch arrays
-          tokens.release(token);
-        });
-    break;
-  }
-  case input_condition::RHO_P_INPUT: {
-    // rho-P input
-    // set frac_vol = 1/nmat
-    // set rho_i = nmat / spvol * frac_mass_i
-    // iterate PTE solver to obtain internal energies
-    // that results in the input P
-    pte_solver_scratch_size = PTESolverRhoTRequiredScratch(nmat);
-    solver_scratch = ScratchV<double>(VAWI("PTE::scratch solver"), scratch_size,
-                                      pte_solver_scratch_size);
-    const std::string rp_name = "PTE::solve (rho,P) input" + perf_nums;
-    portableFor(
-        rp_name.c_str(), 0, ncell, PORTABLE_LAMBDA(const int &iloop) {
-          // cell offset
-          const int i{offsets_v(iloop) - 1};
-          // get "thread-id" like thing with optimization
-          // for small loops
-          const int32_t token{tokens.acquire()};
-          const int32_t tid{small_loop ? iloop : token};
-          double mass_sum{0.0};
-          int npte{0};
-          // initialize values for solver / lookup
-          init_lambda(i, tid, mass_sum, npte, 0.0, 0.0, 1.0);
-          Real sie_tot_true{0.0};
-          const int neq = npte + 1;
-          singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) +
-                                                     neq * (neq + 4) + 2 * npte);
-          if (npte > 1) {
-            // create solver lambda
-            // eos accessor
-            singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
-            PTESolverFixedP<singularity::EOSAccessor_, Real *, Real **> method(
-                npte, eos_inx, 1.0, press_pte(tid, 0), &rho_pte(tid, 0),
-                &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0),
-                &press_pte(tid, 0), cache[0], &solver_scratch(tid, 0));
-            const bool res_{PTESolver(method)};
-            // calculate total sie
-            for (int mp = 0; mp < npte; ++mp) {
-              const int m = pte_mats(tid, mp);
-              sie_tot_true += sie_pte(tid, mp) * frac_mass_v(i, m);
-            }
-          } else {
-            // pure cell (nmat = 1)
-            // calculate sie from single eos
-            auto p_from_t = [&](const Real &t_i) {
-              return eos_v(pte_idxs(tid, 0))
-                  .PressureFromDensityTemperature(rho_pte(tid, 0), t_i, cache[0]);
-            };
-            // calculate sie root bounds
-            Real r_rho{}, r_temp{}, r_sie{}, r_press{}, r_cv{}, r_bmod{};
-            Real r_dpde{}, r_dvdt{};
-            eos_v(pte_idxs(tid, 0))
-                .ValuesAtReferenceState(r_rho, r_temp, r_sie, r_press, r_cv, r_bmod,
-                                        r_dpde, r_dvdt);
-            Real temp_i;
-            RootFinding1D::findRoot(p_from_t, press_v(i), r_temp, 1.e-10 * r_temp,
-                                    1.e10 * r_temp, 1.e-12, 1.e-12, temp_i);
-            sie_pte(tid, 0) = eos_v(pte_idxs(tid, 0))
-                                  .InternalEnergyFromDensityTemperature(rho_pte(tid, 0),
-                                                                        temp_i, cache[0]);
-            sie_tot_true = sie_pte(tid, 0);
-            // set temperature
-            temp_pte(tid, 0) = temp_i;
-          }
-          // sie calculate explicitly already
-          sie_v(i) = sie_tot_true;
-          // assign remaining outputs
-          final_lambda(i, tid, npte, mass_sum, 1.0, 0.0, 0.0, cache);
-          // assign max pressure
-          pmax_v(i) = press_v(i) > pmax_v(i) ? press_v(i) : pmax_v(i);
-          // release the token used for scratch arrays
-          tokens.release(token);
-        });
-    break;
-  }
-  case input_condition::P_T_INPUT: {
-    // P-T input
-    const int pte_solver_scratch_size = nmat * MAX_NUM_LAMBDAS;
-    solver_scratch = ScratchV<double>(VAWI("PTE::scratch solver"), scratch_size,
-                                      pte_solver_scratch_size);
-    const std::string pt_name = "PTE::solve (P,T) input" + perf_nums;
-    portableFor(
-        pt_name.c_str(), 0, ncell, PORTABLE_LAMBDA(const int &iloop) {
-          // cell offset
-          const int i{offsets_v(iloop) - 1};
-          // get "thread-id" like thing with optimization
-          // for small loops
-          const int32_t token{tokens.acquire()};
-          const int32_t tid{small_loop ? iloop : token};
-          // caching mechanism
-          singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0));
-          double mass_sum{0.0};
-          // normalize mass fractions
-          // first find the mass sum
-          // also set idxs as the decrement of the eos offsets
-          // to take into account 1 based indexing in fortran
-          for (int m = 0; m < nmat; ++m) {
-            mass_sum += frac_mass_v(i, m);
-            pte_idxs(tid, m) = eos_offsets_v(m) - 1;
-            pte_mats(tid, m) = m;
-            temp_pte(tid, m) = temp_v(i) * ev2k;
-            press_pte(tid, m) = press_v(i);
-          }
-          for (int m = 0; m < nmat; ++m) {
-            frac_mass_v(i, m) /= mass_sum;
-          }
-          // do r-e of pt for each mat
-          singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
-          Real vfrac_tot{0.0};
-          Real sie_tot{0.0};
-          for (int m = 0; m < nmat; ++m) {
-            // obtain rho and sie from P-T
-            eos_inx[m].DensityEnergyFromPressureTemperature(
-                press_pte(tid, m), temp_pte(tid, m), cache[m], rho_pte(tid, m),
-                sie_pte(tid, m));
-            // assign volume fractions
-            // this is a physical volume
-            vfrac_pte(tid, m) = frac_mass_v(i, m) / rho_pte(tid, m) * mass_sum;
-            vfrac_tot += vfrac_pte(tid, m);
-            // add internal energy component
-            sie_tot += sie_pte(tid, m) * frac_mass_v(i, m);
-          }
-          // assign volume, etc.
-          // total sie is known
-          sie_v(i) = sie_tot;
-          vol_v(i) = vfrac_tot;
-          spvol_v(i) = vol_v(i) / mass_sum;
-          for (int m = 0; m < nmat; ++m) {
-            vfrac_pte(tid, m) /= vfrac_tot;
-          }
-          // assign remaining outputs
-          final_lambda(i, tid, nmat, mass_sum, 0.0, 0.0, 0.0, cache);
-          // assign max pressure
-          pmax_v(i) = press_v(i) > pmax_v(i) ? press_v(i) : pmax_v(i);
-          // release the token used for scratch arrays
-          tokens.release(token);
-        });
-    break;
-  }
-  case input_condition::NORM_RHO_E_INPUT:
-    // rho-sie input
-    // no break so fallthrough to case 1
-  case input_condition::RHO_E_INPUT: {
-    // rho-sie input
-    pte_solver_scratch_size = PTESolverRhoTRequiredScratch(nmat);
-    solver_scratch = ScratchV<double>(VAWI("PTE::scratch solver"), scratch_size,
-                                      pte_solver_scratch_size);
-    const std::string re_name = "PTE::solve (rho,e) input" + perf_nums;
-    portableFor(
-        re_name.c_str(), 0, ncell, PORTABLE_LAMBDA(const int &iloop) {
-          // cell offset
-          const int i{offsets_v(iloop) - 1};
-          // get "thread-id" like thing with optimization
-          // for small loops
-          const int32_t token{tokens.acquire()};
-          const int32_t tid{small_loop ? iloop : token};
-          double mass_sum{0.0};
-          int npte{0};
-          // initialize values for solver / lookup
-          init_lambda(i, tid, mass_sum, npte, 0.0, 1.0, 0.0);
-          // get cache from offsets into scratch
-          const int neq = npte + 1;
-          singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) +
-                                                     neq * (neq + 4) + 2 * npte);
-          if (npte > 1) {
-            // create solver lambda
-            // eos accessor
-            singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
-            // reset inputs
-            PTESolverRhoT<singularity::EOSAccessor_, Real *, Real **> method(
-                npte, eos_inx, 1.0, sie_v(i), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
-                &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache,
-                &solver_scratch(tid, 0));
-            const bool res_{PTESolver(method)};
-          } else {
-            // pure cell (nmat = 1)
-            // calculate sie from single eos
-            Real dpdr_m, dtdr_m, dpde_m, dtde_m;
-            eos_v(pte_idxs(tid, 0))
-                .PTofRE(rho_pte(tid, 0), sie_pte(tid, 0), cache[0], press_pte(tid, 0),
-                        temp_pte(tid, 0), dpdr_m, dpde_m, dtdr_m, dtde_m);
-          }
-          // assign outputs
-          final_lambda(i, tid, npte, mass_sum, 1.0, 0.0, 1.0, cache);
-          // assign max pressure
-          pmax_v(i) = press_v(i) > pmax_v(i) ? press_v(i) : pmax_v(i);
-          // release the token used for scratch arrays
-          tokens.release(token);
-        });
-    break;
-  }
-  }
-
-  Kokkos::fence();
-  // copy results back into local values
-  // there is lots of room for performance optimization
-  // in terms of when to copy and when not necessary
-  // this is the return value (number of solve failures)
-  deep_copy(ret, res);
-  // copy pressure, this is not needed in all cases
-  deep_copy(press_hv, press_v);
-  // return max pressure, this may be needed
-  deep_copy(pmax_hv, pmax_v);
-  // I don't think the volume is necessary
-  deep_copy(vol_hv, vol_v);
-  // specific volume, copy-back not needed in all cases
-  deep_copy(spvol_hv, spvol_v);
-  // internal energy, copy-back not needed in all cases
-  deep_copy(sie_hv, sie_v);
-  // temperature, copy-back not needed in all cases
-  deep_copy(temp_hv, temp_v);
-  // bulk modulus, alwasy copy-back
-  deep_copy(bmod_hv, bmod_v);
-  // dpde, always copy-back
-  deep_copy(dpde_hv, dpde_v);
-  // specific heat, always copy-back
-  deep_copy(cv_hv, cv_v);
-  // volume fractions, always copy-back (maybe not for pure cells)
-  deep_copy(frac_vol_hv, frac_vol_v);
-  // component internal energies, always copy-back (maybe not for pure cells)
-  deep_copy(frac_ie_hv, frac_ie_v);
-  // optionally copy-back the component bmod, dpde, and cv
-  if (do_frac_bmod) {
-    deep_copy(frac_bmod_hv, frac_bmod_v);
-  }
-  if (do_frac_dpde) {
-    deep_copy(frac_dpde_hv, frac_dpde_v);
-  }
-  if (do_frac_cv) {
-    deep_copy(frac_cv_hv, frac_cv_v);
-  }
-#endif // PORTABILITY_STRATEGY_KOKKOS
-  return ret;
-}
-
 int finalize_sg_eos(const int nmat, EOS *&eos, const int own_kokkos) {
   {
     for (int i = 0; i < nmat; ++i)
diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp
index 2382909666..1007e02711 100644
--- a/singularity-eos/eos/singularity_eos.hpp
+++ b/singularity-eos/eos/singularity_eos.hpp
@@ -15,7 +15,6 @@
 #ifndef _SINGULARITY_EOS_EOS_SINGULARITY_EOS_HPP_
 #define _SINGULARITY_EOS_EOS_SINGULARITY_EOS_HPP_
 
-#include <singularity-eos/closure/mixed_cell_models.hpp>
 #include <singularity-eos/eos/eos.hpp>
 
 using singularity::EOS;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index c0fa96f8dd..cc488b7fc2 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -42,11 +42,21 @@ if(SINGULARITY_USE_EOSPAC AND SINGULARITY_TEST_SESAME)
                         PRIVATE Catch2::Catch2 singularity-eos::singularity-eos)
 endif()
 
-if(SINGULARITY_BUILD_CLOSURE AND SINGULARITY_USE_SPINER)
-  add_executable(test_pte test_pte.cpp)
-  target_link_libraries(test_pte PRIVATE Catch2::Catch2
+if(SINGULARITY_BUILD_CLOSURE)
+  if(SINGULARITY_USE_SPINER)
+    add_executable(test_pte test_pte.cpp)
+    target_link_libraries(test_pte PRIVATE Catch2::Catch2
                                          singularity-eos::singularity-eos)
-  add_test(pte test_pte)
+    target_include_directories(test_pte PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+    add_test(pte test_pte)
+  endif()
+  if(SINGULARITY_USE_KOKKOS AND SINGULARITY_USE_FORTRAN)
+    add_executable(test_get_sg_eos test_get_sg_eos.cpp)
+    target_link_libraries(test_get_sg_eos PRIVATE Catch2::Catch2
+                                         singularity-eos::singularity-eos)
+    target_include_directories(test_get_sg_eos PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+    add_test(get_sg_eos test_get_sg_eos)
+  endif()
 endif()
 
 if(SINGULARITY_TEST_STELLAR_COLLAPSE)
diff --git a/test/pte_test_utils.hpp b/test/pte_test_utils.hpp
new file mode 100644
index 0000000000..d61b09747e
--- /dev/null
+++ b/test/pte_test_utils.hpp
@@ -0,0 +1,91 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#ifndef _SINGULARITY_EOS_TEST_PTE_TEST_UTILS_
+#define _SINGULARITY_EOS_TEST_PTE_TEST_UTILS_
+
+#include <stdlib.h>
+
+#include <ports-of-call/portability.hpp>
+#include <ports-of-call/portable_arrays.hpp>
+#include <singularity-eos/eos/eos.hpp>
+
+constexpr int NMAT = 3;
+constexpr int NTRIAL = 100;
+constexpr int NPTS = NTRIAL * NMAT;
+constexpr int HIST_SIZE = 10;
+
+template <typename T>
+class LinearIndexer {
+ public:
+  PORTABLE_FUNCTION LinearIndexer() = default;
+  LinearIndexer(const T &t) : data_(t) {}
+  PORTABLE_INLINE_FUNCTION
+  auto &operator[](const int i) const { return data_(i); }
+
+ private:
+  T data_;
+};
+
+template <typename T>
+class Indexer2D {
+ public:
+  PORTABLE_FUNCTION Indexer2D() = default;
+  PORTABLE_FUNCTION Indexer2D(const int j, const T &t) : j_(j), data_(t) {}
+  Indexer2D(const int j, const T &&t) = delete; // prevents r-value binding
+  PORTABLE_INLINE_FUNCTION
+  auto &operator[](const int i) const { return data_(j_, i); }
+
+ private:
+  const int j_;
+  const T &data_;
+};
+
+template <typename T>
+inline void set_eos(T *eos) {
+  singularity::EOS gr = singularity::Gruneisen(394000.0, 1.489, 0.0, 0.0, 2.02, 0.47,
+                                               8.93, 297.0, 1.0e6, 0.383e7);
+  singularity::EOS dr = singularity::DavisReactants(
+      1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10);
+  singularity::EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419,
+                                                   3.2e10, 0.001072e10, 0.0);
+  eos[0] = gr.GetOnDevice();
+  eos[1] = dr.GetOnDevice();
+  eos[2] = dp.GetOnDevice();
+  return;
+}
+
+template <typename RealIndexer, typename EOSIndexer>
+inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie,
+                      RealIndexer &&temp, EOSIndexer &&eos) {
+
+  rho[0] = 8.93;
+  rho[1] = 1.89;
+  rho[2] = 2.5;
+
+  Real vsum = 0.;
+  for (int i = 0; i < NMAT; i++) {
+    temp[i] = 600.0;
+    sie[i] = eos[i].InternalEnergyFromDensityTemperature(rho[i], temp[i]);
+    vfrac[i] = rand() / (1.0 * RAND_MAX);
+    vsum += vfrac[i];
+  }
+
+  for (int i = 0; i < NMAT; i++)
+    vfrac[i] *= 1.0 / vsum;
+
+  return;
+}
+
+#endif // _SINGULARITY_EOS_TEST_PTE_TEST_UTILS_
diff --git a/test/test_eos_noble_abel.cpp b/test/test_eos_noble_abel.cpp
index eab2e3e2d9..5d6b4264bd 100644
--- a/test/test_eos_noble_abel.cpp
+++ b/test/test_eos_noble_abel.cpp
@@ -532,8 +532,8 @@ SCENARIO("Test NobleAbel Entropy Calls", "[NobleAbel][NobleAbel5]") {
 #endif // PORTABILITY_STRATEGY_KOKKOS
 
       // Gold standard values for a subset of lookups
-      constexpr std::array<Real, num> pressure_true{3.4450499999999993e+07};
-      constexpr std::array<Real, num> temperature_true{1.5320999999999999e+03};
+      // constexpr std::array<Real, num> pressure_true{3.4450499999999993e+07};
+      // constexpr std::array<Real, num> temperature_true{1.5320999999999999e+03};
       constexpr std::array<Real, num> entropy_true{-2.0177705273504910e+08};
 
 #ifdef PORTABILITY_STRATEGY_KOKKOS
diff --git a/test/test_eos_stiff.cpp b/test/test_eos_stiff.cpp
index cdb1f45e8c..548af96a93 100644
--- a/test/test_eos_stiff.cpp
+++ b/test/test_eos_stiff.cpp
@@ -534,8 +534,6 @@ SCENARIO("Test Stiff Gas Entropy Calls", "[StiffGas][StiffGas5]") {
 #endif // PORTABILITY_STRATEGY_KOKKOS
 
       // Gold standard values for a subset of lookups
-      constexpr std::array<Real, num> pressure_true{2.0265000000000000e+06};
-      constexpr std::array<Real, num> temperature_true{4.0000000000000000e+02};
       constexpr std::array<Real, num> entropy_true{-2.0044437857420778e+08};
 
 #ifdef PORTABILITY_STRATEGY_KOKKOS
diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp
index 1ed9499c2f..fd8d3a30c1 100644
--- a/test/test_eos_unit.cpp
+++ b/test/test_eos_unit.cpp
@@ -892,13 +892,14 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS],[DependsRhoT][EOSPAC]")
 
         int nw_ie{0};
 #ifdef PORTABILITY_STRATEGY_KOKKOS
+        Kokkos::fence();
         using atomic_view = Kokkos::MemoryTraits<Kokkos::Atomic>;
         Kokkos::View<int, atomic_view> n_wrong_ie("wrong_ie");
 #else
         PortableMDArray<int> n_wrong_ie(&nw_ie, 1);
 #endif
         portableFor(
-            "calc ie's", 0, 100, PORTABLE_LAMBDA(const int &i) {
+            "calc ie's steel", 0, 100, PORTABLE_LAMBDA(const int &i) {
               const double ie = steelEOS.InternalEnergyFromDensityTemperature(1e0, 1e6);
               if (!isClose(ie, sie_pac)) n_wrong_ie() += 1;
             });
@@ -943,13 +944,14 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS],[DependsRhoT][EOSPAC]")
       REQUIRE(isClose(T, T_pac));
     }
     THEN("Quantities of rho and sie look sane on both host and device") {
-      Real gm1_host = airEOS_host.GruneisenParamFromDensityInternalEnergy(rho, sie);
-      Real T_host = airEOS_host.TemperatureFromDensityInternalEnergy(rho, sie);
-      Real cv_host = airEOS_host.SpecificHeatFromDensityInternalEnergy(rho, sie);
+      const Real gm1_host = airEOS_host.GruneisenParamFromDensityInternalEnergy(rho, sie);
+      const Real T_host = airEOS_host.TemperatureFromDensityInternalEnergy(rho, sie);
+      const Real cv_host = airEOS_host.SpecificHeatFromDensityInternalEnergy(rho, sie);
 
       int nw_gm1{0}, nw_cv{0};
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       using atomic_view = Kokkos::MemoryTraits<Kokkos::Atomic>;
+      Kokkos::fence();
       Kokkos::View<int, atomic_view> n_wrong_gm1("wrong_gm1");
       Kokkos::View<int, atomic_view> n_wrong_cv("wrong_cv");
 #else
@@ -958,8 +960,8 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS],[DependsRhoT][EOSPAC]")
 #endif
       portableFor(
           "calc gm1 and cv", 0, 100, PORTABLE_LAMBDA(const int &i) {
-            const double gm1 = airEOS.GruneisenParamFromDensityInternalEnergy(rho, sie);
-            const double cv = airEOS.SpecificHeatFromDensityInternalEnergy(rho, sie);
+            const Real gm1 = airEOS.GruneisenParamFromDensityInternalEnergy(rho, sie);
+            const Real cv = airEOS.SpecificHeatFromDensityInternalEnergy(rho, sie);
             if (!isClose(gm1, gm1_host)) n_wrong_gm1() += 1;
             if (!isClose(cv, cv_host)) n_wrong_cv() += 1;
           });
@@ -1030,6 +1032,7 @@ SCENARIO("SpinerEOS depends on rho and sie", "[SpinerEOS],[DependsRhoSie]") {
 
       int nw_ie2{0}, nw_te2{0};
 #ifdef PORTABILITY_STRATEGY_KOKKOS
+      Kokkos::fence();
       using atomic_view = Kokkos::MemoryTraits<Kokkos::Atomic>;
       Kokkos::View<int, atomic_view> n_wrong_ie("wrong_ie");
       Kokkos::View<int, atomic_view> n_wrong_te("wrong_te");
@@ -1038,7 +1041,7 @@ SCENARIO("SpinerEOS depends on rho and sie", "[SpinerEOS],[DependsRhoSie]") {
       PortableMDArray<int> n_wrong_te(&nw_te2, 1);
 #endif
       portableFor(
-          "calc ie's", 0, 100, PORTABLE_LAMBDA(const int &i) {
+          "calc ie's steel 2", 0, 100, PORTABLE_LAMBDA(const int &i) {
             const Real ie{steelEOS.InternalEnergyFromDensityTemperature(1e0, 1e6)};
             const Real te{steelEOS.TemperatureFromDensityInternalEnergy(1e0, 1e12)};
             if (!isClose(ie, 4.96415e13)) n_wrong_ie() += 1;
@@ -1082,13 +1085,14 @@ SCENARIO("EOS Builder and SpinerEOS",
         THEN("EOS calls match raw access") {
           int nw_bm{0};
 #ifdef PORTABILITY_STRATEGY_KOKKOS
+          Kokkos::fence();
           using atomic_view = Kokkos::MemoryTraits<Kokkos::Atomic>;
           Kokkos::View<int, atomic_view> n_wrong_bm("wrong_bm");
 #else
           PortableMDArray<int> n_wrong_bm(&nw_bm, 1);
 #endif
           portableFor(
-              "calc ie's", 0, 100, PORTABLE_LAMBDA(const int &i) {
+              "calc ie's steel 3", 0, 100, PORTABLE_LAMBDA(const int &i) {
                 const Real bm{eosDevice.BulkModulusFromDensityTemperature(1e0, 1e6)};
                 if (!isClose(bm, 2.55268e13)) n_wrong_bm() += 1;
               });
diff --git a/test/test_get_sg_eos.cpp b/test/test_get_sg_eos.cpp
new file mode 100644
index 0000000000..85c664ddf6
--- /dev/null
+++ b/test/test_get_sg_eos.cpp
@@ -0,0 +1,154 @@
+//------------------------------------------------------------------------------
+// © 2021-2023. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#include <iostream>
+#include <memory>
+#include <stdlib.h>
+
+#include <ports-of-call/portability.hpp>
+#include <pte_test_utils.hpp>
+#include <singularity-eos/closure/mixed_cell_models.hpp>
+#include <singularity-eos/eos/eos.hpp>
+#include <singularity-eos/eos/singularity_eos.hpp>
+
+using namespace singularity;
+
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+// TODO DAH: when the get_sg_eos function is moved out of sg,
+// this function will have to change how the PTE solutions
+// are obtained.
+int run_sg_get_eos_tests() {
+  int nfails = 0;
+  // initialize inputs outputs
+  static constexpr const double ev2k = 1.160451930280894026e4;
+  EOS eoss[NMAT];
+  set_eos(eoss);
+  // set volume fractions to equipartition
+  // all variables needed to check the pte solve
+  Real mfrac[NMAT];
+  Real P_true = 5.e10, T_true = 800.0;
+  Real T_true_ev = T_true / ev2k;
+  Real mass_tot = 0.0, sie_tot_true = 0.0, v_true = 1.0, rho_tot, spvol;
+  Real pmax, bmod, dpde, cv;
+  mfrac[0] = 0.67;
+  mfrac[1] = 0.14;
+  mfrac[2] = 0.19;
+  for (int i = 0; i < NMAT; ++i) {
+    mass_tot += mfrac[i];
+  }
+  rho_tot = mass_tot / v_true;
+  spvol = 1.0 / rho_tot;
+  int cell_offset = 1;
+  int eos_offset[NMAT];
+  for (int m = 0; m < NMAT; ++m) {
+    eos_offset[m] = m + 1;
+  }
+  // do a rho/e(P,T) solve
+  Real vfrac_true[NMAT], ie_true[NMAT];
+  get_sg_eos(NMAT, 1, 1, -1, eos_offset, eoss, &cell_offset, &P_true, &pmax, &v_true,
+             &spvol, &sie_tot_true, &T_true_ev, &bmod, &dpde, &cv, mfrac, vfrac_true,
+             ie_true, nullptr, nullptr, nullptr);
+  Real sie_tot_check = 0.0;
+  for (int m = 0; m < NMAT; ++m) {
+    const Real r_m = mfrac[m] / vfrac_true[m];
+    const Real sie_m = ie_true[m] / mfrac[m];
+    const Real p_eos = eoss[m].PressureFromDensityInternalEnergy(r_m, sie_m);
+    const Real p_resid = std::abs(P_true - p_eos) / P_true;
+    const Real t_eos = eoss[m].TemperatureFromDensityInternalEnergy(r_m, sie_m);
+    const Real t_resid = std::abs(T_true - t_eos) / T_true;
+    // check internal consitency of P-T input function
+    if (t_resid > 1.e-5 || p_resid > 1.e-5) {
+      printf("P-T: t_resid: %e | p_resid: %e\n", t_resid, p_resid);
+      nfails += 1;
+    }
+    sie_tot_check += ie_true[m];
+  }
+  sie_tot_check /= mass_tot;
+  // further check of internal consintency of P-T input function
+  if (std::abs(sie_tot_check - sie_tot_true) / std::abs(sie_tot_true) > 1.e-5) {
+    printf("P-T: sie_tot_true: %e | sie_tot_check: %e\n", sie_tot_true, sie_tot_check);
+    nfails += 1;
+  }
+  // obtain converged and consistent PTE solution
+  // do rho-T input solve
+  Real p_check, vfrac_check[NMAT], ie_check[NMAT];
+  get_sg_eos(NMAT, 1, 1, -3, eos_offset, eoss, &cell_offset, &p_check, &pmax, &v_true,
+             &spvol, &sie_tot_check, &T_true_ev, &bmod, &dpde, &cv, mfrac, vfrac_check,
+             ie_check, nullptr, nullptr, nullptr);
+  // check output pressure and sie, indicate failure if relative err is too large
+  if (std::abs(P_true - p_check) / std::abs(P_true) > 1.e-5 ||
+      std::abs(sie_tot_true - sie_tot_check) / std::abs(sie_tot_true) > 1.e-5) {
+    printf("r-T: p_true: %e | p_check: %e\n", P_true, p_check);
+    printf("r-T: sie_tot_true: %e | sie_tot_check: %e\n", sie_tot_true, sie_tot_check);
+    nfails += 1;
+  }
+  Real max_vfrac_resid = 0.0;
+  Real max_sie_resid = 0.0;
+  for (int m = 0; m < NMAT; ++m) {
+    max_vfrac_resid = std::max(max_vfrac_resid, std::abs(vfrac_true[m] - vfrac_check[m]) /
+                                                    std::abs(vfrac_true[m]));
+    max_sie_resid = std::max(max_sie_resid,
+                             std::abs(ie_true[m] - ie_check[m]) / std::abs(ie_true[m]));
+  }
+  if (max_vfrac_resid > 1.e-5 || max_sie_resid > 1.e-5) {
+    printf("r-T: vr: %e | sr: %e\n", max_vfrac_resid, max_sie_resid);
+    nfails += 1;
+  }
+  // do rho-P input solve
+  Real t_check;
+  get_sg_eos(NMAT, 1, 1, -2, eos_offset, eoss, &cell_offset, &P_true, &pmax, &v_true,
+             &spvol, &sie_tot_check, &t_check, &bmod, &dpde, &cv, mfrac, vfrac_check,
+             ie_check, nullptr, nullptr, nullptr);
+  // check output temperature and sie, indicate failure if relative err is too large
+  if (std::abs(T_true_ev - t_check) / std::abs(T_true_ev) > 1.e-5 ||
+      std::abs(sie_tot_true - sie_tot_check) / std::abs(sie_tot_true) > 1.e-5) {
+    printf("p-T: t_true: %e | t_check: %e\n", T_true_ev, t_check);
+    printf("p-T: sie_tot_true: %e | sie_tot_check: %e\n", sie_tot_true, sie_tot_check);
+    nfails += 1;
+  }
+  max_vfrac_resid = 0.0;
+  max_sie_resid = 0.0;
+  for (int m = 0; m < NMAT; ++m) {
+    max_vfrac_resid = std::max(max_vfrac_resid, std::abs(vfrac_true[m] - vfrac_check[m]) /
+                                                    std::abs(vfrac_true[m]));
+    max_sie_resid = std::max(max_sie_resid,
+                             std::abs(ie_true[m] - ie_check[m]) / std::abs(ie_true[m]));
+  }
+  if (max_vfrac_resid > 1.e-5 || max_sie_resid > 1.e-5) {
+    printf("p-T: vr: %e | sr: %e\n", max_vfrac_resid, max_sie_resid);
+    nfails += 1;
+  }
+  return nfails;
+}
+#endif
+
+int main(int argc, char *argv[]) {
+
+  int nfails_get_sg_eos = 1;
+#ifdef PORTABILITY_STRATEGY_KOKKOS
+  Kokkos::initialize();
+  {
+    // run the various solver tests:
+    // since this uses the sg_get_eos, only run
+    // if kokkos enable since that function requires it
+    // to run the solvers
+    nfails_get_sg_eos = run_sg_get_eos_tests();
+    if (nfails_get_sg_eos > 0) {
+      printf("nfails of fixed T/P solvers = %i\n", nfails_get_sg_eos);
+    }
+  }
+  Kokkos::finalize();
+#endif // PORTABILITY_STRATEGY_KOKKOS
+  return nfails_get_sg_eos;
+}
diff --git a/test/test_pte.cpp b/test/test_pte.cpp
index 5a8034a4e4..bfed57224f 100644
--- a/test/test_pte.cpp
+++ b/test/test_pte.cpp
@@ -21,203 +21,23 @@
 
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_arrays.hpp>
+#include <pte_test_utils.hpp>
 #include <singularity-eos/closure/mixed_cell_models.hpp>
 #include <singularity-eos/eos/eos.hpp>
-#include <singularity-eos/eos/singularity_eos.hpp>
 #include <spiner/databox.hpp>
 
-constexpr int NMAT = 3;
-constexpr int NTRIAL = 100;
-constexpr int NPTS = NTRIAL * NMAT;
-constexpr int HIST_SIZE = 10;
-
 using namespace singularity;
 
 using DataBox = Spiner::DataBox<Real>;
 
-template <typename T>
-class LinearIndexer {
- public:
-  PORTABLE_FUNCTION LinearIndexer() = default;
-  LinearIndexer(const T &t) : data_(t) {}
-  PORTABLE_INLINE_FUNCTION
-  auto &operator[](const int i) const { return data_(i); }
-
- private:
-  T data_;
-};
-
-template <typename T>
-class Indexer2D {
- public:
-  PORTABLE_FUNCTION Indexer2D() = default;
-  PORTABLE_FUNCTION Indexer2D(const int j, const T &t) : j_(j), data_(t) {}
-  Indexer2D(const int j, const T &&t) = delete; // prevents r-value binding
-  PORTABLE_INLINE_FUNCTION
-  auto &operator[](const int i) const { return data_(j_, i); }
-
- private:
-  const int j_;
-  const T &data_;
-};
-
-template <typename EOSIndexer>
-void set_eos(EOSIndexer &&eos) {
-  eos[0] = Gruneisen(394000.0, 1.489, 0.0, 0.0, 2.02, 0.47, 8.93, 297.0, 1.0e6, 0.383e7);
-  eos[1] = DavisReactants(1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0,
-                          0.4265, 0.001074e10);
-  eos[2] =
-      DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, 3.2e10, 0.001072e10, 0.0);
-}
-
-template <typename RealIndexer, typename EOSIndexer>
-void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie,
-               RealIndexer &&temp, EOSIndexer &&eos) {
-
-  rho[0] = 8.93;
-  rho[1] = 1.89;
-  rho[2] = 2.5;
-
-  Real vsum = 0.;
-  for (int i = 0; i < NMAT; i++) {
-    temp[i] = 600.0;
-    sie[i] = eos[i].InternalEnergyFromDensityTemperature(rho[i], temp[i]);
-    vfrac[i] = rand() / (1.0 * RAND_MAX);
-    vsum += vfrac[i];
-  }
-
-  for (int i = 0; i < NMAT; i++)
-    vfrac[i] *= 1.0 / vsum;
-}
-
-#ifdef PORTABILITY_STRATEGY_KOKKOS
-// TODO DAH: when the get_sg_eos function is moved out of sg,
-// this function will have to change how the PTE solutions
-// are obtained.
-int run_sg_get_eos_tests() {
-  int nfails = 0;
-  // initialize inputs outputs
-  static constexpr const double ev2k = 1.160451930280894026e4;
-  EOS eoss[NMAT];
-  set_eos(eoss);
-  // set volume fractions to equipartition
-  // all variables needed to check the pte solve
-  Real mfrac[NMAT];
-  Real P_true = 5.e10, T_true = 800.0;
-  Real T_true_ev = T_true / ev2k;
-  Real mass_tot = 0.0, sie_tot_true = 0.0, v_true = 1.0, rho_tot, spvol;
-  Real pmax, bmod, dpde, cv;
-  mfrac[0] = 0.67;
-  mfrac[1] = 0.14;
-  mfrac[2] = 0.19;
-  for (int i = 0; i < NMAT; ++i) {
-    mass_tot += mfrac[i];
-  }
-  rho_tot = mass_tot / v_true;
-  spvol = 1.0 / rho_tot;
-  int cell_offset = 1;
-  int eos_offset[NMAT];
-  for (int m = 0; m < NMAT; ++m) {
-    eos_offset[m] = m + 1;
-  }
-  // do a rho/e(P,T) solve
-  Real vfrac_true[NMAT], ie_true[NMAT];
-  get_sg_eos(NMAT, 1, 1, -1, eos_offset, eoss, &cell_offset, &P_true, &pmax, &v_true,
-             &spvol, &sie_tot_true, &T_true_ev, &bmod, &dpde, &cv, mfrac, vfrac_true,
-             ie_true, nullptr, nullptr, nullptr);
-  Real sie_tot_check = 0.0;
-  for (int m = 0; m < NMAT; ++m) {
-    const Real r_m = mfrac[m] / vfrac_true[m];
-    const Real sie_m = ie_true[m] / mfrac[m];
-    const Real p_eos = eoss[m].PressureFromDensityInternalEnergy(r_m, sie_m);
-    const Real p_resid = std::abs(P_true - p_eos) / P_true;
-    const Real t_eos = eoss[m].TemperatureFromDensityInternalEnergy(r_m, sie_m);
-    const Real t_resid = std::abs(T_true - t_eos) / T_true;
-    // check internal consitency of P-T input function
-    if (t_resid > 1.e-5 || p_resid > 1.e-5) {
-      printf("P-T: t_resid: %e | p_resid: %e\n", t_resid, p_resid);
-      nfails += 1;
-    }
-    sie_tot_check += ie_true[m];
-  }
-  sie_tot_check /= mass_tot;
-  // further check of internal consintency of P-T input function
-  if (std::abs(sie_tot_check - sie_tot_true) / std::abs(sie_tot_true) > 1.e-5) {
-    printf("P-T: sie_tot_true: %e | sie_tot_check: %e\n", sie_tot_true, sie_tot_check);
-    nfails += 1;
-  }
-  // obtain converged and consistent PTE solution
-  // do rho-T input solve
-  Real p_check, vfrac_check[NMAT], ie_check[NMAT];
-  get_sg_eos(NMAT, 1, 1, -3, eos_offset, eoss, &cell_offset, &p_check, &pmax, &v_true,
-             &spvol, &sie_tot_check, &T_true_ev, &bmod, &dpde, &cv, mfrac, vfrac_check,
-             ie_check, nullptr, nullptr, nullptr);
-  // check output pressure and sie, indicate failure if relative err is too large
-  if (std::abs(P_true - p_check) / std::abs(P_true) > 1.e-5 ||
-      std::abs(sie_tot_true - sie_tot_check) / std::abs(sie_tot_true) > 1.e-5) {
-    printf("r-T: p_true: %e | p_check: %e\n", P_true, p_check);
-    printf("r-T: sie_tot_true: %e | sie_tot_check: %e\n", sie_tot_true, sie_tot_check);
-    nfails += 1;
-  }
-  Real max_vfrac_resid = 0.0;
-  Real max_sie_resid = 0.0;
-  for (int m = 0; m < NMAT; ++m) {
-    max_vfrac_resid = std::max(max_vfrac_resid, std::abs(vfrac_true[m] - vfrac_check[m]) /
-                                                    std::abs(vfrac_true[m]));
-    max_sie_resid = std::max(max_sie_resid,
-                             std::abs(ie_true[m] - ie_check[m]) / std::abs(ie_true[m]));
-  }
-  if (max_vfrac_resid > 1.e-5 || max_sie_resid > 1.e-5) {
-    printf("r-T: vr: %e | sr: %e\n", max_vfrac_resid, max_sie_resid);
-    nfails += 1;
-  }
-  // do rho-P input solve
-  Real t_check;
-  get_sg_eos(NMAT, 1, 1, -2, eos_offset, eoss, &cell_offset, &P_true, &pmax, &v_true,
-             &spvol, &sie_tot_check, &t_check, &bmod, &dpde, &cv, mfrac, vfrac_check,
-             ie_check, nullptr, nullptr, nullptr);
-  // check output temperature and sie, indicate failure if relative err is too large
-  if (std::abs(T_true_ev - t_check) / std::abs(T_true_ev) > 1.e-5 ||
-      std::abs(sie_tot_true - sie_tot_check) / std::abs(sie_tot_true) > 1.e-5) {
-    printf("p-T: t_true: %e | t_check: %e\n", T_true_ev, t_check);
-    printf("p-T: sie_tot_true: %e | sie_tot_check: %e\n", sie_tot_true, sie_tot_check);
-    nfails += 1;
-  }
-  max_vfrac_resid = 0.0;
-  max_sie_resid = 0.0;
-  for (int m = 0; m < NMAT; ++m) {
-    max_vfrac_resid = std::max(max_vfrac_resid, std::abs(vfrac_true[m] - vfrac_check[m]) /
-                                                    std::abs(vfrac_true[m]));
-    max_sie_resid = std::max(max_sie_resid,
-                             std::abs(ie_true[m] - ie_check[m]) / std::abs(ie_true[m]));
-  }
-  if (max_vfrac_resid > 1.e-5 || max_sie_resid > 1.e-5) {
-    printf("p-T: vr: %e | sr: %e\n", max_vfrac_resid, max_sie_resid);
-    nfails += 1;
-  }
-  return nfails;
-}
-#endif
-
 int main(int argc, char *argv[]) {
 
   int nsuccess = 0;
-  int nfails_get_sg_eos = 0;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
   Kokkos::initialize();
 #endif
   {
-    // run the various solver tests:
-    // since this uses the sg_get_eos, only run
-    // if kokkos enable since that function requires it
-    // to run the solvers
-#ifdef PORTABILITY_STRATEGY_KOKKOS
-    nfails_get_sg_eos = run_sg_get_eos_tests();
-    if (nfails_get_sg_eos > 0) {
-      printf("nfails of fixed T/P solvers = %i\n", nfails_get_sg_eos);
-    }
-#endif // PORTABILITY_STRATEGY_KOKKOS
-       // EOS
+    // EOS
 #ifdef PORTABILITY_STRATEGY_KOKKOS
     Kokkos::View<EOS *> eos_v("eos", NMAT);
     auto eos_hv = Kokkos::create_mirror_view(eos_v);
@@ -228,7 +48,7 @@ int main(int argc, char *argv[]) {
 #endif
 
     LinearIndexer<decltype(eos_hv)> eos_h(eos_hv);
-    set_eos(eos_h);
+    set_eos(eos_hv.data());
 
 #ifdef PORTABILITY_STRATEGY_KOKKOS
     Kokkos::deep_copy(eos_v, eos_hv);
@@ -384,9 +204,5 @@ int main(int argc, char *argv[]) {
 #endif
 
   // poor-man's ctest integration
-  if (nsuccess >= 0.5 * NTRIAL) {
-    return 0 + nfails_get_sg_eos; // exit success
-  } else {
-    return 1 + nfails_get_sg_eos; // exit failure
-  }
+  return (nsuccess >= 0.5 * NTRIAL) ? 0 : 1;
 }
diff --git a/utils/cuda_compatibility.patch b/utils/cuda_compatibility.patch
index 3852721321..7c415277ff 100644
--- a/utils/cuda_compatibility.patch
+++ b/utils/cuda_compatibility.patch
@@ -1,7 +1,16 @@
 diff --git a/include/mpark/variant.hpp b/include/mpark/variant.hpp
-index 2fb2ac549..aeef9dcdd 100644
+index 2fb2ac549..76f58530a 100644
 --- a/include/mpark/variant.hpp
 +++ b/include/mpark/variant.hpp
+@@ -245,7 +245,7 @@ namespace mpark {
+   };
+ 
+   [[noreturn]] inline void throw_bad_variant_access() {
+-#ifdef MPARK_EXCEPTIONS
++#if defined(MPARK_EXCEPTIONS) && !defined(__CUDACC__)
+     throw bad_variant_access{};
+ #else
+     std::terminate();
 @@ -1021,9 +1021,15 @@ namespace mpark {
      friend struct access::recursive_union;                                 \
    }
@@ -19,7 +28,13 @@ index 2fb2ac549..aeef9dcdd 100644
                                    ~recursive_union() {});
      MPARK_VARIANT_RECURSIVE_UNION(Trait::Unavailable,
                                    ~recursive_union() = delete;);
-@@ -1133,9 +1139,9 @@ namespace mpark {
+@@ -1128,14 +1134,14 @@ namespace mpark {
+     MPARK_VARIANT_DESTRUCTOR(
+         Trait::TriviallyAvailable,
+         ~destructor() = default;,
+-        inline void destroy() noexcept {
++        inline constexpr void destroy() noexcept {
+           this->index_ = static_cast<index_t<Ts...>>(-1);
          });
  
      MPARK_VARIANT_DESTRUCTOR(
@@ -57,20 +72,15 @@ index 2fb2ac549..aeef9dcdd 100644
              lib::all<std::is_nothrow_move_constructible<Ts>::value...>::value)
              : move_constructor(valueless_t{}) {
            this->generic_construct(*this, lib::move(that));
-@@ -1256,11 +1262,11 @@ namespace mpark {
-   }
- 
-     MPARK_VARIANT_COPY_CONSTRUCTOR(
--        Trait::TriviallyAvailable,
-+        Trait::TriviallyAvailable,
-         copy_constructor(const copy_constructor &that) = default;);
+@@ -1261,7 +1267,7 @@ namespace mpark {
  
      MPARK_VARIANT_COPY_CONSTRUCTOR(
--        Trait::Available,
-+        Trait::Available,V_GPU_FUNCTION
-         copy_constructor(const copy_constructor &that)
+         Trait::Available,
+-        copy_constructor(const copy_constructor &that)
++        constexpr copy_constructor(const copy_constructor &that)
              : copy_constructor(valueless_t{}) {
            this->generic_construct(*this, that);
+         });
 @@ -1281,7 +1287,7 @@ namespace mpark {
        using super::operator=;
  
@@ -98,20 +108,15 @@ index 2fb2ac549..aeef9dcdd 100644
          if (this->valueless_by_exception() && that.valueless_by_exception()) {
            // do nothing.
          } else if (that.valueless_by_exception()) {
-@@ -1376,11 +1382,11 @@ namespace mpark {
-   }
- 
-     MPARK_VARIANT_MOVE_ASSIGNMENT(
--        Trait::TriviallyAvailable,
-+        Trait::TriviallyAvailable,
-         move_assignment &operator=(move_assignment &&that) = default;);
+@@ -1381,7 +1387,7 @@ namespace mpark {
  
      MPARK_VARIANT_MOVE_ASSIGNMENT(
--        Trait::Available,
-+        Trait::Available,V_GPU_FUNCTION
-         move_assignment &
+         Trait::Available,
+-        move_assignment &
++        constexpr move_assignment &
          operator=(move_assignment &&that) noexcept(
              lib::all<(std::is_nothrow_move_constructible<Ts>::value &&
+                       std::is_nothrow_move_assignable<Ts>::value)...>::value) {
 @@ -1394,6 +1400,7 @@ namespace mpark {
          move_assignment &operator=(move_assignment &&) = delete;);
  
@@ -129,14 +134,95 @@ index 2fb2ac549..aeef9dcdd 100644
            this->generic_assign(that);
            return *this;
          });
+@@ -1447,7 +1454,7 @@ namespace mpark {
+       impl &operator=(impl &&) = default;
+ 
+       template <std::size_t I, typename Arg>
+-      inline void assign(Arg &&arg) {
++      inline constexpr void assign(Arg &&arg) {
+         this->assign_alt(access::base::get_alt<I>(*this),
+                          lib::forward<Arg>(arg));
+       }
+@@ -1499,7 +1506,7 @@ namespace mpark {
+ #ifndef MPARK_GENERIC_LAMBDAS
+       struct swapper {
+         template <typename ThisAlt, typename ThatAlt>
+-        inline void operator()(ThisAlt &this_alt, ThatAlt &that_alt) const {
++        inline constexpr void operator()(ThisAlt &this_alt, ThatAlt &that_alt) const {
+           using std::swap;
+           swap(this_alt.value, that_alt.value);
+         }
 @@ -1694,8 +1701,8 @@ namespace mpark {
  
      ~variant() = default;
  
 -    variant &operator=(const variant &) = default;
 -    variant &operator=(variant &&) = default;
-+    /*V_GPU_FUNCTION*/constexpr variant &operator=(const variant &) = default;
-+    /*V_GPU_FUNCTION*/constexpr variant &operator=(variant &&) = default;
++    constexpr variant &operator=(const variant &) = default;
++    constexpr variant &operator=(variant &&) = default;
  
      template <typename Arg,
                lib::enable_if_t<!std::is_same<lib::decay_t<Arg>, variant>::value,
+@@ -1705,7 +1712,7 @@ namespace mpark {
+               lib::enable_if_t<(std::is_assignable<T &, Arg>::value &&
+                                 std::is_constructible<T, Arg>::value),
+                                int> = 0>
+-    inline variant &operator=(Arg &&arg) noexcept(
++    inline constexpr variant &operator=(Arg &&arg) noexcept(
+         (std::is_nothrow_assignable<T &, Arg>::value &&
+          std::is_nothrow_constructible<T, Arg>::value)) {
+       impl_.template assign<I>(lib::forward<Arg>(arg));
+@@ -1717,7 +1724,7 @@ namespace mpark {
+         typename... Args,
+         typename T = lib::type_pack_element_t<I, Ts...>,
+         lib::enable_if_t<std::is_constructible<T, Args...>::value, int> = 0>
+-    inline T &emplace(Args &&... args) {
++    inline constexpr T &emplace(Args &&... args) {
+       return impl_.template emplace<I>(lib::forward<Args>(args)...);
+     }
+ 
+@@ -1730,7 +1737,7 @@ namespace mpark {
+                                                std::initializer_list<Up> &,
+                                                Args...>::value,
+                          int> = 0>
+-    inline T &emplace(std::initializer_list<Up> il, Args &&... args) {
++    inline constexpr T &emplace(std::initializer_list<Up> il, Args &&... args) {
+       return impl_.template emplace<I>(il, lib::forward<Args>(args)...);
+     }
+ 
+@@ -1739,7 +1746,7 @@ namespace mpark {
+         typename... Args,
+         std::size_t I = detail::find_index_sfinae<T, Ts...>::value,
+         lib::enable_if_t<std::is_constructible<T, Args...>::value, int> = 0>
+-    inline T &emplace(Args &&... args) {
++    inline constexpr T &emplace(Args &&... args) {
+       return impl_.template emplace<I>(lib::forward<Args>(args)...);
+     }
+ 
+@@ -1752,7 +1759,7 @@ namespace mpark {
+                                                std::initializer_list<Up> &,
+                                                Args...>::value,
+                          int> = 0>
+-    inline T &emplace(std::initializer_list<Up> il, Args &&... args) {
++    inline constexpr T &emplace(std::initializer_list<Up> il, Args &&... args) {
+       return impl_.template emplace<I>(il, lib::forward<Args>(args)...);
+     }
+ 
+@@ -1772,7 +1779,7 @@ namespace mpark {
+                             lib::dependent_type<lib::is_swappable<Ts>,
+                                                 Dummy>::value)...>::value,
+                   int> = 0>
+-    inline void swap(variant &that) noexcept(
++    inline constexpr void swap(variant &that) noexcept(
+         lib::all<(std::is_nothrow_move_constructible<Ts>::value &&
+                   lib::is_nothrow_swappable<Ts>::value)...>::value) {
+       impl_.swap(that.impl_);
+@@ -2093,7 +2100,7 @@ namespace mpark {
+ #endif
+ 
+   template <typename... Ts>
+-  inline auto swap(variant<Ts...> &lhs,
++  inline constexpr auto swap(variant<Ts...> &lhs,
+                    variant<Ts...> &rhs) noexcept(noexcept(lhs.swap(rhs)))
+       -> decltype(lhs.swap(rhs)) {
+     lhs.swap(rhs);