Give the SSHIFT volume shift to the saturations

flowexperimental/comp/wells/CompWellFlash.hpp

@@ -85,7 +85,11 @@ void flashWellboreFluidState(CompositionalFluidState<T, FluidSystem>& fluid_stat
             fluid_state.setLvalue(L);
         }
     }
-    T So = max((L * Z_L / (L * Z_L + (1 - L) * Z_V)), 0.0);
+    // Use the same translated molar volumes for saturations and densities.
+    // The phase-label check above uses the unshifted EOS compressibility factors.
+    const auto Vm_L = param_cache.correctedMolarVolume(FluidSystem::oilPhaseIdx);
+    const auto Vm_V = param_cache.correctedMolarVolume(FluidSystem::gasPhaseIdx);
+    T So = max((L * Vm_L / (L * Vm_L + (1 - L) * Vm_V)), 0.0);
     T Sg = max(1 - So, 0.0);
     T sumS = So + Sg;
     So /= sumS;

opm/models/ptflash/flashintensivequantities.hh

@@ -221,6 +221,8 @@ public:
         paramCache.updatePhase(fluidState_, FluidSystem::oilPhaseIdx);
 
         const Scalar R = Opm::Constants<Scalar>::R;
+        // Store compressibility factors from the unshifted EOS roots.
+        // Saturations below use the translated molar volumes, as density does.
         const Evaluation Z_L = (paramCache.molarVolume(FluidSystem::oilPhaseIdx) *
                                 fluidState_.pressure(FluidSystem::oilPhaseIdx)) /
                                (R * fluidState_.temperature(FluidSystem::oilPhaseIdx));
@@ -235,7 +237,9 @@ public:
             Sw = priVars.makeEvaluation(water0Idx, timeIdx);
         }
         const Evaluation L = fluidState_.L();
-        Evaluation So = max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
+        const Evaluation Vm_L = paramCache.correctedMolarVolume(FluidSystem::oilPhaseIdx);
+        const Evaluation Vm_V = paramCache.correctedMolarVolume(FluidSystem::gasPhaseIdx);
+        Evaluation So = max((1 - Sw) * (L * Vm_L / ( L * Vm_L + (1 - L) * Vm_V)), 0.0);
         Evaluation Sg = max(1 - So - Sw, 0.0);
         const Scalar sumS = getValue(So) + getValue(Sg) + getValue(Sw);
         So /= sumS;
@@ -255,8 +259,8 @@ public:
         if (flashVerbosity >= 5) {
              std::cout << "So = " << So << std::endl;
              std::cout << "Sg = " << Sg << std::endl;
-             std::cout << "Z_L = " << Z_L << std::endl;
-             std::cout << "Z_V = " << Z_V << std::endl;
+             std::cout << "Vm_L = " << Vm_L << std::endl;
+             std::cout << "Vm_V = " << Vm_V << std::endl;
          }
 
         /////////////

regressionTests.cmake

@@ -59,6 +59,25 @@ add_test_compareECLFiles(
     compositional
 )
 
