From 2f14f7b062da5a392500f698a194b39f821bd0a2 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 8 Sep 2026 14:10:32 +0200 Subject: [PATCH 1/5] Apply SSHIFT to compositional density and viscosity Pass the parsed component SSHIFT coefficients to the fluid system and cache each phase's volume translation, sum_c x_c s_c b_c. Use the corrected molar volume for density and the molar-density input to the LBC correlation. Reject corrected volumes that are non-positive or NaN. Keep the cached EOS molar volume unshifted for fugacity calculations. At equal phase pressure and temperature, the translation factors cancel from the equilibrium ratios, preserving the phase split. Add tests for density, viscosity, equilibrium ratios, AD derivatives, parsed SSHIFT values and invalid volumes. Saturation updates are handled in opm-simulators PR #7391. --- CMakeLists_files.cmake | 1 + .../GenericOilGasWaterFluidSystem.hpp | 42 +- .../fluidsystems/PTFlashParameterCache.hpp | 64 +++ .../viscositymodels/ViscosityModels.hpp | 19 +- tests/material/test_volume_shift.cpp | 456 ++++++++++++++++++ 5 files changed, 572 insertions(+), 10 deletions(-) create mode 100644 tests/material/test_volume_shift.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 5b684b52fb9..44980e8e074 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -717,6 +717,7 @@ list(APPEND DUNE_TEST_SOURCE_FILES tests/material/test_ptflash_ssi_newton_fallback.cpp tests/material/test_tabulation.cpp tests/material/test_threecomponents_ptflash.cpp + tests/material/test_volume_shift.cpp ) if(dune-common_FOUND) diff --git a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp index 7d01b3f9a51..1de1de1683e 100644 --- a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp +++ b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp @@ -93,15 +93,18 @@ namespace Opm { Scalar critic_pres; // unit: parscal Scalar critic_vol; // unit: m^3/kmol Scalar acentric_factor; // unit: dimension less + Scalar volume_shift; // dimensionless SSHIFT coefficient ComponentParam(const std::string_view name_, const Scalar molar_mass_, const Scalar critic_temp_, - const Scalar critic_pres_, const Scalar critic_vol_, const Scalar acentric_factor_) + const Scalar critic_pres_, const Scalar critic_vol_, const Scalar acentric_factor_, + const Scalar volume_shift_ = 0.0) : name(name_), molar_mass(molar_mass_), critic_temp(critic_temp_), critic_pres(critic_pres_), critic_vol(critic_vol_), - acentric_factor(acentric_factor_) + acentric_factor(acentric_factor_), + volume_shift(volume_shift_) {} }; @@ -177,7 +180,10 @@ namespace Opm { static_cast(eos_props.critical_temperature[c]), static_cast(eos_props.critical_pressure[c]), static_cast(eos_props.critical_volume[c] * 1.e3), - static_cast(eos_props.acentric_factors[c])}); + static_cast(eos_props.acentric_factors[c]), + c < eos_props.volume_shifts.size() + ? static_cast(eos_props.volume_shifts[c]) + : Scalar{0}}); } const auto& bic = eos_props.binary_interaction_coefficient; @@ -220,6 +226,23 @@ namespace Opm { return component_param_[compIdx].acentric_factor; } + + /*! + * \brief Dimensionless volume-translation coefficient (SSHIFT). + * + * The component volume correction is s_c b_c, where b_c is the + * equation-of-state covolume. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar volumeShift(unsigned compIdx) + { + assert(isConsistent()); + assert(compIdx < numComponents); + + return component_param_[compIdx].volume_shift; + } + /*! * \brief Critical temperature of a component [K]. * @@ -325,7 +348,10 @@ namespace Opm { assert(phaseIdx < numPhases); if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) { - return decay(fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx)); + // Density uses the translated volume; fugacity coefficients + // use the unshifted equation-of-state root. + const auto Vm = paramCache.correctedMolarVolume(phaseIdx); + return decay(fluidState.averageMolarMass(phaseIdx) / Vm); } else { const LhsEval& p = decay(fluidState.pressure(phaseIdx)); @@ -348,8 +374,12 @@ namespace Opm { assert(phaseIdx < numPhases); if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) { - // Use LBC method to calculate viscosity - return decay(ViscosityModel::LBC(fluidState, paramCache, phaseIdx)); + // LBC is a reduced-density correlation, so use the physical + // molar density after applying SSHIFT. + const auto rho = density(fluidState, paramCache, phaseIdx); + const auto molarDensity = rho / fluidState.averageMolarMass(phaseIdx); + return decay( + ViscosityModel::LBCWithMolarDensity(fluidState, molarDensity, phaseIdx)); } else { const LhsEval& p = decay(fluidState.pressure(phaseIdx)); diff --git a/opm/material/fluidsystems/PTFlashParameterCache.hpp b/opm/material/fluidsystems/PTFlashParameterCache.hpp index 85b2fb2a524..9c700713a64 100644 --- a/opm/material/fluidsystems/PTFlashParameterCache.hpp +++ b/opm/material/fluidsystems/PTFlashParameterCache.hpp @@ -30,6 +30,8 @@ #ifndef OPM_PTFlash_PARAMETER_CACHE_HPP #define OPM_PTFlash_PARAMETER_CACHE_HPP +#include + #include #include #include @@ -37,6 +39,8 @@ #include +#include + #include namespace Opm { @@ -80,8 +84,10 @@ class PTFlashParameterCache { VmUpToDate_[oilPhaseIdx] = false; Valgrind::SetUndefined(Vm_[oilPhaseIdx]); + Valgrind::SetUndefined(volumeShift_[oilPhaseIdx]); VmUpToDate_[gasPhaseIdx] = false; Valgrind::SetUndefined(Vm_[gasPhaseIdx]); + Valgrind::SetUndefined(volumeShift_[gasPhaseIdx]); oilPhaseParams_.setEOSType(eos_type); gasPhaseParams_.setEOSType(eos_type); @@ -286,6 +292,44 @@ class PTFlashParameterCache return Vm_[phaseIdx]; } + /*! + * \brief Phase volume translation, sum_c x_c s_c b_c [m^3/mol]. + * + * Call updatePhase() before reading this cached value. The translation + * is stored separately from the equation-of-state root in molarVolume(). + * correctedMolarVolume() subtracts it to obtain the physical molar volume. + * + * \param phaseIdx The fluid phase of interest + */ + Scalar phaseVolumeShift(unsigned phaseIdx) const + { + return volumeShift_[phaseIdx]; + } + + /*! + * \brief Molar volume after applying SSHIFT [m^3/mol]. + * + * Density, viscosity and phase saturations use this physical volume. + * Fugacity coefficients and equation-of-state compressibility factors + * use the unshifted molarVolume(). At equal phase pressure and temperature, + * the volume-translation factors cancel from the equilibrium ratios. + * + * \param phaseIdx The fluid phase of interest + */ + Scalar correctedMolarVolume(unsigned phaseIdx) const + { + const Scalar Vm = molarVolume(phaseIdx) - phaseVolumeShift(phaseIdx); + + // Reject non-positive or NaN volumes before computing fluid properties. + if (!(scalarValue(Vm) > 0)) { + throw NumericalProblem( + fmt::format("The SSHIFT volume shift of phase {} leaves a corrected " + "molar volume of {}, which is not positive.", + phaseIdx, scalarValue(Vm))); + } + return Vm; + } + /*! * \brief Returns the Peng-Robinson mixture parameters for the oil @@ -381,6 +425,7 @@ class PTFlashParameterCache unsigned phaseIdx) { VmUpToDate_[phaseIdx] = true; + volumeShift_[phaseIdx] = computeVolumeShift_(fluidState, phaseIdx); // calculate molar volume of the phase (we will need this for the // fugacity coefficients and the density anyway) @@ -415,8 +460,27 @@ class PTFlashParameterCache }; } + //! \brief Compute the phase volume translation, sum_c x_c s_c b_c [m^3/mol]. + template + Scalar computeVolumeShift_(const FluidState& fluidState, unsigned phaseIdx) const + { + // b_c from the dimensionless B_c = b_c p / (R T). + const Scalar T = decay(fluidState.temperature(phaseIdx)); + const Scalar p = decay(fluidState.pressure(phaseIdx)); + const Scalar RT_p = Constants::R * T / p; + + Scalar shift = 0; + for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { + const Scalar b = decay(Bi(phaseIdx, compIdx)) * RT_p; + shift += decay(fluidState.moleFraction(phaseIdx, compIdx)) + * FluidSystem::volumeShift(compIdx) * b; + } + return shift; + } + bool VmUpToDate_[numPhases]; Scalar Vm_[numPhases]; + Scalar volumeShift_[numPhases]; OilPhaseParams oilPhaseParams_; GasPhaseParams gasPhaseParams_; diff --git a/opm/material/viscositymodels/ViscosityModels.hpp b/opm/material/viscositymodels/ViscosityModels.hpp index 83a8793c953..3121cbf9a63 100644 --- a/opm/material/viscositymodels/ViscosityModels.hpp +++ b/opm/material/viscositymodels/ViscosityModels.hpp @@ -62,11 +62,24 @@ class ViscosityModels const Params& /*paramCache*/, unsigned phaseIdx) { - const Scalar MPa_atm = 0.101325; const Scalar R = Opm::Constants::R; const auto& T = Opm::decay(fluidState.temperature(phaseIdx)); const auto& P = Opm::decay(fluidState.pressure(phaseIdx)); const auto& Z = Opm::decay(fluidState.compressFactor(phaseIdx)); + const LhsEval molarDensity = P / (R * T * Z); + + return LBCWithMolarDensity(fluidState, molarDensity, phaseIdx); + } + + // LBC correlation at the supplied physical molar density [mol/m^3]. + template + static LhsEval LBCWithMolarDensity(const FluidState& fluidState, + const MolarDensity& molarDensity, + unsigned phaseIdx) + { + const Scalar MPa_atm = 0.101325; + const auto& T = Opm::decay(fluidState.temperature(phaseIdx)); LhsEval sumVolume = 0.0; for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { @@ -76,9 +89,7 @@ class ViscosityModels } LhsEval rho_pc = 1.0 / sumVolume; - LhsEval V = (R * T * Z)/P; - LhsEval rho = 1.0 / V; - LhsEval rho_r = rho / rho_pc; + const LhsEval rho_r = Opm::decay(molarDensity) / rho_pc; LhsEval xsum_T_c = 0.0; // mixture pseudocritical temperature LhsEval xsum_Mm = 0.0; // mixture molar mass diff --git a/tests/material/test_volume_shift.cpp b/tests/material/test_volume_shift.cpp new file mode 100644 index 00000000000..94018769a48 --- /dev/null +++ b/tests/material/test_volume_shift.cpp @@ -0,0 +1,456 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2026 SINTEF Digital + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief Tests SSHIFT density, viscosity, equilibrium ratios and derivatives. + * + * Component properties come from COMP_EQUIL_1D_VERTICAL_EQLNUM7. Its + * reference gas density is 165.87 kg/m3 at 272.599 bar and 393.15 K; + * the unshifted calculation gives 172.84 kg/m3. The viscosity reference + * comes from VGAS at report step 1 of the same run. + */ +#include "config.h" + +#define BOOST_TEST_MODULE VolumeShift +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +namespace { + +using Scalar = double; +constexpr int numComponents = 7; + +using FluidSystem = Opm::GenericOilGasWaterFluidSystem; +// Changing enableWater gives the reference its own static component data. +// It uses the same component properties with zero shifts; only oil and gas +// properties are evaluated in these comparisons. +using UnshiftedSystem = Opm::GenericOilGasWaterFluidSystem; +using CompVec = std::array; + +constexpr auto eosType = Opm::CompositionalConfig::EOSType::PR; + +// Reference conditions at the top of the column: RTEMP = 120 degC. +constexpr Scalar temperature = 393.15; +constexpr Scalar pressure = 272.599e5; + +// A seven-component reservoir fluid: methane, CO2, ethane, propane and three +// lumped fractions. The shifts are the SSHIFT values of the same deck. +const CompVec z{0.88583, 0.0437, 0.034, 0.0189, 0.0152, 0.0023, 0.00007}; +const CompVec shift{-0.1595, -0.0817, -0.1134, -0.0863, -0.0243568, + 0.10784125400986, 0.206892761354429}; + +struct Component +{ + const char* name; + Scalar molarMass; // kg/mol + Scalar criticalT; // K + Scalar criticalP; // Pa + Scalar criticalV; // m^3/kmol + Scalar acentric; +}; + +// Deck component properties converted to the units listed above. +const std::array components{{ + {"CH4", 0.016043, 190.60, 45.40e5, 0.09900000, 0.00800}, + {"CO2", 0.044010, 304.20, 72.80e5, 0.09400000, 0.22500}, + {"C2H6", 0.030070, 305.40, 48.20e5, 0.14800000, 0.09800}, + {"C3H8", 0.044097, 369.80, 41.90e5, 0.20300000, 0.15200}, + {"C4-C9", 0.076900, 474.70, 32.40e5, 0.32321204, 0.25490}, + {"C10-C20", 0.163000, 646.00, 20.42e5, 0.63255728, 0.55100}, + {"C21+", 0.310850, 812.00, 14.50e5, 0.97988801, 0.90900}, +}}; + +// Initialize both systems once for the tests that share their static data. +struct Fixture +{ + Fixture() + { + using CompParam = typename FluidSystem::ComponentParam; + FluidSystem::init(); + for (int c = 0; c < numComponents; ++c) { + const auto& p = components[c]; + FluidSystem::addComponent(CompParam{p.name, p.molarMass, p.criticalT, + p.criticalP, p.criticalV, p.acentric, + shift[c]}); + } + + using UnshiftedParam = typename UnshiftedSystem::ComponentParam; + UnshiftedSystem::init(); + for (int c = 0; c < numComponents; ++c) { + const auto& p = components[c]; + UnshiftedSystem::addComponent(UnshiftedParam{p.name, p.molarMass, p.criticalT, + p.criticalP, p.criticalV, + p.acentric, 0.0}); + } + } +}; + +/// Gas density, viscosity and fugacity coefficients at the reference conditions. +struct PhaseState +{ + Scalar density{}; + Scalar viscosity{}; + CompVec fugacityCoefficient{}; +}; + +PhaseState gasState() +{ + Opm::CompositionalFluidState fs; + fs.setTemperature(temperature); + fs.setPressure(FluidSystem::oilPhaseIdx, pressure); + fs.setPressure(FluidSystem::gasPhaseIdx, pressure); + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, z[c]); + fs.setMoleFraction(FluidSystem::oilPhaseIdx, c, z[c]); + } + + typename FluidSystem::template ParameterCache paramCache(eosType); + paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx); + + // Store the unshifted EOS compressibility factor, as the simulator does. + const Scalar Z = paramCache.molarVolume(FluidSystem::gasPhaseIdx) * pressure + / (Opm::Constants::R * temperature); + fs.setCompressFactor(FluidSystem::gasPhaseIdx, Z); + + PhaseState state; + state.density = FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx); + state.viscosity = FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx); + for (int c = 0; c < numComponents; ++c) { + state.fugacityCoefficient[c] = + FluidSystem::fugacityCoefficient(fs, paramCache, FluidSystem::gasPhaseIdx, c); + } + return state; +} + +// Evaluate the same property calculation with scalar or AD inputs. +template +std::pair shiftedGasProperties(const Eval& p, const Eval& xCH4) +{ + Opm::CompositionalFluidState fs; + fs.setTemperature(Eval{temperature}); + fs.setPressure(FluidSystem::oilPhaseIdx, p); + fs.setPressure(FluidSystem::gasPhaseIdx, p); + fs.setMoleFraction(FluidSystem::gasPhaseIdx, 0, xCH4); + for (int c = 1; c < numComponents; ++c) { + fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, Eval{z[c]}); + } + + typename FluidSystem::template ParameterCache paramCache(eosType); + paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx); + + return {FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx), + FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx)}; +} + +} // Anonymous namespace + +BOOST_GLOBAL_FIXTURE(Fixture); + +BOOST_AUTO_TEST_CASE(ShiftMovesTheDensityOntoTheReference) +{ + const auto gas = gasState(); + + // COMPVD supplies the test composition; the reference uses the equilibrated + // composition at the top of the column. The 1% tolerance allows for this + // difference while excluding the unshifted result of 172.84 kg/m3. + BOOST_CHECK_CLOSE(gas.density, 165.87, 1.0); +} + +BOOST_AUTO_TEST_CASE(ShiftLeavesTheEquilibriumRatiosAlone) +{ + // At equal phase pressure and temperature, the Peneloux factor + // exp(-s_c b_c p / (R T)) cancels from phi_c^liquid / phi_c^vapour. + // Compare this ratio with a separate system whose shifts are zero. + Opm::CompositionalFluidState fs; + Opm::CompositionalFluidState fsRef; + fs.setTemperature(temperature); + fsRef.setTemperature(temperature); + for (int ph : {FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}) { + fs.setPressure(ph, pressure); + fsRef.setPressure(ph, pressure); + } + // Use distinct liquid and vapour compositions so the ratio does not + // trivially compare identical phase properties. These are prescribed + // compositions; this test does not perform a flash calculation. + const CompVec x{0.55, 0.09, 0.12, 0.10, 0.11, 0.028, 0.002}; + const CompVec y{0.95, 0.03, 0.012, 0.005, 0.0025, 0.0004, 0.0001}; + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(FluidSystem::oilPhaseIdx, c, x[c]); + fsRef.setMoleFraction(UnshiftedSystem::oilPhaseIdx, c, x[c]); + fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, y[c]); + fsRef.setMoleFraction(UnshiftedSystem::gasPhaseIdx, c, y[c]); + } + + typename FluidSystem::template ParameterCache pc(eosType); + typename UnshiftedSystem::template ParameterCache pcRef(eosType); + for (int ph : {FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}) { + pc.updatePhase(fs, ph); + pcRef.updatePhase(fsRef, ph); + } + + // Confirm that the two systems have different volume translations. + BOOST_CHECK_GT(std::abs(pc.phaseVolumeShift(FluidSystem::gasPhaseIdx)), 0.0); + BOOST_CHECK_SMALL(pcRef.phaseVolumeShift(UnshiftedSystem::gasPhaseIdx), 1.0e-30); + + for (int c = 0; c < numComponents; ++c) { + const Scalar K = + FluidSystem::fugacityCoefficient(fs, pc, FluidSystem::oilPhaseIdx, c) / + FluidSystem::fugacityCoefficient(fs, pc, FluidSystem::gasPhaseIdx, c); + const Scalar KRef = + UnshiftedSystem::fugacityCoefficient(fsRef, pcRef, UnshiftedSystem::oilPhaseIdx, c) / + UnshiftedSystem::fugacityCoefficient(fsRef, pcRef, UnshiftedSystem::gasPhaseIdx, c); + BOOST_CHECK_CLOSE(K, KRef, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(ShiftDoesNotReachTheCachedVolume) +{ + // Check that the corrected volume subtracts the cached translation and + // that density uses the resulting volume. + Opm::CompositionalFluidState fs; + fs.setTemperature(temperature); + fs.setPressure(FluidSystem::oilPhaseIdx, pressure); + fs.setPressure(FluidSystem::gasPhaseIdx, pressure); + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, z[c]); + fs.setMoleFraction(FluidSystem::oilPhaseIdx, c, z[c]); + } + typename FluidSystem::template ParameterCache pc(eosType); + pc.updatePhase(fs, FluidSystem::gasPhaseIdx); + + const Scalar Vm = pc.molarVolume(FluidSystem::gasPhaseIdx); + const Scalar shift = pc.phaseVolumeShift(FluidSystem::gasPhaseIdx); + BOOST_CHECK_CLOSE(pc.correctedMolarVolume(FluidSystem::gasPhaseIdx), + Vm - shift, 1.0e-10); + // Density uses the translated molar volume. + BOOST_CHECK_CLOSE(FluidSystem::density(fs, pc, FluidSystem::gasPhaseIdx), + fs.averageMolarMass(FluidSystem::gasPhaseIdx) / (Vm - shift), + 1.0e-10); +} + +BOOST_AUTO_TEST_CASE(ShiftMovesTheLiquidDensityToo) +{ + // Check the oil-phase path against an independently calculated translation + // using the component SSHIFT values and b_c = Omega_b R Tc/pc. + Opm::CompositionalFluidState fs; + fs.setTemperature(temperature); + fs.setPressure(FluidSystem::oilPhaseIdx, pressure); + fs.setPressure(FluidSystem::gasPhaseIdx, pressure); + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(FluidSystem::oilPhaseIdx, c, z[c]); + fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, z[c]); + } + typename FluidSystem::template ParameterCache pc(eosType); + pc.updatePhase(fs, FluidSystem::oilPhaseIdx); + + constexpr Scalar OmegaB = 0.0777960739; // Peng-Robinson + constexpr Scalar R = 8.31446261815324; + Scalar expectedShift = 0.0; + for (int c = 0; c < numComponents; ++c) { + const auto& p = components[c]; + const Scalar b = OmegaB * R * p.criticalT / p.criticalP; + expectedShift += z[c] * shift[c] * b; + } + // Allow for rounding in the independently specified EOS constants. + BOOST_CHECK_CLOSE(pc.phaseVolumeShift(FluidSystem::oilPhaseIdx), expectedShift, 1.0e-3); + + const Scalar Vm = pc.molarVolume(FluidSystem::oilPhaseIdx); + const Scalar rho = FluidSystem::density(fs, pc, FluidSystem::oilPhaseIdx); + BOOST_CHECK_CLOSE(rho, + fs.averageMolarMass(FluidSystem::oilPhaseIdx) / (Vm - expectedShift), + 1.0e-4); + // The negative mixture translation increases the volume and lowers density. + BOOST_CHECK_LT(rho, fs.averageMolarMass(FluidSystem::oilPhaseIdx) / Vm); +} + +BOOST_AUTO_TEST_CASE(ShiftMovesTheViscosityOntoTheReference) +{ + // Reference VGAS for cell 1 at report step 1 is 0.0223926 cP. The 1% + // tolerance allows for the different pressure (272.64 bar) and equilibrated + // composition at that step. The unshifted result is about 2.3% higher. + BOOST_CHECK_CLOSE(gasState().viscosity, 2.23926e-5, 1.0); +} + +BOOST_AUTO_TEST_CASE(WithoutAShiftTheTwoViscosityPathsAgree) +{ + // With zero shifts, the explicit molar-density and compressibility-factor + // entry points must produce the same LBC viscosity. + Opm::CompositionalFluidState fs; + fs.setTemperature(temperature); + fs.setPressure(UnshiftedSystem::gasPhaseIdx, pressure); + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(UnshiftedSystem::gasPhaseIdx, c, z[c]); + } + + typename UnshiftedSystem::template ParameterCache paramCache(eosType); + paramCache.updatePhase(fs, UnshiftedSystem::gasPhaseIdx); + const Scalar Z = paramCache.molarVolume(UnshiftedSystem::gasPhaseIdx) * pressure + / (Opm::Constants::R * temperature); + fs.setCompressFactor(UnshiftedSystem::gasPhaseIdx, Z); + + const Scalar viaMolarDensity = + UnshiftedSystem::viscosity(fs, paramCache, UnshiftedSystem::gasPhaseIdx); + const Scalar viaCompressFactor = + Opm::ViscosityModels::LBC(fs, paramCache, + UnshiftedSystem::gasPhaseIdx); + + BOOST_CHECK_CLOSE(viaMolarDensity, viaCompressFactor, 1.0e-8); +} + +BOOST_AUTO_TEST_CASE(TheShiftedPathCarriesItsDerivatives) +{ + // Compare pressure and composition derivatives with central differences. + // The composition perturbation also exercises the derivative of the + // mixture volume translation. + using Evaluation = Opm::DenseAd::Evaluation; + constexpr unsigned pIdx = 0; + constexpr unsigned xIdx = 1; + + const auto ad = shiftedGasProperties( + Evaluation::createVariable(pressure, pIdx), + Evaluation::createVariable(z[0], xIdx)); + + const auto plain = shiftedGasProperties(pressure, z[0]); + BOOST_CHECK_CLOSE(ad.first.value(), plain.first, 1.0e-8); + BOOST_CHECK_CLOSE(ad.second.value(), plain.second, 1.0e-8); + + const Scalar dp = 1.0e-6 * pressure; + const auto pUp = shiftedGasProperties(pressure + dp, z[0]); + const auto pDown = shiftedGasProperties(pressure - dp, z[0]); + BOOST_CHECK_CLOSE(ad.first.derivative(pIdx), + (pUp.first - pDown.first) / (2 * dp), 1.0e-4); + BOOST_CHECK_CLOSE(ad.second.derivative(pIdx), + (pUp.second - pDown.second) / (2 * dp), 1.0e-4); + + const Scalar dx = 1.0e-6; + const auto xUp = shiftedGasProperties(pressure, z[0] + dx); + const auto xDown = shiftedGasProperties(pressure, z[0] - dx); + BOOST_CHECK_CLOSE(ad.first.derivative(xIdx), + (xUp.first - xDown.first) / (2 * dx), 1.0e-4); + BOOST_CHECK_CLOSE(ad.second.derivative(xIdx), + (xUp.second - xDown.second) / (2 * dx), 1.0e-4); +} + +// Separate static component data for the deck-initialization test. +using DeckSystem = Opm::GenericOilGasWaterFluidSystem; + +BOOST_AUTO_TEST_CASE(ParsedShiftReachesTheFluidSystem) +{ + // Verify the handoff from parsed SSHIFT values to the component parameters + // through initFromState(). The other tests register components directly. + const auto deck = Opm::Parser{}.parseString(R"( +RUNSPEC +METRIC +DIMENS + 1 1 1 / +COMPS +3 / +TABDIMS + 1 / +OIL +GAS +WATER +GRID +DXV + 1 / +DYV + 1 / +DZV + 1 / +DEPTHZ + 4*2000 / +PROPS +CNAMES + C1 C10 CO2 / +TCRIT + 190.6 617.7 304.2 / +PCRIT + 46.0 21.1 73.8 / +VCRIT + 0.0990 0.6240 0.0940 / +MW + 16.043 142.285 44.010 / +ACF + 0.008 0.4885 0.225 / +SSHIFT + -0.1595 0.10784 -0.0817 / +SOLUTION +SCHEDULE +END +)"); + + const auto eclState = Opm::EclipseState{ deck }; + const auto schedule = Opm::Schedule{ deck, eclState }; + + BOOST_REQUIRE_NO_THROW(DeckSystem::initFromState(eclState, schedule)); + + const std::array expected{-0.1595, 0.10784, -0.0817}; + for (unsigned c = 0; c < 3; ++c) { + BOOST_CHECK_CLOSE(DeckSystem::volumeShift(c), expected[c], 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(AnOverlargeShiftIsRejectedRatherThanReturned) +{ + // A translation larger than the EOS molar volume must raise an error + // before density is evaluated from a non-positive volume. + using BigShift = Opm::GenericOilGasWaterFluidSystem; + using CompParam = typename BigShift::ComponentParam; + BigShift::init(); + BigShift::addComponent(CompParam{"C1", 0.016043, 190.60, 45.40e5, 0.099, 0.008, + /*volume_shift=*/1.0e3}); + + Opm::CompositionalFluidState fs; + fs.setTemperature(Scalar{400}); + fs.setPressure(BigShift::oilPhaseIdx, Scalar{100e5}); + fs.setPressure(BigShift::gasPhaseIdx, Scalar{100e5}); + fs.setMoleFraction(BigShift::gasPhaseIdx, 0, Scalar{1}); + fs.setMoleFraction(BigShift::oilPhaseIdx, 0, Scalar{1}); + + typename BigShift::template ParameterCache pc(eosType); + pc.updatePhase(fs, BigShift::gasPhaseIdx); + + BOOST_CHECK_THROW(BigShift::density(fs, pc, BigShift::gasPhaseIdx), + Opm::NumericalProblem); +} From 05adfce244a7918cb84aae5508855840b65b4781 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 8 Sep 2026 14:01:49 +0200 Subject: [PATCH 2/5] Reset component registrations when reinitializing the fluid system Clear the static component list in init() so a later initialization can register a different configuration. Without this reset, addComponent() rejects new registrations once the list is full, leaving the previous component properties and SSHIFT coefficients active. --- opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp index 1de1de1683e..9332acd6bc8 100644 --- a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp +++ b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp @@ -205,6 +205,9 @@ namespace Opm { static void init() { waterPvt_ = std::make_shared(); + // Discard the previous configuration so subsequent component + // registrations replace it. + component_param_.clear(); component_param_.reserve(numComponents); } From bc077a6c5900663e9f81b73eee1395375751b8b2 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 8 Sep 2026 14:47:26 +0200 Subject: [PATCH 3/5] Preserve PTFlash cache compatibility and fix API test initialization Default to zero translation for fluid systems without volumeShift(). This restores compilation of the CO2/brine flash, cubic-EOS finite-guard and SSI/Newton fallback tests. Let checkFluidSystem() accept an initializer and run it before creating the state and cache. Register the generic fluid system's components for each scalar/AD check, preserving the reset in init(). Validation: six focused opm-common tests and both compositional well tests pass. SIMPLE_COMP_SSHIFT completes with flow_comp3_2p. --- opm/material/checkFluidSystem.hpp | 17 ++++++-- .../fluidsystems/PTFlashParameterCache.hpp | 21 +++++----- tests/material/test_fluidsystems.cpp | 39 ++++++++++--------- 3 files changed, 46 insertions(+), 31 deletions(-) diff --git a/opm/material/checkFluidSystem.hpp b/opm/material/checkFluidSystem.hpp index 2206bd4feaa..061eed143f7 100644 --- a/opm/material/checkFluidSystem.hpp +++ b/opm/material/checkFluidSystem.hpp @@ -291,14 +291,18 @@ ParameterCache initParamCache() /*! * \brief Checks whether a fluid system adheres to the specification. */ -template -void checkFluidSystem() +template +void checkFluidSystem(Initializer initializeFluidSystem) { std::cout << "Testing fluid system '" << Opm::getDemangledType() << ", RhsEval = " << Opm::getDemangledType() << ", LhsEval = " << Opm::getDemangledType() << "'\n"; + // Initialize before creating the fluid state or its cache. Configurable + // systems need to register their components after init() clears them. + initializeFluidSystem(); + // make sure the fluid system provides the number of phases and // the number of components static constexpr int numPhases = FluidSystem::numPhases; @@ -354,7 +358,6 @@ void checkFluidSystem() val = 2*val; // get rid of GCC warning (only occurs with paranoid warning flags) // actually check the fluid system API - try { FluidSystem::init(); } catch (...) {}; for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { fs.restrictToPhase(static_cast(phaseIdx)); fs.allowPressure(FluidSystem::isCompressible(phaseIdx)); @@ -405,4 +408,12 @@ void checkFluidSystem() } } +template +void checkFluidSystem() +{ + checkFluidSystem([] { + try { FluidSystem::init(); } catch (...) {}; + }); +} + #endif diff --git a/opm/material/fluidsystems/PTFlashParameterCache.hpp b/opm/material/fluidsystems/PTFlashParameterCache.hpp index 9c700713a64..8b4a5deefcf 100644 --- a/opm/material/fluidsystems/PTFlashParameterCache.hpp +++ b/opm/material/fluidsystems/PTFlashParameterCache.hpp @@ -464,16 +464,19 @@ class PTFlashParameterCache template Scalar computeVolumeShift_(const FluidState& fluidState, unsigned phaseIdx) const { - // b_c from the dimensionless B_c = b_c p / (R T). - const Scalar T = decay(fluidState.temperature(phaseIdx)); - const Scalar p = decay(fluidState.pressure(phaseIdx)); - const Scalar RT_p = Constants::R * T / p; - Scalar shift = 0; - for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { - const Scalar b = decay(Bi(phaseIdx, compIdx)) * RT_p; - shift += decay(fluidState.moleFraction(phaseIdx, compIdx)) - * FluidSystem::volumeShift(compIdx) * b; + // Fluid systems without SSHIFT support use a zero translation. + if constexpr (requires { FluidSystem::volumeShift(0u); }) { + // b_c from the dimensionless B_c = b_c p / (R T). + const Scalar T = decay(fluidState.temperature(phaseIdx)); + const Scalar p = decay(fluidState.pressure(phaseIdx)); + const Scalar RT_p = Constants::R * T / p; + + for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { + const Scalar b = decay(Bi(phaseIdx, compIdx)) * RT_p; + shift += decay(fluidState.moleFraction(phaseIdx, compIdx)) + * FluidSystem::volumeShift(compIdx) * b; + } } return shift; } diff --git a/tests/material/test_fluidsystems.cpp b/tests/material/test_fluidsystems.cpp index 4743d8bc084..28aea9de2b7 100644 --- a/tests/material/test_fluidsystems.cpp +++ b/tests/material/test_fluidsystems.cpp @@ -393,19 +393,17 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GenericFluidSystem, Scalar, ScalarTypes) using Evaluation = Opm::DenseAd::Evaluation; using FluidSystem = Opm::GenericOilGasWaterFluidSystem; - registerComponent>(); - registerComponent>(); - registerComponent>(); - registerComponent>(); - - // initialize water pvt - using WaterPvt = typename FluidSystem::WaterPvt; - std::shared_ptr waterPvt; - FluidSystem::setWaterPvt(waterPvt); - - checkFluidSystem(); - checkFluidSystem(); - checkFluidSystem(); + const auto initializeFluidSystem = [] { + FluidSystem::init(); + registerComponent>(); + registerComponent>(); + registerComponent>(); + registerComponent>(); + }; + + checkFluidSystem(initializeFluidSystem); + checkFluidSystem(initializeFluidSystem); + checkFluidSystem(initializeFluidSystem); } BOOST_AUTO_TEST_CASE_TEMPLATE(GenericFluidSystemNoWater, Scalar, ScalarTypes) @@ -413,11 +411,14 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GenericFluidSystemNoWater, Scalar, ScalarTypes) using Evaluation = Opm::DenseAd::Evaluation; using FluidSystem = Opm::GenericOilGasWaterFluidSystem; - registerComponent>(); - registerComponent>(); - registerComponent>(); + const auto initializeFluidSystem = [] { + FluidSystem::init(); + registerComponent>(); + registerComponent>(); + registerComponent>(); + }; - checkFluidSystem(); - checkFluidSystem(); - checkFluidSystem(); + checkFluidSystem(initializeFluidSystem); + checkFluidSystem(initializeFluidSystem); + checkFluidSystem(initializeFluidSystem); } From 4e2ceb3015e9af3b83014018bec039072aa73736 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 8 Sep 2026 22:07:40 +0200 Subject: [PATCH 4/5] Reset compositional interaction and viscosity coefficients Clear BIC values and restore the default LBC coefficients in init() so manual component registration cannot inherit settings from an earlier deck. Add a regression test that loads custom coefficients and checks the defaults after reinitialization. Document the unshifted fugacity convention used by the equal-pressure PT flash and distinguish it from absolute translated fugacity. Validation: six focused opm-common tests pass. Restoring the previous reset behavior makes the new regression test fail. --- .../GenericOilGasWaterFluidSystem.hpp | 12 +++- tests/material/test_volume_shift.cpp | 59 +++++++++++++++++-- 2 files changed, 66 insertions(+), 5 deletions(-) diff --git a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp index 9332acd6bc8..ae88039f7f7 100644 --- a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp +++ b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp @@ -209,6 +209,8 @@ namespace Opm { // registrations replace it. component_param_.clear(); component_param_.reserve(numComponents); + interaction_coefficients_.clear(); + lbc_coefficients_ = ViscosityModel::defaultLBCCoefficients(); } /*! @@ -391,7 +393,15 @@ namespace Opm { } } - //! \copydoc BaseFluidSystem::fugacityCoefficient + /*! + * \copydoc BaseFluidSystem::fugacityCoefficient + * + * Returns the coefficient of the unshifted cubic EOS. The PT flash + * uses these coefficients at equal phase pressure and temperature, + * where the SSHIFT translation factors cancel from their ratios. + * An absolute coefficient of the translated EOS additionally requires + * the factor exp(-p s_c b_c / (R T)), which is not included here. + */ template static LhsEval fugacityCoefficient(const FluidState& fluidState, const ParameterCache& paramCache, diff --git a/tests/material/test_volume_shift.cpp b/tests/material/test_volume_shift.cpp index 94018769a48..78fd8bcbb66 100644 --- a/tests/material/test_volume_shift.cpp +++ b/tests/material/test_volume_shift.cpp @@ -375,11 +375,11 @@ BOOST_AUTO_TEST_CASE(TheShiftedPathCarriesItsDerivatives) // Separate static component data for the deck-initialization test. using DeckSystem = Opm::GenericOilGasWaterFluidSystem; -BOOST_AUTO_TEST_CASE(ParsedShiftReachesTheFluidSystem) +namespace { + +Opm::Deck deckWithVolumeShift() { - // Verify the handoff from parsed SSHIFT values to the component parameters - // through initFromState(). The other tests register components directly. - const auto deck = Opm::Parser{}.parseString(R"( + return Opm::Parser{}.parseString(R"( RUNSPEC METRIC DIMENS @@ -413,13 +413,25 @@ MW 16.043 142.285 44.010 / ACF 0.008 0.4885 0.225 / +BIC + 0.12 0.23 0.34 / +LBCCOEF + 0.2 0.03 0.06 -0.04 0.009 / SSHIFT -0.1595 0.10784 -0.0817 / SOLUTION SCHEDULE END )"); +} +} // Anonymous namespace + +BOOST_AUTO_TEST_CASE(ParsedShiftReachesTheFluidSystem) +{ + // Verify the handoff from parsed SSHIFT values to the component parameters + // through initFromState(). The other tests register components directly. + const auto deck = deckWithVolumeShift(); const auto eclState = Opm::EclipseState{ deck }; const auto schedule = Opm::Schedule{ deck, eclState }; @@ -431,6 +443,45 @@ END } } +BOOST_AUTO_TEST_CASE(ReinitializationRestoresDefaultCoefficients) +{ + const auto deck = deckWithVolumeShift(); + const auto eclState = Opm::EclipseState{ deck }; + const auto schedule = Opm::Schedule{ deck, eclState }; + DeckSystem::initFromState(eclState, schedule); + + const std::array deckLbc{0.2, 0.03, 0.06, -0.04, 0.009}; + const auto& loadedLbc = DeckSystem::lbcCoefficients(); + BOOST_REQUIRE_EQUAL_COLLECTIONS(loadedLbc.begin(), loadedLbc.end(), + deckLbc.begin(), deckLbc.end()); + BOOST_REQUIRE_EQUAL(DeckSystem::interactionCoefficient(0, 1), 0.12); + BOOST_REQUIRE_EQUAL(DeckSystem::interactionCoefficient(0, 2), 0.23); + BOOST_REQUIRE_EQUAL(DeckSystem::interactionCoefficient(1, 2), 0.34); + + // Manual registration after init() must use defaults for properties that + // are not supplied by ComponentParam, regardless of the previous deck. + DeckSystem::init(); + const auto& config = eclState.compositionalConfig(); + const auto& props = config.eosProps(0); + for (unsigned c = 0; c < DeckSystem::numComponents; ++c) { + DeckSystem::addComponent(DeckSystem::ComponentParam{ + config.compName()[c], props.molecular_weights[c], + props.critical_temperature[c], props.critical_pressure[c], + props.critical_volume[c] * 1.e3, props.acentric_factors[c]}); + } + + for (unsigned c = 0; c < DeckSystem::numComponents; ++c) { + BOOST_CHECK_EQUAL(DeckSystem::volumeShift(c), 0.0); + for (unsigned d = 0; d < DeckSystem::numComponents; ++d) { + BOOST_CHECK_EQUAL(DeckSystem::interactionCoefficient(c, d), 0.0); + } + } + const auto defaultLbc = DeckSystem::ViscosityModel::defaultLBCCoefficients(); + const auto& resetLbc = DeckSystem::lbcCoefficients(); + BOOST_CHECK_EQUAL_COLLECTIONS(resetLbc.begin(), resetLbc.end(), + defaultLbc.begin(), defaultLbc.end()); +} + BOOST_AUTO_TEST_CASE(AnOverlargeShiftIsRejectedRatherThanReturned) { // A translation larger than the EOS molar volume must raise an error From 8ad0c1e62cf5f15f76b407a01553325598dc4ae0 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 8 Sep 2026 23:06:40 +0200 Subject: [PATCH 5/5] Translate the fugacity coefficient with SSHIFT The cubic EOS is evaluated on the unshifted root, and the translation enters as the per-component factor exp(-p s_c b_c / (R T)), written as exp(-s_c B_c) since B_c is that group. Without it the coefficient was the untranslated one, which is only usable through ratios. The factor is equal in both phases at equal phase pressure and temperature, so it cancels from the equilibrium ratios and the flash, which converges on a fugacity ratio, is unaffected. A deck without SSHIFT multiplies by exp(0), so its results are unchanged bit for bit. Add a test on an individual coefficient, since the ratio test cannot see the factor it cancels. --- .../GenericOilGasWaterFluidSystem.hpp | 18 ++++++--- tests/material/test_volume_shift.cpp | 39 +++++++++++++++++++ 2 files changed, 51 insertions(+), 6 deletions(-) diff --git a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp index ae88039f7f7..d78bb79b884 100644 --- a/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp +++ b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp @@ -396,11 +396,12 @@ namespace Opm { /*! * \copydoc BaseFluidSystem::fugacityCoefficient * - * Returns the coefficient of the unshifted cubic EOS. The PT flash - * uses these coefficients at equal phase pressure and temperature, - * where the SSHIFT translation factors cancel from their ratios. - * An absolute coefficient of the translated EOS additionally requires - * the factor exp(-p s_c b_c / (R T)), which is not included here. + * The cubic EOS is evaluated on the unshifted root, and the SSHIFT + * translation enters as the per-component factor exp(-p s_c b_c / + * (R T)), written here as exp(-s_c B_c) since B_c = b_c p / (R T). + * The factor is the same in both phases at equal phase pressure and + * temperature, so it cancels from the equilibrium ratios and leaves + * the phase split untouched. */ template static LhsEval fugacityCoefficient(const FluidState& fluidState, @@ -415,7 +416,12 @@ namespace Opm { assert(phaseIdx < numPhases); assert(compIdx < numComponents); - return decay(CubicEOS::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx)); + const auto fugCoeff = + CubicEOS::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); + const auto translation = + exp(-volumeShift(compIdx) * paramCache.Bi(phaseIdx, compIdx)); + + return decay(fugCoeff * translation); } // TODO: the following interfaces are needed by function checkFluidSystem() diff --git a/tests/material/test_volume_shift.cpp b/tests/material/test_volume_shift.cpp index 78fd8bcbb66..e89246ec291 100644 --- a/tests/material/test_volume_shift.cpp +++ b/tests/material/test_volume_shift.cpp @@ -244,6 +244,45 @@ BOOST_AUTO_TEST_CASE(ShiftLeavesTheEquilibriumRatiosAlone) } } +BOOST_AUTO_TEST_CASE(ShiftTranslatesTheFugacityCoefficient) +{ + // The ratio test below only sees the factor cancel. This pins the factor + // itself: the coefficient of the translated EOS is the unshifted one times + // exp(-s_c B_c), so a component with a nonzero shift must carry it. + Opm::CompositionalFluidState fs; + Opm::CompositionalFluidState fsRef; + fs.setTemperature(temperature); + fsRef.setTemperature(temperature); + fs.setPressure(FluidSystem::gasPhaseIdx, pressure); + fsRef.setPressure(UnshiftedSystem::gasPhaseIdx, pressure); + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(FluidSystem::gasPhaseIdx, c, z[c]); + fsRef.setMoleFraction(UnshiftedSystem::gasPhaseIdx, c, z[c]); + } + + typename FluidSystem::template ParameterCache pc(eosType); + typename UnshiftedSystem::template ParameterCache pcRef(eosType); + pc.updatePhase(fs, FluidSystem::gasPhaseIdx); + pcRef.updatePhase(fsRef, UnshiftedSystem::gasPhaseIdx); + + bool sawAShift = false; + for (int c = 0; c < numComponents; ++c) { + const Scalar shifted = + FluidSystem::fugacityCoefficient(fs, pc, FluidSystem::gasPhaseIdx, c); + const Scalar unshifted = + UnshiftedSystem::fugacityCoefficient(fsRef, pcRef, UnshiftedSystem::gasPhaseIdx, c); + const Scalar factor = + std::exp(-FluidSystem::volumeShift(c) * pc.Bi(FluidSystem::gasPhaseIdx, c)); + + BOOST_CHECK_CLOSE(shifted, unshifted * factor, 1.0e-8); + if (std::abs(factor - 1.0) > 1.0e-6) { + sawAShift = true; + } + } + // Otherwise every factor is one and the check above proves nothing. + BOOST_CHECK(sawAShift); +} + BOOST_AUTO_TEST_CASE(ShiftDoesNotReachTheCachedVolume) { // Check that the corrected volume subtracts the cached translation and