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/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/GenericOilGasWaterFluidSystem.hpp b/opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp index 7d01b3f9a51..d78bb79b884 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; @@ -199,7 +205,12 @@ 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); + interaction_coefficients_.clear(); + lbc_coefficients_ = ViscosityModel::defaultLBCCoefficients(); } /*! @@ -220,6 +231,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 +353,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 +379,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)); @@ -358,7 +393,16 @@ namespace Opm { } } - //! \copydoc BaseFluidSystem::fugacityCoefficient + /*! + * \copydoc BaseFluidSystem::fugacityCoefficient + * + * 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, const ParameterCache& paramCache, @@ -372,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/opm/material/fluidsystems/PTFlashParameterCache.hpp b/opm/material/fluidsystems/PTFlashParameterCache.hpp index 85b2fb2a524..8b4a5deefcf 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,30 @@ 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 + { + Scalar shift = 0; + // 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; + } + 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_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); } diff --git a/tests/material/test_volume_shift.cpp b/tests/material/test_volume_shift.cpp new file mode 100644 index 00000000000..e89246ec291 --- /dev/null +++ b/tests/material/test_volume_shift.cpp @@ -0,0 +1,546 @@ +// -*- 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(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 + // 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; + +namespace { + +Opm::Deck deckWithVolumeShift() +{ + return 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 / +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 }; + + 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(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 + // 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); +}