+add_test_compareECLFiles(
+  CASENAME
+    sshift_compositional
+  FILENAME
+    SIMPLE_COMP_SSHIFT
+  SIMULATOR
+    flow_comp
+  DEV_SIMULATOR
+    flow_comp3_2p
+  REFERENCE_SIMULATOR
+    flow_comp
+  ABS_TOL
+    ${abs_tol}
+  REL_TOL
+    ${rel_tol}
+  DIR
+    compositional
+)
+
 add_test_compareECLFiles(
   CASENAME
     spe12

tests/test_compwell_jacobian.cpp

@@ -68,10 +68,14 @@ using FluidSystem = Opm::GenericOilGasWaterFluidSystem<Scalar, 3, /*enableWater=
 
 constexpr int numComponents = FluidSystem::numComponents; // 3
 
+// Dimensionless SSHIFT coefficients for CO2, C1 and C10.
+constexpr std::array<Scalar, 3> noVolumeShift{0.0, 0.0, 0.0};
+constexpr std::array<Scalar, 3> volumeShift{-0.0817, -0.1540, 0.0855};
+
 // Register the fixed CO2/Methane/Decane composition with the generic fluid
 // system. The component data is shared static state, so this must run before any
 // flash.
-void registerFluidSystemComponents()
+void registerFluidSystemComponents(const std::array<Scalar, 3>& sshift = noVolumeShift)
 {
     using CO2 = Opm::SimpleCO2<Scalar>;
     using C1  = Opm::C1<Scalar>;
@@ -80,11 +84,14 @@ void registerFluidSystemComponents()
 
     FluidSystem::init();
     FluidSystem::addComponent(CompParam{CO2::name(), CO2::molarMass(), CO2::criticalTemperature(),
-                                        CO2::criticalPressure(), CO2::criticalVolume(), CO2::acentricFactor()});
+                                        CO2::criticalPressure(), CO2::criticalVolume(), CO2::acentricFactor(),
+                                        sshift[0]});
     FluidSystem::addComponent(CompParam{C1::name(), C1::molarMass(), C1::criticalTemperature(),
-                                        C1::criticalPressure(), C1::criticalVolume(), C1::acentricFactor()});
+                                        C1::criticalPressure(), C1::criticalVolume(), C1::acentricFactor(),
+                                        sshift[1]});
     FluidSystem::addComponent(CompParam{C10::name(), C10::molarMass(), C10::criticalTemperature(),
-                                        C10::criticalPressure(), C10::criticalVolume(), C10::acentricFactor()});
+                                        C10::criticalPressure(), C10::criticalVolume(), C10::acentricFactor(),
+                                        sshift[2]});
 }
 
 // The wellbore primary variables that the component masses depend on are the
@@ -172,23 +179,29 @@ computeWellboreQuantities(const T& pressure,
     return q;
 }
 
-} // anonymous namespace
+// A state inside the two-phase region, away from phase transitions that would
+// make the central-difference comparison unreliable.
+constexpr Scalar temperature = 300.0;            // K
+constexpr Scalar wellbore_volume = 21.6e-3;      // m^3 (matches CompWell)
+constexpr Scalar p0 = 10.0e5;                    // Pa
+constexpr Scalar z0_0 = 0.5;
+constexpr Scalar z1_0 = 0.3;
 
-BOOST_AUTO_TEST_CASE(WellboreFlashDerivatives)
-{
-    registerFluidSystemComponents();
+// Resolve the flash more accurately than the finite-difference perturbations.
+constexpr Scalar flash_tolerance = 1.e-8;
 
-    // A composition/pressure that sits comfortably inside the two-phase region,
-    // so the flash is smooth and the central differences are meaningful.
-    const Scalar temperature = 300.0;            // K
-    const Scalar wellbore_volume = 21.6e-3;      // m^3 (matches CompWell)
-    const Scalar p0 = 10.0e5;                    // Pa
-    const Scalar z0_0 = 0.5;
-    const Scalar z1_0 = 0.3;
+// Wellbore mixture density at the test state using the registered components.
+Scalar baseWellboreDensity()
+{
+    const std::array<Scalar, numComponents> z{z0_0, z1_0, 1.0 - z0_0 - z1_0};
+    return computeWellboreQuantities<Scalar>(p0, z, temperature,
+                                             wellbore_volume, flash_tolerance).fluid_density;
+}
 
-    // Tighten the flash tolerance well below the perturbation-induced signal so
-    // the finite differences are not swamped by flash-convergence noise.
-    const Scalar flash_tolerance = 1.e-8;
+// Compare wellbore AD derivatives with central differences for the current
+// component configuration.
+void checkWellboreFlashDerivatives()
+{
 
     // --- Analytical (AD) quantities at the base point ---------------------
     Evaluation P  = Evaluation::createVariable(p0, pIdx);
@@ -271,3 +284,30 @@ BOOST_AUTO_TEST_CASE(WellboreFlashDerivatives)
                    "d(fluid_density)/dx[" + std::to_string(s) + "]");
     }
 }
+
+} // anonymous namespace
+
+BOOST_AUTO_TEST_CASE(WellboreFlashDerivatives)
+{
+    registerFluidSystemComponents();
+    checkWellboreFlashDerivatives();
+}
+
+BOOST_AUTO_TEST_CASE(WellboreFlashDerivativesWithVolumeShift)
+{
+    // Exercise the derivatives with SSHIFT applied to densities and saturations.
+    registerFluidSystemComponents();
+    const Scalar rho_unshifted = baseWellboreDensity();
+
+    registerFluidSystemComponents(volumeShift);
+    const Scalar rho_shifted = baseWellboreDensity();
+
+    // Require a measurable density change so ignoring SSHIFT cannot pass the
+    // derivative check. At this state, the change is about 0.0146 kg/m3;
+    // the 0.001 kg/m3 threshold is well below that change.
+    BOOST_TEST_MESSAGE("wellbore density unshifted = " << rho_unshifted
+                       << ", shifted = " << rho_shifted);
+    BOOST_REQUIRE_GT(std::abs(rho_shifted - rho_unshifted), 1.0e-3);
+
+    checkWellboreFlashDerivatives();
+}