From df6bdbe70087d3d849215432d33a7b346a0416f5 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Wed, 8 Jan 2025 13:08:17 +0100 Subject: [PATCH 1/9] Added dummy water phase to compositional simulator. Note: compositional is now always three phase! --- flowexperimental/comp/flowexp_comp.hpp | 4 +- opm/models/ptflash/flashindices.hh | 14 +++-- .../ptflash/flashintensivequantities.hh | 21 ++++---- opm/models/ptflash/flashlocalresidual.hh | 52 ++++++++++++++----- opm/models/ptflash/flashnewtonmethod.hh | 6 +++ opm/models/ptflash/flashprimaryvariables.hh | 16 +++++- opm/simulators/flow/FlowProblemComp.hpp | 7 +++ .../flow/GenericOutputBlackoilModule.cpp | 4 +- 8 files changed, 92 insertions(+), 32 deletions(-) diff --git a/flowexperimental/comp/flowexp_comp.hpp b/flowexperimental/comp/flowexp_comp.hpp index 213d63c07..44591c560 100644 --- a/flowexperimental/comp/flowexp_comp.hpp +++ b/flowexperimental/comp/flowexp_comp.hpp @@ -20,7 +20,7 @@ #define FLOWEXP_COMP_HPP #include -#include +#include #include #include @@ -169,7 +169,7 @@ private: static constexpr int num_comp = getPropValue(); public: - using type = Opm::GenericOilGasFluidSystem; + using type = Opm::GenericOilGasWaterFluidSystem; }; template struct EnableMech> { diff --git a/opm/models/ptflash/flashindices.hh b/opm/models/ptflash/flashindices.hh index 578b62d76..9027e4334 100644 --- a/opm/models/ptflash/flashindices.hh +++ b/opm/models/ptflash/flashindices.hh @@ -51,13 +51,16 @@ class FlashIndices using EnergyIndices = Opm::EnergyIndices; public: - static constexpr bool waterEnabled = false; + //! All phases active (note: immiscible/"dummy" water phase) + static constexpr bool waterEnabled = true; static constexpr bool gasEnabled = true; static constexpr bool oilEnabled = true; - static constexpr int waterPhaseIdx = -1; - static constexpr int numPhases = 2; + + //! number of active phases + static constexpr int numPhases = 3; + //! number of equations/primary variables - static const int numEq = numComponents + EnergyIndices::numEq_; + static const int numEq = numComponents + EnergyIndices::numEq_ + 1; // Primary variable indices @@ -66,6 +69,9 @@ public: //! Index of the molefraction of the first component static constexpr int z0Idx = pressure0Idx + 1; + + //! Index of water saturation + static constexpr int water0Idx = z0Idx + numComponents - 1; // equation indices diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh index 7e71ef5a3..a7038898d 100644 --- a/opm/models/ptflash/flashintensivequantities.hh +++ b/opm/models/ptflash/flashintensivequantities.hh @@ -76,6 +76,7 @@ class FlashIntensiveQuantities enum { enableEnergy = getPropValue() }; enum { dimWorld = GridView::dimensionworld }; enum { pressure0Idx = Indices::pressure0Idx }; + enum { water0Idx = Indices::water0Idx}; using Scalar = GetPropType; using Evaluation = GetPropType; @@ -215,20 +216,22 @@ public: (R * fluidState_.temperature(FluidSystem::gasPhaseIdx)); - // Update saturation - // \Note: the current implementation assume oil-gas system. + // Update oil/gas saturation; water saturation is a primary variable + Evaluation Sw = priVars.makeEvaluation(water0Idx, timeIdx); Evaluation L = fluidState_.L(); - Evaluation So = Opm::max((L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0); - Evaluation Sg = Opm::max(1 - So, 0.0); - Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg); + Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0); + Evaluation Sg = Opm::max(1 - So - Sw, 0.0); + Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg) + Opm::getValue(Sw); So /= sumS; Sg /= sumS; + Sw /= sumS; - fluidState_.setSaturation(0, So); - fluidState_.setSaturation(1, Sg); + fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So); + fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg); + fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw); - fluidState_.setCompressFactor(0, Z_L); - fluidState_.setCompressFactor(1, Z_V); + fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L); + fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V); // Print saturation if (flashVerbosity >= 5) { diff --git a/opm/models/ptflash/flashlocalresidual.hh b/opm/models/ptflash/flashlocalresidual.hh index c528189e7..8d6172fc1 100644 --- a/opm/models/ptflash/flashlocalresidual.hh +++ b/opm/models/ptflash/flashlocalresidual.hh @@ -49,12 +49,16 @@ class FlashLocalResidual: public GetPropType; using IntensiveQuantities = GetPropType; using ElementContext = GetPropType; + using FluidSystem = GetPropType; enum { numEq = getPropValue() }; enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; + enum { water0Idx = Indices::water0Idx }; enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + enum { enableDiffusion = getPropValue() }; using DiffusionModule = Opm::DiffusionModule; @@ -77,15 +81,25 @@ public: const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); const auto& fs = intQuants.fluidState(); - // compute storage term of all components within all phases - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - unsigned eqIdx = conti0EqIdx + compIdx; - storage[eqIdx] += - Toolbox::template decay(fs.massFraction(phaseIdx, compIdx)) - * Toolbox::template decay(fs.density(phaseIdx)) + // compute water storage term + if (phaseIdx == waterPhaseIdx) { + unsigned eqIdx = conti0EqIdx + numComponents; + storage[eqIdx] = + Toolbox::template decay(fs.density(phaseIdx)) * Toolbox::template decay(fs.saturation(phaseIdx)) * Toolbox::template decay(intQuants.porosity()); } + else { + // compute storage term of all components within oil/gas phases + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + unsigned eqIdx = conti0EqIdx + compIdx; + storage[eqIdx] += + Toolbox::template decay(fs.massFraction(phaseIdx, compIdx)) + * Toolbox::template decay(fs.density(phaseIdx)) + * Toolbox::template decay(fs.saturation(phaseIdx)) + * Toolbox::template decay(intQuants.porosity()); + } + } EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx); } @@ -146,19 +160,31 @@ public: up.fluidState().density(phaseIdx) * extQuants.volumeFlux(phaseIdx); - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - flux[conti0EqIdx + compIdx] += - tmp*up.fluidState().massFraction(phaseIdx, compIdx); + if (phaseIdx == waterPhaseIdx) { + unsigned eqIdx = conti0EqIdx + numComponents; + flux[eqIdx] = tmp; + } + else { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + flux[conti0EqIdx + compIdx] += + tmp*up.fluidState().massFraction(phaseIdx, compIdx); + } } } else { Evaluation tmp = Toolbox::value(up.fluidState().density(phaseIdx)) * extQuants.volumeFlux(phaseIdx); - - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - flux[conti0EqIdx + compIdx] += - tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx)); + + if (phaseIdx == waterPhaseIdx) { + unsigned eqIdx = conti0EqIdx + numComponents; + flux[eqIdx] = tmp; + } + else { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + flux[conti0EqIdx + compIdx] += + tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx)); + } } } } diff --git a/opm/models/ptflash/flashnewtonmethod.hh b/opm/models/ptflash/flashnewtonmethod.hh index 32500f38b..093f4a6f7 100644 --- a/opm/models/ptflash/flashnewtonmethod.hh +++ b/opm/models/ptflash/flashnewtonmethod.hh @@ -125,6 +125,12 @@ protected: for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol); } + + // limit change in water saturation to 0.2 + constexpr Scalar dSwMax = 0.2; + if (update[Indices::water0Idx] > dSwMax) { + nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax; + } } private: void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const diff --git a/opm/models/ptflash/flashprimaryvariables.hh b/opm/models/ptflash/flashprimaryvariables.hh index 2df1c67c3..7f9914327 100644 --- a/opm/models/ptflash/flashprimaryvariables.hh +++ b/opm/models/ptflash/flashprimaryvariables.hh @@ -61,6 +61,12 @@ class FlashPrimaryVariables : public FvBasePrimaryVariables // primary variable indices enum { z0Idx = Indices::z0Idx }; enum { pressure0Idx = Indices::pressure0Idx }; + enum { water0Idx = Indices::water0Idx }; + + // phase indices + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; + enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; @@ -108,10 +114,15 @@ public: // the energy module EnergyModule::setPriVarTemperatures(*this, fluidState); + // assign components total fraction for (int i = 0; i < numComponents - 1; ++i) (*this)[z0Idx + i] = getValue(fluidState.moleFraction(i)); - (*this)[pressure0Idx] = getValue(fluidState.pressure(0)); + // assign pressure + (*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx)); + + // assign water saturation + (*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx)); } /*! @@ -121,12 +132,13 @@ public: */ void print(std::ostream& os = std::cout) const { - os << "(p_" << FluidSystem::phaseName(0) << " = " + os << "(p_" << FluidSystem::phaseName(FluidSystem::oilPhaseIdx) << " = " << this->operator[](pressure0Idx); for (unsigned compIdx = 0; compIdx < numComponents - 2; ++compIdx) { os << ", z_" << FluidSystem::componentName(compIdx) << " = " << this->operator[](z0Idx + compIdx); } + os << ", S_w = " << this->operator[](water0Idx); os << ")" << std::flush; } }; diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index be6a2616b..575423d32 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -362,6 +362,9 @@ public: Dune::FieldVector z(0.0); Scalar sumMoles = 0.0; for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (phaseIdx == waterPhaseIdx){ + continue; + } const auto saturation = fs.saturation(phaseIdx); for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { Scalar tmp = fs.molarity(phaseIdx, compIdx) * saturation; @@ -527,6 +530,10 @@ protected: - waterSaturationData[dofIdx] - gasSaturationData[dofIdx]); } + if (water_active) { + dofFluidState.setSaturation(FluidSystem::waterPhaseIdx, + waterSaturationData[dofIdx]); + } ////// // set phase pressures diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index abd8f3cc3..2d6847fad 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -30,7 +30,7 @@ #include #include -#include +#include #include #include @@ -1677,7 +1677,7 @@ INSTANTIATE_TYPE(float) #endif #define INSTANTIATE_COMP(NUM) \ - template using FS##NUM = GenericOilGasFluidSystem; \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ template class GenericOutputBlackoilModule>; INSTANTIATE_COMP(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput From fb403d05ae99981401286690fa72e3c510263bca Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Wed, 8 Jan 2025 13:12:54 +0100 Subject: [PATCH 2/9] Changed to three phase using new GenericOilGasWaterFluidSystem --- examples/problems/co2ptflashproblem.hh | 45 +++++++++++++++++++------- 1 file changed, 34 insertions(+), 11 deletions(-) diff --git a/examples/problems/co2ptflashproblem.hh b/examples/problems/co2ptflashproblem.hh index 789e2c40a..a75036acb 100644 --- a/examples/problems/co2ptflashproblem.hh +++ b/examples/problems/co2ptflashproblem.hh @@ -38,7 +38,8 @@ #include #include #include -#include +#include +#include #include #include #include @@ -106,7 +107,7 @@ private: static constexpr int num_comp = getPropValue(); public: - using type = Opm::GenericOilGasFluidSystem; + using type = Opm::GenericOilGasWaterFluidSystem; }; // Set the material Law @@ -118,8 +119,8 @@ private: enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; using Scalar = GetPropType; - using Traits = Opm::TwoPhaseMaterialTraits; @@ -195,6 +196,8 @@ class CO2PTProblem : public GetPropType enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { pressure0Idx = Indices::pressure0Idx }; + enum { z0Idx = Indices::z0Idx }; enum { numComponents = getPropValue() }; enum { enableEnergy = getPropValue() }; enum { enableDiffusion = getPropValue() }; @@ -237,6 +240,21 @@ public: porosity_ = 0.1; } + void initWaterPVT() + { + using WaterPvt = typename FluidSystem::WaterPvt; + std::shared_ptr waterPvt = std::make_shared(); + waterPvt->setApproach(WaterPvtApproach::ConstantCompressibilityWater); + auto& ccWaterPvt = waterPvt->template getRealPvt(); + ccWaterPvt.setNumRegions(/*numPvtRegions=*/1); + Scalar rhoRefW = 1037.0; // [kg] + ccWaterPvt.setReferenceDensities(/*regionIdx=*/0, /*rhoRefO=*/Scalar(0.0), /*rhoRefG=*/Scalar(0.0), rhoRefW); + ccWaterPvt.setViscosity(/*regionIdx=*/0, 9.6e-4); + ccWaterPvt.setCompressibility(/*regionIdx=*/0, 1.450377e-10); + waterPvt->initEnd(); + FluidSystem::setWaterPvt(std::move(waterPvt)); + } + template const DimVector& gravity([[maybe_unused]]const Context& context, @@ -262,8 +280,12 @@ public: void finishInit() { ParentType::finishInit(); + // initialize fixed parameters; temperature, permeability, porosity initPetrophysics(); + + // Initialize water pvt + initWaterPVT(); } /*! @@ -358,8 +380,8 @@ public: // Calculate storage terms PrimaryVariables storageL, storageG; - this->model().globalPhaseStorage(storageL, /*phaseIdx=*/0); - this->model().globalPhaseStorage(storageG, /*phaseIdx=*/1); + this->model().globalPhaseStorage(storageL, /*phaseIdx=*/oilPhaseIdx); + this->model().globalPhaseStorage(storageG, /*phaseIdx=*/gasPhaseIdx); // Write mass balance information for rank 0 // if (this->gridView().comm().rank() == 0) { @@ -456,12 +478,12 @@ private: int prod = Parameters::Get() - 1; int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx); ComponentVector comp; - comp[0] = Evaluation::createVariable(0.5, 1); - comp[1] = Evaluation::createVariable(0.3, 2); + comp[0] = Evaluation::createVariable(0.5, z0Idx); + comp[1] = Evaluation::createVariable(0.3, z0Idx + 1); comp[2] = 1. - comp[0] - comp[1]; if (spatialIdx == inj) { - comp[0] = Evaluation::createVariable(0.99, 1); - comp[1] = Evaluation::createVariable(0.01 - 1e-3, 2); + comp[0] = Evaluation::createVariable(0.99, z0Idx); + comp[1] = Evaluation::createVariable(0.01 - 1e-3, z0Idx + 1); comp[2] = 1. - comp[0] - comp[1]; } @@ -475,10 +497,11 @@ private: if (spatialIdx == prod) { p0 *= 0.5; } - Evaluation p_init = Evaluation::createVariable(p0, 0); + Evaluation p_init = Evaluation::createVariable(p0, pressure0Idx); fs.setPressure(FluidSystem::oilPhaseIdx, p_init); fs.setPressure(FluidSystem::gasPhaseIdx, p_init); + fs.setPressure(FluidSystem::waterPhaseIdx, p_init); fs.setTemperature(temperature_); From 2f0443947d998579304ec9d477f3a2e0d732c92b Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Mon, 20 Jan 2025 12:51:23 +0100 Subject: [PATCH 3/9] Two- and three-phase compositional simulators. Choice of simulator is based on if WATER is present in RUNSPEC. --- CMakeLists.txt | 1 + flowexperimental/comp/flowexp_comp.cpp | 22 ++-- flowexperimental/comp/flowexp_comp.hpp | 113 +++++++++++---------- flowexperimental/comp/flowexp_comp2.cpp | 4 +- flowexperimental/comp/flowexp_comp2_2p.cpp | 36 +++++++ flowexperimental/comp/flowexp_comp3.cpp | 4 +- flowexperimental/comp/flowexp_comp3_2p.cpp | 36 +++++++ flowexperimental/comp/flowexp_comp4.cpp | 4 +- flowexperimental/comp/flowexp_comp4_2p.cpp | 36 +++++++ flowexperimental/comp/flowexp_comp5.cpp | 4 +- flowexperimental/comp/flowexp_comp5_2p.cpp | 36 +++++++ flowexperimental/comp/flowexp_comp6.cpp | 4 +- flowexperimental/comp/flowexp_comp6_2p.cpp | 36 +++++++ flowexperimental/comp/flowexp_comp7.cpp | 4 +- flowexperimental/comp/flowexp_comp7_2p.cpp | 36 +++++++ 15 files changed, 306 insertions(+), 70 deletions(-) create mode 100644 flowexperimental/comp/flowexp_comp2_2p.cpp create mode 100644 flowexperimental/comp/flowexp_comp3_2p.cpp create mode 100644 flowexperimental/comp/flowexp_comp4_2p.cpp create mode 100644 flowexperimental/comp/flowexp_comp5_2p.cpp create mode 100644 flowexperimental/comp/flowexp_comp6_2p.cpp create mode 100644 flowexperimental/comp/flowexp_comp7_2p.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2d0ffa034..45c8fa34e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -625,6 +625,7 @@ opm_add_test(flowexp_blackoil set(FLOWEXP_COMPONENTS_SOURCES) foreach(component IN LISTS OPM_COMPILE_COMPONENTS) list(APPEND FLOWEXP_COMPONENTS_SOURCES flowexperimental/comp/flowexp_comp${component}.cpp) + list(APPEND FLOWEXP_COMPONENTS_SOURCES flowexperimental/comp/flowexp_comp${component}_2p.cpp) endforeach() opm_add_test(flowexp_comp diff --git a/flowexperimental/comp/flowexp_comp.cpp b/flowexperimental/comp/flowexp_comp.cpp index 878811398..afda32676 100644 --- a/flowexperimental/comp/flowexp_comp.cpp +++ b/flowexperimental/comp/flowexp_comp.cpp @@ -36,10 +36,13 @@ template std::tuple -runComponent(int runtimeComponent, int argc, char** argv) +runComponent(int runtimeComponent, bool water, int argc, char** argv) { if (runtimeComponent == compileTimeComponent) { - return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); + if (water) + return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); + else + return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); } return std::make_tuple(false, EXIT_FAILURE); } @@ -63,18 +66,21 @@ runComponent(int runtimeComponent, int argc, char** argv) */ template std::tuple -runComponent(int runtimecomponent, int argc, char** argv) +runComponent(int runtimecomponent, bool water, int argc, char** argv) { if (currentCompileTimeComponent == runtimecomponent) { - return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); + if (water) + return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); + else + return std::make_tuple(true, Opm::dispatchFlowExpComp(argc, argv)); } - return runComponent(runtimecomponent, argc, argv); + return runComponent(runtimecomponent, water, argc, argv); } int main(int argc, char** argv) { - using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<0>; + using TypeTag = Opm::Properties::TTag::FlowExpCompProblem<0, true>; Opm::registerEclTimeSteppingParameters(); // At the moment, this is probably as optimal as can be. @@ -90,9 +96,11 @@ main(int argc, char** argv) = Opm::Parser {}.parseFile(inputFilename, Opm::ParseContext {}, std::vector {Opm::Ecl::SectionType::RUNSPEC}); const auto runspec = Opm::Runspec(deck); const auto numComps = runspec.numComps(); + const auto& phases = runspec.phases(); + const auto wat = phases.active(Opm::Phase::WATER); auto [componentSupported, executionStatus] - = runComponent(numComps, argc, argv); + = runComponent(numComps, wat, argc, argv); if (!componentSupported) { fmt::print("Deck has {} components, not supported. In this build of the simulator, we support the " diff --git a/flowexperimental/comp/flowexp_comp.hpp b/flowexperimental/comp/flowexp_comp.hpp index 44591c560..02920258e 100644 --- a/flowexperimental/comp/flowexp_comp.hpp +++ b/flowexperimental/comp/flowexp_comp.hpp @@ -39,7 +39,7 @@ // suggestTimeStep is taken from newton solver in problem.limitTimestep namespace Opm { -template +template int dispatchFlowExpComp(int argc, char** argv); } @@ -47,15 +47,15 @@ int dispatchFlowExpComp(int argc, char** argv); namespace Opm::Properties { namespace TTag { -template +template struct FlowExpCompProblem { using InheritsFrom = std::tuple; }; } -template -struct SparseMatrixAdapter> +template +struct SparseMatrixAdapter> { private: using Scalar = GetPropType; @@ -109,30 +109,30 @@ struct LocalLinearizerSplice #endif // Set the problem property -template -struct Problem> +template +struct Problem> { using type = FlowProblemComp; }; -template -struct AquiferModel> { +template +struct AquiferModel> { using type = EmptyModel; }; -template -struct WellModel> { +template +struct WellModel> { using type = EmptyModel; }; -template -struct TracerModel> { +template +struct TracerModel> { using type = EmptyModel; }; -template -struct FlashSolver> { +template +struct FlashSolver> { private: using Scalar = GetPropType; using FluidSystem = GetPropType; @@ -147,10 +147,18 @@ template struct NumComp { using type = UndefinedProperty; }; // TODO: this is unfortunate, have to check why we need to hard-code it -template -struct NumComp> { +template +struct NumComp> { static constexpr int value = NumComp_; }; + +template +struct EnableDummyWater { using type = UndefinedProperty; }; + +template +struct EnableDummyWater> { + static constexpr bool value = EnableWater_; +}; #if 0 struct Temperature { using type = UndefinedProperty; }; @@ -161,26 +169,29 @@ struct Temperature { using type = UndefinedProperty; }; }; #endif -template -struct FluidSystem> +template +struct FluidSystem> { private: using Scalar = GetPropType; static constexpr int num_comp = getPropValue(); + static constexpr bool enable_water = getPropValue(); public: - using type = Opm::GenericOilGasWaterFluidSystem; + using type = Opm::GenericOilGasWaterFluidSystem; }; -template -struct EnableMech> { +template +struct EnableMech> { static constexpr bool value = false; }; -template -struct EnableDisgasInWater> { static constexpr bool value = false; }; +template +struct EnableDisgasInWater> { + static constexpr bool value = false; +}; -template -struct Stencil> +template +struct Stencil> { private: using Scalar = GetPropType; @@ -190,62 +201,62 @@ public: using type = EcfvStencil; }; -template -struct EnableApiTracking> { +template +struct EnableApiTracking> { static constexpr bool value = false; }; -template -struct EnableTemperature> { +template +struct EnableTemperature> { static constexpr bool value = false; }; -template -struct EnableSaltPrecipitation> { +template +struct EnableSaltPrecipitation> { static constexpr bool value = false; }; -template -struct EnablePolymerMW> { +template +struct EnablePolymerMW> { static constexpr bool value = false; }; -template -struct EnablePolymer> { +template +struct EnablePolymer> { static constexpr bool value = false; }; -template -struct EnableDispersion> { +template +struct EnableDispersion> { static constexpr bool value = false; }; -template -struct EnableBrine> { +template +struct EnableBrine> { static constexpr bool value = false; }; -template -struct EnableVapwat> { +template +struct EnableVapwat> { static constexpr bool value = false; }; -template -struct EnableSolvent> { +template +struct EnableSolvent> { static constexpr bool value = false; }; -template -struct EnableEnergy> { +template +struct EnableEnergy> { static constexpr bool value = false; }; -template -struct EnableFoam> { +template +struct EnableFoam> { static constexpr bool value = false; }; -template -struct EnableExtbo> { +template +struct EnableExtbo> { static constexpr bool value = false; }; -template -struct EnableMICP> { +template +struct EnableMICP> { static constexpr bool value = false; }; diff --git a/flowexperimental/comp/flowexp_comp2.cpp b/flowexperimental/comp/flowexp_comp2.cpp index 6886e0870..1da968da6 100644 --- a/flowexperimental/comp/flowexp_comp2.cpp +++ b/flowexperimental/comp/flowexp_comp2.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<2>(int argc, char** argv) +int dispatchFlowExpComp<2, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp2_2p.cpp b/flowexperimental/comp/flowexp_comp2_2p.cpp new file mode 100644 index 000000000..984bdb1a7 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp2_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, 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 3 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 . +*/ + +#include "config.h" + +#include + +#include + +#include "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<2, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp3.cpp b/flowexperimental/comp/flowexp_comp3.cpp index 0e5725092..0afadad14 100644 --- a/flowexperimental/comp/flowexp_comp3.cpp +++ b/flowexperimental/comp/flowexp_comp3.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<3>(int argc, char** argv) +int dispatchFlowExpComp<3, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp3_2p.cpp b/flowexperimental/comp/flowexp_comp3_2p.cpp new file mode 100644 index 000000000..d89d51533 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp3_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, 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 3 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 . +*/ + +#include "config.h" + +#include + +#include + +#include "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<3, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp4.cpp b/flowexperimental/comp/flowexp_comp4.cpp index ddb834075..7c018af27 100644 --- a/flowexperimental/comp/flowexp_comp4.cpp +++ b/flowexperimental/comp/flowexp_comp4.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<4>(int argc, char** argv) +int dispatchFlowExpComp<4, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp4_2p.cpp b/flowexperimental/comp/flowexp_comp4_2p.cpp new file mode 100644 index 000000000..68446bc24 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp4_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, 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 3 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 . +*/ + +#include "config.h" + +#include + +#include + +#include "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<4, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp5.cpp b/flowexperimental/comp/flowexp_comp5.cpp index a92539b26..d7b9cf5b7 100644 --- a/flowexperimental/comp/flowexp_comp5.cpp +++ b/flowexperimental/comp/flowexp_comp5.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<5>(int argc, char** argv) +int dispatchFlowExpComp<5, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp5_2p.cpp b/flowexperimental/comp/flowexp_comp5_2p.cpp new file mode 100644 index 000000000..c2c034d8a --- /dev/null +++ b/flowexperimental/comp/flowexp_comp5_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, 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 3 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 . +*/ + +#include "config.h" + +#include + +#include + +#include "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<5, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp6.cpp b/flowexperimental/comp/flowexp_comp6.cpp index b274591dd..95ee13066 100644 --- a/flowexperimental/comp/flowexp_comp6.cpp +++ b/flowexperimental/comp/flowexp_comp6.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<6>(int argc, char** argv) +int dispatchFlowExpComp<6, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp6_2p.cpp b/flowexperimental/comp/flowexp_comp6_2p.cpp new file mode 100644 index 000000000..3d6b88210 --- /dev/null +++ b/flowexperimental/comp/flowexp_comp6_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, 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 3 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 . +*/ + +#include "config.h" + +#include + +#include + +#include "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<6, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} diff --git a/flowexperimental/comp/flowexp_comp7.cpp b/flowexperimental/comp/flowexp_comp7.cpp index a40d2d23b..0c5fcc7f2 100644 --- a/flowexperimental/comp/flowexp_comp7.cpp +++ b/flowexperimental/comp/flowexp_comp7.cpp @@ -28,9 +28,9 @@ namespace Opm { template<> -int dispatchFlowExpComp<7>(int argc, char** argv) +int dispatchFlowExpComp<7, true>(int argc, char** argv) { - return start>(argc, argv, false); + return start>(argc, argv, false); } } diff --git a/flowexperimental/comp/flowexp_comp7_2p.cpp b/flowexperimental/comp/flowexp_comp7_2p.cpp new file mode 100644 index 000000000..0dba55e3b --- /dev/null +++ b/flowexperimental/comp/flowexp_comp7_2p.cpp @@ -0,0 +1,36 @@ +/* + Copyright 2024, 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 3 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 . +*/ + +#include "config.h" + +#include + +#include + +#include "flowexp_comp.hpp" + +namespace Opm { + +template<> +int dispatchFlowExpComp<7, false>(int argc, char** argv) +{ + return start>(argc, argv, false); +} + +} From f59834d3e65a97230ee7078e68e2fb2094819785 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Mon, 20 Jan 2025 12:52:54 +0100 Subject: [PATCH 4/9] Accommodate for two or three phases in compositional simulations. --- examples/problems/co2ptflashproblem.hh | 11 +++++++- opm/models/common/multiphasebaseproperties.hh | 2 ++ opm/models/ptflash/flashindices.hh | 9 +++--- .../ptflash/flashintensivequantities.hh | 16 ++++++++--- opm/models/ptflash/flashlocalresidual.hh | 8 ++++-- opm/models/ptflash/flashmodel.hh | 4 +++ opm/models/ptflash/flashnewtonmethod.hh | 12 +++++--- opm/models/ptflash/flashprimaryvariables.hh | 10 +++++-- opm/simulators/flow/FlowProblemComp.hpp | 4 +-- .../flow/GenericOutputBlackoilModule.cpp | 28 +++++++++++++------ 10 files changed, 75 insertions(+), 29 deletions(-) diff --git a/examples/problems/co2ptflashproblem.hh b/examples/problems/co2ptflashproblem.hh index a75036acb..e7df55c89 100644 --- a/examples/problems/co2ptflashproblem.hh +++ b/examples/problems/co2ptflashproblem.hh @@ -77,6 +77,14 @@ struct NumComp { static constexpr int value = 3; }; +template +struct EnableDummyWater { using type = UndefinedProperty; }; + +template +struct EnableDummyWater { + static constexpr bool value = true; +}; + // Set the grid type: --->2D template struct Grid { using type = Dune::YaspGrid; }; @@ -105,9 +113,10 @@ struct FluidSystem private: using Scalar = GetPropType; static constexpr int num_comp = getPropValue(); + static constexpr bool enable_water = getPropValue(); public: - using type = Opm::GenericOilGasWaterFluidSystem; + using type = Opm::GenericOilGasWaterFluidSystem; }; // Set the material Law diff --git a/opm/models/common/multiphasebaseproperties.hh b/opm/models/common/multiphasebaseproperties.hh index f067612f9..968a750ba 100644 --- a/opm/models/common/multiphasebaseproperties.hh +++ b/opm/models/common/multiphasebaseproperties.hh @@ -83,6 +83,8 @@ struct EnableDispersion { using type = UndefinedProperty; }; //! Enable convective mixing? template struct EnableConvectiveMixing { using type = UndefinedProperty; }; +template +struct EnableWater { using type = UndefinedProperty; }; } // namespace Opm::Properties diff --git a/opm/models/ptflash/flashindices.hh b/opm/models/ptflash/flashindices.hh index 9027e4334..ac4c0ebf4 100644 --- a/opm/models/ptflash/flashindices.hh +++ b/opm/models/ptflash/flashindices.hh @@ -48,19 +48,20 @@ class FlashIndices { static constexpr int numComponents = getPropValue(); enum { enableEnergy = getPropValue() }; + enum { enableWater = getPropValue() }; using EnergyIndices = Opm::EnergyIndices; public: //! All phases active (note: immiscible/"dummy" water phase) - static constexpr bool waterEnabled = true; + static constexpr bool waterEnabled = enableWater; static constexpr bool gasEnabled = true; static constexpr bool oilEnabled = true; //! number of active phases - static constexpr int numPhases = 3; + static constexpr int numPhases = enableWater ? 3 : 2; //! number of equations/primary variables - static const int numEq = numComponents + EnergyIndices::numEq_ + 1; + static const int numEq = numComponents + EnergyIndices::numEq_ + (enableWater ? 1 : 0); // Primary variable indices @@ -71,7 +72,7 @@ public: static constexpr int z0Idx = pressure0Idx + 1; //! Index of water saturation - static constexpr int water0Idx = z0Idx + numComponents - 1; + static constexpr int water0Idx = enableWater ? z0Idx + numComponents - 1 : -1000; // equation indices diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh index a7038898d..4fad4b706 100644 --- a/opm/models/ptflash/flashintensivequantities.hh +++ b/opm/models/ptflash/flashintensivequantities.hh @@ -78,6 +78,8 @@ class FlashIntensiveQuantities enum { pressure0Idx = Indices::pressure0Idx }; enum { water0Idx = Indices::water0Idx}; + static constexpr bool waterEnabled = Indices::waterEnabled; + using Scalar = GetPropType; using Evaluation = GetPropType; using FluidSystem = GetPropType; @@ -216,8 +218,12 @@ public: (R * fluidState_.temperature(FluidSystem::gasPhaseIdx)); - // Update oil/gas saturation; water saturation is a primary variable - Evaluation Sw = priVars.makeEvaluation(water0Idx, timeIdx); + // Update saturation + Evaluation Sw = 0.0; + if constexpr (waterEnabled) { + Sw = priVars.makeEvaluation(water0Idx, timeIdx); + fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw); + } Evaluation L = fluidState_.L(); Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0); Evaluation Sg = Opm::max(1 - So - Sw, 0.0); @@ -228,7 +234,6 @@ public: fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So); fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg); - fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw); fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L); fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V); @@ -253,7 +258,10 @@ public: // set the phase viscosity and density for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - paramCache.updatePhase(fluidState_, phaseIdx); + if (phaseIdx == static_cast(FluidSystem::oilPhaseIdx) + || phaseIdx == static_cast(FluidSystem::gasPhaseIdx)) { + paramCache.updatePhase(fluidState_, phaseIdx); + } const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx); diff --git a/opm/models/ptflash/flashlocalresidual.hh b/opm/models/ptflash/flashlocalresidual.hh index 8d6172fc1..29ac021db 100644 --- a/opm/models/ptflash/flashlocalresidual.hh +++ b/opm/models/ptflash/flashlocalresidual.hh @@ -65,6 +65,8 @@ class FlashLocalResidual: public GetPropType() }; using EnergyModule = Opm::EnergyModule; + static const bool waterEnabled = Indices::waterEnabled; + using Toolbox = Opm::MathToolbox; public: @@ -82,7 +84,7 @@ public: const auto& fs = intQuants.fluidState(); // compute water storage term - if (phaseIdx == waterPhaseIdx) { + if (waterEnabled && phaseIdx == static_cast(waterPhaseIdx)) { unsigned eqIdx = conti0EqIdx + numComponents; storage[eqIdx] = Toolbox::template decay(fs.density(phaseIdx)) @@ -160,7 +162,7 @@ public: up.fluidState().density(phaseIdx) * extQuants.volumeFlux(phaseIdx); - if (phaseIdx == waterPhaseIdx) { + if (waterEnabled && phaseIdx == static_cast(waterPhaseIdx)) { unsigned eqIdx = conti0EqIdx + numComponents; flux[eqIdx] = tmp; } @@ -176,7 +178,7 @@ public: Toolbox::value(up.fluidState().density(phaseIdx)) * extQuants.volumeFlux(phaseIdx); - if (phaseIdx == waterPhaseIdx) { + if (waterEnabled && phaseIdx == static_cast(waterPhaseIdx)) { unsigned eqIdx = conti0EqIdx + numComponents; flux[eqIdx] = tmp; } diff --git a/opm/models/ptflash/flashmodel.hh b/opm/models/ptflash/flashmodel.hh index f348523af..8a192f50e 100644 --- a/opm/models/ptflash/flashmodel.hh +++ b/opm/models/ptflash/flashmodel.hh @@ -135,6 +135,10 @@ template struct EnableEnergy { static constexpr bool value = false; }; +template +struct EnableWater +{ static constexpr int value = GetPropType::waterEnabled; }; + } // namespace Opm::Properties namespace Opm { diff --git a/opm/models/ptflash/flashnewtonmethod.hh b/opm/models/ptflash/flashnewtonmethod.hh index 093f4a6f7..5b6fa1ddb 100644 --- a/opm/models/ptflash/flashnewtonmethod.hh +++ b/opm/models/ptflash/flashnewtonmethod.hh @@ -63,6 +63,8 @@ class FlashNewtonMethod : public GetPropType() }; + static constexpr bool waterEnabled = Indices::waterEnabled; + public: /*! * \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& ) @@ -126,10 +128,12 @@ protected: clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol); } - // limit change in water saturation to 0.2 - constexpr Scalar dSwMax = 0.2; - if (update[Indices::water0Idx] > dSwMax) { - nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax; + if constexpr (waterEnabled) { + // limit change in water saturation to 0.2 + constexpr Scalar dSwMax = 0.2; + if (update[Indices::water0Idx] > dSwMax) { + nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax; + } } } private: diff --git a/opm/models/ptflash/flashprimaryvariables.hh b/opm/models/ptflash/flashprimaryvariables.hh index 7f9914327..897c322ff 100644 --- a/opm/models/ptflash/flashprimaryvariables.hh +++ b/opm/models/ptflash/flashprimaryvariables.hh @@ -63,6 +63,8 @@ class FlashPrimaryVariables : public FvBasePrimaryVariables enum { pressure0Idx = Indices::pressure0Idx }; enum { water0Idx = Indices::water0Idx }; + static constexpr bool waterEnabled = Indices::waterEnabled; + // phase indices enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; @@ -122,7 +124,9 @@ public: (*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx)); // assign water saturation - (*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx)); + if constexpr (waterEnabled) { + (*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx)); + } } /*! @@ -138,7 +142,9 @@ public: os << ", z_" << FluidSystem::componentName(compIdx) << " = " << this->operator[](z0Idx + compIdx); } - os << ", S_w = " << this->operator[](water0Idx); + if constexpr (waterEnabled) { + os << ", S_w = " << this->operator[](water0Idx); + } os << ")" << std::flush; } }; diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index 575423d32..d6e6c817f 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -362,7 +362,7 @@ public: Dune::FieldVector z(0.0); Scalar sumMoles = 0.0; for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (phaseIdx == waterPhaseIdx){ + if (Indices::waterEnabled && phaseIdx == static_cast(waterPhaseIdx)){ continue; } const auto saturation = fs.saturation(phaseIdx); @@ -494,7 +494,7 @@ protected: const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx); const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx); - if (water_active && Indices::numPhases > 1) + if (water_active && Indices::numPhases > 2) waterSaturationData = fp.get_double("SWAT"); else waterSaturationData.resize(numDof); diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 2d6847fad..19e11d4cd 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -1676,17 +1676,27 @@ INSTANTIATE_TYPE(double) INSTANTIATE_TYPE(float) #endif -#define INSTANTIATE_COMP(NUM) \ - template using FS##NUM = GenericOilGasWaterFluidSystem; \ +#define INSTANTIATE_COMP_THREEPHASE(NUM) \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ template class GenericOutputBlackoilModule>; -INSTANTIATE_COMP(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput +INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput +INSTANTIATE_COMP_THREEPHASE(2) +INSTANTIATE_COMP_THREEPHASE(3) +INSTANTIATE_COMP_THREEPHASE(4) +INSTANTIATE_COMP_THREEPHASE(5) +INSTANTIATE_COMP_THREEPHASE(6) +INSTANTIATE_COMP_THREEPHASE(7) -INSTANTIATE_COMP(2) -INSTANTIATE_COMP(3) -INSTANTIATE_COMP(4) -INSTANTIATE_COMP(5) -INSTANTIATE_COMP(6) -INSTANTIATE_COMP(7) +#define INSTANTIATE_COMP_TWOPHASE(NUM) \ + template using GFS##NUM = GenericOilGasWaterFluidSystem; \ + template class GenericOutputBlackoilModule>; + +INSTANTIATE_COMP_TWOPHASE(2) +INSTANTIATE_COMP_TWOPHASE(3) +INSTANTIATE_COMP_TWOPHASE(4) +INSTANTIATE_COMP_TWOPHASE(5) +INSTANTIATE_COMP_TWOPHASE(6) +INSTANTIATE_COMP_TWOPHASE(7) } // namespace Opm From 4be8a5c673c3acc363a59c0507254a522e68ccab Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Wed, 29 Jan 2025 12:44:48 +0100 Subject: [PATCH 5/9] Use fluid system phase indices --- opm/simulators/flow/FlowProblemCompProperties.hpp | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/opm/simulators/flow/FlowProblemCompProperties.hpp b/opm/simulators/flow/FlowProblemCompProperties.hpp index 12ecb2319..7b838a94c 100644 --- a/opm/simulators/flow/FlowProblemCompProperties.hpp +++ b/opm/simulators/flow/FlowProblemCompProperties.hpp @@ -69,17 +69,10 @@ private: using Scalar = GetPropType; using FluidSystem = GetPropType; - // using Traits = ThreePhaseMaterialTraits; - - // TODO: We should be able to use FluidSystem here and using Indices to handle the active phases - // some more development is needed using Traits = ThreePhaseMaterialTraits; + /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx, + /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx, + /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>; public: using EclMaterialLawManager = ::Opm::EclMaterialLawManager; From 20ca15427852bd412370e08aa68b1d003e8c95a1 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Wed, 5 Feb 2025 08:35:18 +0100 Subject: [PATCH 6/9] Change to new compositional fluid system --- opm/simulators/flow/FIPContainer.cpp | 31 +++++++++++++++++++--------- 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/opm/simulators/flow/FIPContainer.cpp b/opm/simulators/flow/FIPContainer.cpp index 3ae9b46d9..00557a1cd 100644 --- a/opm/simulators/flow/FIPContainer.cpp +++ b/opm/simulators/flow/FIPContainer.cpp @@ -27,7 +27,7 @@ #include #include -#include +#include #include @@ -438,16 +438,27 @@ INSTANTIATE_TYPE(double) INSTANTIATE_TYPE(float) #endif -#define INSTANTIATE_COMP(NUM) \ - template using FS##NUM = GenericOilGasFluidSystem; \ +#define INSTANTIATE_COMP_THREEPHASE(NUM) \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ template class FIPContainer>; -INSTANTIATE_COMP(0) -INSTANTIATE_COMP(2) -INSTANTIATE_COMP(3) -INSTANTIATE_COMP(4) -INSTANTIATE_COMP(5) -INSTANTIATE_COMP(6) -INSTANTIATE_COMP(7) +INSTANTIATE_COMP_THREEPHASE(0) +INSTANTIATE_COMP_THREEPHASE(2) +INSTANTIATE_COMP_THREEPHASE(3) +INSTANTIATE_COMP_THREEPHASE(4) +INSTANTIATE_COMP_THREEPHASE(5) +INSTANTIATE_COMP_THREEPHASE(6) +INSTANTIATE_COMP_THREEPHASE(7) + +#define INSTANTIATE_COMP_TWOPHASE(NUM) \ + template using GFS##NUM = GenericOilGasWaterFluidSystem; \ + template class FIPContainer>; + +INSTANTIATE_COMP_TWOPHASE(2) +INSTANTIATE_COMP_TWOPHASE(3) +INSTANTIATE_COMP_TWOPHASE(4) +INSTANTIATE_COMP_TWOPHASE(5) +INSTANTIATE_COMP_TWOPHASE(6) +INSTANTIATE_COMP_TWOPHASE(7) } // namespace Opm From 21321231a774ea8e918029b435e9320a1f6ed8a1 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Wed, 5 Feb 2025 09:27:02 +0100 Subject: [PATCH 7/9] Simplify macro calls --- opm/simulators/flow/FIPContainer.cpp | 25 ++++++++----------- .../flow/GenericOutputBlackoilModule.cpp | 25 ++++++++----------- 2 files changed, 22 insertions(+), 28 deletions(-) diff --git a/opm/simulators/flow/FIPContainer.cpp b/opm/simulators/flow/FIPContainer.cpp index 00557a1cd..5bc7a9ff5 100644 --- a/opm/simulators/flow/FIPContainer.cpp +++ b/opm/simulators/flow/FIPContainer.cpp @@ -442,23 +442,20 @@ INSTANTIATE_TYPE(float) template using FS##NUM = GenericOilGasWaterFluidSystem; \ template class FIPContainer>; -INSTANTIATE_COMP_THREEPHASE(0) -INSTANTIATE_COMP_THREEPHASE(2) -INSTANTIATE_COMP_THREEPHASE(3) -INSTANTIATE_COMP_THREEPHASE(4) -INSTANTIATE_COMP_THREEPHASE(5) -INSTANTIATE_COMP_THREEPHASE(6) -INSTANTIATE_COMP_THREEPHASE(7) - #define INSTANTIATE_COMP_TWOPHASE(NUM) \ template using GFS##NUM = GenericOilGasWaterFluidSystem; \ template class FIPContainer>; -INSTANTIATE_COMP_TWOPHASE(2) -INSTANTIATE_COMP_TWOPHASE(3) -INSTANTIATE_COMP_TWOPHASE(4) -INSTANTIATE_COMP_TWOPHASE(5) -INSTANTIATE_COMP_TWOPHASE(6) -INSTANTIATE_COMP_TWOPHASE(7) +#define INSTANTIATE_COMP(NUM) \ + INSTANTIATE_COMP_THREEPHASE(NUM) \ + INSTANTIATE_COMP_TWOPHASE(NUM) + +INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput +INSTANTIATE_COMP(2) +INSTANTIATE_COMP(3) +INSTANTIATE_COMP(4) +INSTANTIATE_COMP(5) +INSTANTIATE_COMP(6) +INSTANTIATE_COMP(7) } // namespace Opm diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 19e11d4cd..e882934c8 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -1680,23 +1680,20 @@ INSTANTIATE_TYPE(float) template using FS##NUM = GenericOilGasWaterFluidSystem; \ template class GenericOutputBlackoilModule>; -INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput -INSTANTIATE_COMP_THREEPHASE(2) -INSTANTIATE_COMP_THREEPHASE(3) -INSTANTIATE_COMP_THREEPHASE(4) -INSTANTIATE_COMP_THREEPHASE(5) -INSTANTIATE_COMP_THREEPHASE(6) -INSTANTIATE_COMP_THREEPHASE(7) - #define INSTANTIATE_COMP_TWOPHASE(NUM) \ template using GFS##NUM = GenericOilGasWaterFluidSystem; \ template class GenericOutputBlackoilModule>; -INSTANTIATE_COMP_TWOPHASE(2) -INSTANTIATE_COMP_TWOPHASE(3) -INSTANTIATE_COMP_TWOPHASE(4) -INSTANTIATE_COMP_TWOPHASE(5) -INSTANTIATE_COMP_TWOPHASE(6) -INSTANTIATE_COMP_TWOPHASE(7) +#define INSTANTIATE_COMP(NUM) \ + INSTANTIATE_COMP_THREEPHASE(NUM) \ + INSTANTIATE_COMP_TWOPHASE(NUM) + +INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput +INSTANTIATE_COMP(2) +INSTANTIATE_COMP(3) +INSTANTIATE_COMP(4) +INSTANTIATE_COMP(5) +INSTANTIATE_COMP(6) +INSTANTIATE_COMP(7) } // namespace Opm From c9011448ebe7d2d01676b2de15f20c44edfca8a3 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Fri, 7 Feb 2025 14:08:05 +0100 Subject: [PATCH 8/9] std::move is done inside the fluid system function --- examples/problems/co2ptflashproblem.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/problems/co2ptflashproblem.hh b/examples/problems/co2ptflashproblem.hh index e7df55c89..07eda2d92 100644 --- a/examples/problems/co2ptflashproblem.hh +++ b/examples/problems/co2ptflashproblem.hh @@ -261,7 +261,7 @@ public: ccWaterPvt.setViscosity(/*regionIdx=*/0, 9.6e-4); ccWaterPvt.setCompressibility(/*regionIdx=*/0, 1.450377e-10); waterPvt->initEnd(); - FluidSystem::setWaterPvt(std::move(waterPvt)); + FluidSystem::setWaterPvt(waterPvt); } template From e7a2d954fe4bd0bbc66d31c38e586c6ebb10c6b0 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Thu, 13 Feb 2025 11:08:14 +0100 Subject: [PATCH 9/9] Update Sw in fluid state after normalization --- opm/models/ptflash/flashintensivequantities.hh | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh index 4fad4b706..4237d19cb 100644 --- a/opm/models/ptflash/flashintensivequantities.hh +++ b/opm/models/ptflash/flashintensivequantities.hh @@ -222,7 +222,6 @@ public: Evaluation Sw = 0.0; if constexpr (waterEnabled) { Sw = priVars.makeEvaluation(water0Idx, timeIdx); - fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw); } Evaluation L = fluidState_.L(); Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0); @@ -230,10 +229,13 @@ public: Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg) + Opm::getValue(Sw); So /= sumS; Sg /= sumS; - Sw /= sumS; fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So); fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg); + if constexpr (waterEnabled) { + Sw /= sumS; + fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw); + } fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L); fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V);