Fix bulk modulus calculation in spiner EOS

CHANGELOG.md

@@ -12,6 +12,7 @@
 - [[PR#357]](https://github.com/lanl/singularity-eos/pull/357) Added support for C++17 (e.g., needed when using newer Kokkos).
 
 ### Fixed (Repair bugs, etc)
+- [[PR370]](https://github.com/lanl/singularity-eos/pull/370) Fix bulk modulus calculation in spiner EOS
 - [[PR343]](https://github.com/lanl/singularity-eos/pull/343) Add chemical potentials to stellar collapse gold files
 - [[PR342]](https://github.com/lanl/singularity-eos/pull/342) Fix missing using statement in stellar collapse root finding routines
 - [[PR341]](https://github.com/lanl/singularity-eos/pull/341) Short-circuit HDF5 machinery when cray-wrappers used in-tree
@@ -36,7 +37,7 @@ Date: 11/28/2023
 - [[PR278]](https://github.com/lanl/singularity-eos/pull/278) Fixed EOSPAC unit conversion errors for scalar lookups
 - [[PR316]](https://github.com/lanl/singularity-eos/pull/316) removed `fmax-errors=3` from `singularity-eos` compile flags
 - [[PR296]](https://github.com/lanl/singularity-eos/pull/296) changed `CMAKE_SOURCE_DIR` to `PROJECT_SOURCE_DIR` to fix downstream submodule build
-- [[PR291]](https://github.com/lanl/singularity-eos/pull/291) package.py updates to reflect new CMake options 
+- [[PR291]](https://github.com/lanl/singularity-eos/pull/291) package.py updates to reflect new CMake options
 - [[PR290]](https://github.com/lanl/singularity-eos/pull/290) Added target guards on export config
 - [[PR288]](https://github.com/lanl/singularity-eos/pull/288) Don't build tests that depend on spiner when spiner is disabled
 - [[PR287]](https://github.com/lanl/singularity-eos/pull/287) Fix testing logic with new HDF5 options

singularity-eos/eos/eos_spiner.hpp

@@ -870,15 +870,15 @@ inline void SpinerEOSDependsRhoT::fixBulkModulus_() {
     for (int i = 0; i < numT_; i++) {
       Real lT = bMod_.range(0).x(i);
       Real press = P_.interpToReal(lRho, lT);
-      Real DPDR = dPdRho_.interpToReal(lRho, lT);
-      Real DPDE = dPdE_.interpToReal(lRho, lT);
-      Real DEDR = dEdRho_.interpToReal(lRho, lT);
-      Real DTDE = dTdE_.interpToReal(lRho, lT);
+      Real DPDR_E = dPdRho_.interpToReal(lRho, lT);
+      Real DPDE_R = dPdE_.interpToReal(lRho, lT);
+      Real DEDR_T = dEdRho_.interpToReal(lRho, lT);
+      Real DPDR_T = DPDR_E + DPDE_R * DEDR_T;
       Real bMod;
-      if (DPDE > 0.0 && rho > 0.0) {
-        bMod = rho * DPDR + DPDE * (press / rho - rho * DEDR);
+      if (DPDE_R > 0.0 && rho > 0.0) {
+        bMod = rho * DPDR_E + DPDE_R * (press / rho);
       } else if (rho > 0.0) {
-        bMod = std::max(rho * DPDR, 0.0);
+        bMod = std::max(rho * DPDR_T, 0.0);
       } else {
         bMod = 0.0;
       }
@@ -1637,15 +1637,15 @@ inline void SpinerEOSDependsRhoSie::calcBMod_(SP5Tables &tables) {
     Real rho = fromLog_(lRho, lRhoOffset_);
     for (int i = 0; i < tables.bMod.dim(1); i++) {
       Real press = tables.P(j, i);
-      Real DPDR = tables.dPdRho(j, i);
-      Real DPDE = tables.dPdE(j, i);
-      Real DEDR = tables.dEdRho(j, i);
-      Real DTDE = tables.dTdE(j, i);
+      Real DPDR_E = tables.dPdRho(j, i);
+      Real DPDE_R = tables.dPdE(j, i);
+      Real DEDR_T = tables.dEdRho(j, i);
+      Real DPDR_T = DPDR_E + DPDE_R * DEDR_T;
       Real bMod;
-      if (DPDE > 0.0 && rho > 0.0) {
-        bMod = rho * DPDR + DPDE * (press / rho - rho * DEDR);
+      if (DPDE_R > 0.0 && rho > 0.0) {
+        bMod = rho * DPDR_E + DPDE_R * (press / rho);
       } else if (rho > 0.0) {
-        bMod = std::max(rho * DPDR, 0.0);
+        bMod = std::max(rho * DPDR_T, 0.0);
       } else {
         bMod = 0.0;
       }