From 22dadf5928c397ff45316ce719acbda86f38979d Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 4 Feb 2025 09:43:24 +0100 Subject: [PATCH 01/36] component indices are now int inactive components are represented using negative values in the index(traits) class so we should not be using unsigned values in the fluidsystem while at it make activeToCanonicalComponentIndex usable at compile time --- opm/models/blackoil/blackoilindices.hh | 4 ++-- opm/models/blackoil/blackoilonephaseindices.hh | 8 +++++--- opm/models/blackoil/blackoiltwophaseindices.hh | 8 ++++---- .../wells/MultisegmentWellPrimaryVariables.cpp | 4 ++-- .../wells/MultisegmentWellPrimaryVariables.hpp | 4 ++-- opm/simulators/wells/RatioCalculator.cpp | 6 +++--- opm/simulators/wells/RatioCalculator.hpp | 12 ++++++------ opm/simulators/wells/StandardWellConnections.cpp | 2 +- .../wells/StandardWellPrimaryVariables.cpp | 4 ++-- .../wells/StandardWellPrimaryVariables.hpp | 2 +- opm/simulators/wells/WellInterfaceIndices.cpp | 2 +- opm/simulators/wells/WellInterfaceIndices.hpp | 2 +- 12 files changed, 30 insertions(+), 28 deletions(-) diff --git a/opm/models/blackoil/blackoilindices.hh b/opm/models/blackoil/blackoilindices.hh index 8350373d0..53f4ed6bb 100644 --- a/opm/models/blackoil/blackoilindices.hh +++ b/opm/models/blackoil/blackoilindices.hh @@ -94,10 +94,10 @@ struct BlackOilIndices numEnergy + numFoam + numBrine + numMICPs; //! \brief returns the index of "active" component - static constexpr unsigned canonicalToActiveComponentIndex(unsigned compIdx) + static constexpr int canonicalToActiveComponentIndex(const int compIdx) { return compIdx; } - static constexpr unsigned activeToCanonicalComponentIndex(unsigned compIdx) + static constexpr int activeToCanonicalComponentIndex(const int compIdx) { return compIdx; } //////// diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh index c2abedb02..863c8b00e 100644 --- a/opm/models/blackoil/blackoilonephaseindices.hh +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -30,6 +30,8 @@ #include +#include + namespace Opm { /*! @@ -179,15 +181,15 @@ struct BlackOilOnePhaseIndices ////////////////////// //! \brief returns the index of "active" component - static constexpr unsigned canonicalToActiveComponentIndex(unsigned /*compIdx*/) + static constexpr int canonicalToActiveComponentIndex(const int /*compIdx*/) { return 0; } - static unsigned activeToCanonicalComponentIndex([[maybe_unused]] unsigned compIdx) + static constexpr int activeToCanonicalComponentIndex([[maybe_unused]] const int compIdx) { // assumes canonical oil = 0, water = 1, gas = 2; - assert(compIdx == 0); + constexpr_assert(compIdx == 0); if (gasEnabled) { return 2; } else if (waterEnabled) { diff --git a/opm/models/blackoil/blackoiltwophaseindices.hh b/opm/models/blackoil/blackoiltwophaseindices.hh index e4eefa50d..9d112a148 100644 --- a/opm/models/blackoil/blackoiltwophaseindices.hh +++ b/opm/models/blackoil/blackoiltwophaseindices.hh @@ -181,7 +181,7 @@ struct BlackOilTwoPhaseIndices ////////////////////// //! \brief returns the index of "active" component - static constexpr unsigned canonicalToActiveComponentIndex(const unsigned compIdx) + static constexpr int canonicalToActiveComponentIndex(const int compIdx) { // assumes canonical oil = 0, water = 1, gas = 2; if (!gasEnabled) { @@ -201,10 +201,10 @@ struct BlackOilTwoPhaseIndices return compIdx - 1; } - static unsigned activeToCanonicalComponentIndex(unsigned compIdx) + static constexpr int activeToCanonicalComponentIndex(const int compIdx) { // assumes canonical oil = 0, water = 1, gas = 2; - assert(compIdx < 2); + constexpr_assert(compIdx < 2); if (!gasEnabled) { // oil = 0, water = 1 return compIdx; @@ -212,7 +212,7 @@ struct BlackOilTwoPhaseIndices // oil = 0, gas = 1 return compIdx * 2; } else { - assert(!oilEnabled); + constexpr_assert(!oilEnabled); } // water = 0, gas = 1; diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index 1c2a99507..f3d410b24 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -515,7 +515,7 @@ template typename MultisegmentWellPrimaryVariables::EvalWell MultisegmentWellPrimaryVariables:: volumeFraction(const int seg, - const unsigned compIdx) const + const int compIdx) const { if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { return evaluation_[seg][WFrac]; @@ -584,7 +584,7 @@ typename MultisegmentWellPrimaryVariables::EvalWell MultisegmentWellPrimaryVariables:: getSegmentRateUpwinding(const int seg, const int seg_upwind, - const std::size_t comp_idx) const + const int comp_idx) const { // the result will contain the derivative with respect to WQTotal in segment seg, // and the derivatives with respect to WFrac GFrac in segment seg_upwind. diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp index 3dd05a61c..951d34b4c 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp @@ -115,7 +115,7 @@ public: //! \brief Returns upwinding rate for a component in a segment. EvalWell getSegmentRateUpwinding(const int seg, const int seg_upwind, - const std::size_t comp_idx) const; + const int comp_idx) const; //! \brief Get bottomhole pressure. EvalWell getBhp() const; @@ -154,7 +154,7 @@ private: //! \brief Returns volume fraction for component in a segment. EvalWell volumeFraction(const int seg, - const unsigned compIdx) const; + const int compIdx) const; //! \brief The values for the primary variables //! \details Based on different solution strategies, the wells can have different primary variables diff --git a/opm/simulators/wells/RatioCalculator.cpp b/opm/simulators/wells/RatioCalculator.cpp index d64990b7c..72d8ea5cd 100644 --- a/opm/simulators/wells/RatioCalculator.cpp +++ b/opm/simulators/wells/RatioCalculator.cpp @@ -51,9 +51,9 @@ namespace Opm { template RatioCalculator:: -RatioCalculator(unsigned gasCompIdx, - unsigned oilCompIdx, - unsigned waterCompIdx, +RatioCalculator(int gasCompIdx, + int oilCompIdx, + int waterCompIdx, std::string_view name) : gasComp_{gasCompIdx} , oilComp_(oilCompIdx) diff --git a/opm/simulators/wells/RatioCalculator.hpp b/opm/simulators/wells/RatioCalculator.hpp index 8121c3bb1..985881b88 100644 --- a/opm/simulators/wells/RatioCalculator.hpp +++ b/opm/simulators/wells/RatioCalculator.hpp @@ -39,9 +39,9 @@ class RatioCalculator public: using Scalar = decltype(getValue(Value{})); - RatioCalculator(unsigned gasCompIdx, - unsigned oilCompIdx, - unsigned waterCompIdx, + RatioCalculator(int gasCompIdx, + int oilCompIdx, + int waterCompIdx, std::string_view name); void disOilVapWatVolumeRatio(Value& volumeRatio, @@ -91,9 +91,9 @@ public: const bool isProducer) const; private: - unsigned gasComp_; - unsigned oilComp_; - unsigned waterComp_; + int gasComp_; + int oilComp_; + int waterComp_; std::string name_; }; diff --git a/opm/simulators/wells/StandardWellConnections.cpp b/opm/simulators/wells/StandardWellConnections.cpp index 3b4c17f20..d333ee0fc 100644 --- a/opm/simulators/wells/StandardWellConnections.cpp +++ b/opm/simulators/wells/StandardWellConnections.cpp @@ -731,7 +731,7 @@ connectionRateFoam(const std::vector& cq_s, } case Phase::SOLVENT: { if constexpr (Indices::enableSolvent) - return static_cast(Indices::contiSolventEqIdx); + return Indices::contiSolventEqIdx; else OPM_DEFLOG_THROW(std::runtime_error, "Foam transport phase is SOLVENT but SOLVENT is not activated.", deferred_logger); } diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp index 17410e2f3..3f7369ef9 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.cpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.cpp @@ -446,7 +446,7 @@ copyToWellStatePolyMW(WellState& well_state) const template typename StandardWellPrimaryVariables::EvalWell StandardWellPrimaryVariables:: -volumeFraction(const unsigned compIdx) const +volumeFraction(const int compIdx) const { if (FluidSystem::numActivePhases() == 1) { return EvalWell(numWellEq_ + Indices::numEq, 1.0); @@ -456,7 +456,7 @@ volumeFraction(const unsigned compIdx) const return evaluation_[GFrac]; } - if (Indices::enableSolvent && compIdx == (unsigned)Indices::contiSolventEqIdx) { + if (Indices::enableSolvent && compIdx == Indices::contiSolventEqIdx) { return evaluation_[SFrac]; } diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index f482af12e..7c0c551ac 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp @@ -157,7 +157,7 @@ private: DeferredLogger& deferred_logger) const; //! \brief Returns volume fraction for a component. - EvalWell volumeFraction(const unsigned compIdx) const; + EvalWell volumeFraction(const int compIdx) const; //! \brief Handle non-reasonable fractions due to numerical overshoot. void processFractions(); diff --git a/opm/simulators/wells/WellInterfaceIndices.cpp b/opm/simulators/wells/WellInterfaceIndices.cpp index 1150bdc18..983318488 100644 --- a/opm/simulators/wells/WellInterfaceIndices.cpp +++ b/opm/simulators/wells/WellInterfaceIndices.cpp @@ -78,7 +78,7 @@ flowPhaseToModelCompIdx(const int phaseIdx) const template int WellInterfaceIndices:: -modelCompIdxToFlowCompIdx(const unsigned compIdx) const +modelCompIdxToFlowCompIdx(const int compIdx) const { const auto& pu = this->phaseUsage(); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == compIdx) diff --git a/opm/simulators/wells/WellInterfaceIndices.hpp b/opm/simulators/wells/WellInterfaceIndices.hpp index 6ac6138e1..16146c9a9 100644 --- a/opm/simulators/wells/WellInterfaceIndices.hpp +++ b/opm/simulators/wells/WellInterfaceIndices.hpp @@ -41,7 +41,7 @@ public: using ModelParameters = typename WellInterfaceFluidSystem::ModelParameters; int flowPhaseToModelCompIdx(const int phaseIdx) const; - int modelCompIdxToFlowCompIdx(const unsigned compIdx) const; + int modelCompIdxToFlowCompIdx(const int compIdx) const; Scalar scalingFactor(const int phaseIdx) const; template From df6bdbe70087d3d849215432d33a7b346a0416f5 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Wed, 8 Jan 2025 13:08:17 +0100 Subject: [PATCH 02/36] 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 03/36] 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 04/36] 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 05/36] 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 06/36] 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 07/36] 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 08/36] 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 09/36] 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 8b84bc6b4b6861987b41425e550e72f9a6511453 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Sat, 8 Feb 2025 12:50:56 +0100 Subject: [PATCH 10/36] Using the CriticalException in PreconditionerFactory and ISTLSolver --- opm/simulators/linalg/ISTLSolver.hpp | 8 +++++--- opm/simulators/linalg/PreconditionerFactory_impl.hpp | 7 ++++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index 053a6fbc6..93ef60c75 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -377,10 +378,11 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, void prepare(const Matrix& M, Vector& b) { OPM_TIMEBLOCK(istlSolverPrepare); + OPM_BEGIN_TRY_CATCH_RETHROW_AS_CRITICAL_ERROR(); + initPrepare(M,b); - initPrepare(M,b); - - prepareFlexibleSolver(); + prepareFlexibleSolver(); + OPM_END_TRY_CATCH_RETHROW_AS_CRITICAL_ERROR(); } diff --git a/opm/simulators/linalg/PreconditionerFactory_impl.hpp b/opm/simulators/linalg/PreconditionerFactory_impl.hpp index 44287b6eb..b6dca5f0a 100644 --- a/opm/simulators/linalg/PreconditionerFactory_impl.hpp +++ b/opm/simulators/linalg/PreconditionerFactory_impl.hpp @@ -20,6 +20,7 @@ #include +#include #include #include @@ -776,7 +777,7 @@ PreconditionerFactory::create(const Operator& op, const std::function& weightsCalculator, std::size_t pressureIndex) { - return instance().doCreate(op, prm, weightsCalculator, pressureIndex); + OPM_TRY_THROW_AS_CRITICAL_ERROR(return instance().doCreate(op, prm, weightsCalculator, pressureIndex)); } template @@ -787,7 +788,7 @@ PreconditionerFactory::create(const Operator& op, const Comm& comm, std::size_t pressureIndex) { - return instance().doCreate(op, prm, weightsCalculator, pressureIndex, comm); + OPM_TRY_THROW_AS_CRITICAL_ERROR(return instance().doCreate(op, prm, weightsCalculator, pressureIndex, comm)); } @@ -798,7 +799,7 @@ PreconditionerFactory::create(const Operator& op, const Comm& comm, std::size_t pressureIndex) { - return instance().doCreate(op, prm, std::function(), pressureIndex, comm); + OPM_TRY_THROW_AS_CRITICAL_ERROR(return instance().doCreate(op, prm, std::function(), pressureIndex, comm)); } template From 264a8d44d39c63c5ee8e237b29cab28c96cf2bc7 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Sat, 8 Feb 2025 19:39:15 +0100 Subject: [PATCH 11/36] Removed rethrow from PreconditionerFactory::create --- opm/simulators/linalg/PreconditionerFactory_impl.hpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/opm/simulators/linalg/PreconditionerFactory_impl.hpp b/opm/simulators/linalg/PreconditionerFactory_impl.hpp index b6dca5f0a..44287b6eb 100644 --- a/opm/simulators/linalg/PreconditionerFactory_impl.hpp +++ b/opm/simulators/linalg/PreconditionerFactory_impl.hpp @@ -20,7 +20,6 @@ #include -#include #include #include @@ -777,7 +776,7 @@ PreconditionerFactory::create(const Operator& op, const std::function& weightsCalculator, std::size_t pressureIndex) { - OPM_TRY_THROW_AS_CRITICAL_ERROR(return instance().doCreate(op, prm, weightsCalculator, pressureIndex)); + return instance().doCreate(op, prm, weightsCalculator, pressureIndex); } template @@ -788,7 +787,7 @@ PreconditionerFactory::create(const Operator& op, const Comm& comm, std::size_t pressureIndex) { - OPM_TRY_THROW_AS_CRITICAL_ERROR(return instance().doCreate(op, prm, weightsCalculator, pressureIndex, comm)); + return instance().doCreate(op, prm, weightsCalculator, pressureIndex, comm); } @@ -799,7 +798,7 @@ PreconditionerFactory::create(const Operator& op, const Comm& comm, std::size_t pressureIndex) { - OPM_TRY_THROW_AS_CRITICAL_ERROR(return instance().doCreate(op, prm, std::function(), pressureIndex, comm)); + return instance().doCreate(op, prm, std::function(), pressureIndex, comm); } template From 791e5974f240c283b3ec03ec5788c16817f4391a Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Mon, 10 Feb 2025 16:08:43 +0100 Subject: [PATCH 12/36] Now providing a message for CritcalError when rethrowing in ISTLSolver --- opm/simulators/linalg/ISTLSolver.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index 93ef60c75..e118fb647 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -378,11 +378,11 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, void prepare(const Matrix& M, Vector& b) { OPM_TIMEBLOCK(istlSolverPrepare); - OPM_BEGIN_TRY_CATCH_RETHROW_AS_CRITICAL_ERROR(); + try { initPrepare(M,b); prepareFlexibleSolver(); - OPM_END_TRY_CATCH_RETHROW_AS_CRITICAL_ERROR(); + } OPM_CATCH_AND_RETHROW_AS_CRITICAL_ERROR("This is likely due to a faulty linear solver JSON specification. Check for errors related to missing nodes."); } From 66fa971079b51a2ddba5816f8b7c67bfdc3b175a Mon Sep 17 00:00:00 2001 From: Tobias Meyer Andersen Date: Thu, 5 Dec 2024 10:47:05 +0100 Subject: [PATCH 13/36] use cudagraphs in gpu DILU and ILU0 apply expected to give performance gain on current nvidia cards, no significant change on AMD cards for cases of reasonable size uses new GPU convenence objects --- opm/simulators/linalg/gpuistl/GpuBuffer.cpp | 1 + opm/simulators/linalg/gpuistl/GpuBuffer.hpp | 1 + opm/simulators/linalg/gpuistl/GpuDILU.cpp | 75 +++++- opm/simulators/linalg/gpuistl/GpuDILU.hpp | 19 +- opm/simulators/linalg/gpuistl/GpuVector.cpp | 13 + opm/simulators/linalg/gpuistl/GpuVector.hpp | 1 + opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp | 77 ++++-- opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp | 16 +- .../preconditionerKernels/DILUKernels.cu | 245 +++++++++++++----- .../preconditionerKernels/DILUKernels.hpp | 12 +- .../preconditionerKernels/ILU0Kernels.cu | 94 +++---- .../preconditionerKernels/ILU0Kernels.hpp | 12 +- 12 files changed, 417 insertions(+), 149 deletions(-) diff --git a/opm/simulators/linalg/gpuistl/GpuBuffer.cpp b/opm/simulators/linalg/gpuistl/GpuBuffer.cpp index b4132ce09..49ce121e1 100644 --- a/opm/simulators/linalg/gpuistl/GpuBuffer.cpp +++ b/opm/simulators/linalg/gpuistl/GpuBuffer.cpp @@ -188,6 +188,7 @@ GpuBuffer::copyFromHost(const std::vector& data) { copyFromHost(data.data(), data.size()); } + template void GpuBuffer::copyToHost(std::vector& data) const diff --git a/opm/simulators/linalg/gpuistl/GpuBuffer.hpp b/opm/simulators/linalg/gpuistl/GpuBuffer.hpp index 369ea3ba2..1ba5cb861 100644 --- a/opm/simulators/linalg/gpuistl/GpuBuffer.hpp +++ b/opm/simulators/linalg/gpuistl/GpuBuffer.hpp @@ -27,6 +27,7 @@ #include #include #include +#include namespace Opm::gpuistl diff --git a/opm/simulators/linalg/gpuistl/GpuDILU.cpp b/opm/simulators/linalg/gpuistl/GpuDILU.cpp index d311386b0..2011fb2a2 100644 --- a/opm/simulators/linalg/gpuistl/GpuDILU.cpp +++ b/opm/simulators/linalg/gpuistl/GpuDILU.cpp @@ -37,6 +37,7 @@ #include #include #include + namespace Opm::gpuistl { @@ -93,7 +94,8 @@ GpuDILU::GpuDILU(const M& A, bool splitMatrix, bool tuneKernels, int } } - computeDiagAndMoveReorderedData(m_moveThreadBlockSize, m_DILUFactorizationThreadBlockSize); + reorderAndSplitMatrix(m_moveThreadBlockSize); + computeDiagonal(m_DILUFactorizationThreadBlockSize); if (m_tuneThreadBlockSizes) { tuneThreadBlockSizes(); @@ -110,10 +112,31 @@ template void GpuDILU::apply(X& v, const Y& d) { + // ensure that this stream only starts doing work when main stream is completed up to this point + OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0)); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0)); + OPM_TIMEBLOCK(prec_apply); { - apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); + const auto ptrs = std::make_pair(v.data(), d.data()); + + auto it = m_apply_graphs.find(ptrs); + + if (it == m_apply_graphs.end()) { + OPM_GPU_SAFE_CALL(cudaStreamBeginCapture(m_stream.get(), cudaStreamCaptureModeGlobal)); + + // The apply functions contains lots of small function calls which call a kernel each + apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); + + OPM_GPU_SAFE_CALL(cudaStreamEndCapture(m_stream.get(), &m_apply_graphs[ptrs].get())); + OPM_GPU_SAFE_CALL(cudaGraphInstantiate(&m_executableGraphs[ptrs].get(), m_apply_graphs[ptrs].get(), nullptr, nullptr, 0)); + } + OPM_GPU_SAFE_CALL(cudaGraphLaunch(m_executableGraphs[ptrs].get(), 0)); } + + // ensure that main stream only continues after this stream is completed + OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get())); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0)); } template @@ -135,7 +158,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int m_gpuDInvFloat->data(), d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) { detail::DILU::solveLowerLevelSetSplit( m_gpuMatrixReorderedLowerFloat->getNonZeroValues().data(), @@ -147,7 +171,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int m_gpuDInv.data(), d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } else { detail::DILU::solveLowerLevelSetSplit( m_gpuMatrixReorderedLower->getNonZeroValues().data(), @@ -159,7 +184,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int m_gpuDInv.data(), d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } } else { detail::DILU::solveLowerLevelSet( @@ -172,7 +198,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int m_gpuDInv.data(), d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } levelStartIdx += numOfRowsInLevel; } @@ -193,7 +220,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int numOfRowsInLevel, m_gpuDInvFloat->data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG){ detail::DILU::solveUpperLevelSetSplit( m_gpuMatrixReorderedUpperFloat->getNonZeroValues().data(), @@ -204,7 +232,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int numOfRowsInLevel, m_gpuDInv.data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } else { detail::DILU::solveUpperLevelSetSplit( m_gpuMatrixReorderedUpper->getNonZeroValues().data(), @@ -215,7 +244,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int numOfRowsInLevel, m_gpuDInv.data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } } else { detail::DILU::solveUpperLevelSet( @@ -227,7 +257,8 @@ GpuDILU::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, int numOfRowsInLevel, m_gpuDInv.data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } } } @@ -252,7 +283,9 @@ GpuDILU::update() { OPM_TIMEBLOCK(prec_update); { - update(m_moveThreadBlockSize, m_DILUFactorizationThreadBlockSize); + m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu + reorderAndSplitMatrix(m_moveThreadBlockSize); + computeDiagonal(m_DILUFactorizationThreadBlockSize); } } @@ -260,13 +293,22 @@ template void GpuDILU::update(int moveThreadBlockSize, int factorizationBlockSize) { - m_gpuMatrix.updateNonzeroValues(m_cpuMatrix, true); // send updated matrix to the gpu - computeDiagAndMoveReorderedData(moveThreadBlockSize, factorizationBlockSize); + // ensure that this stream only starts doing work when main stream is completed up to this point + OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0)); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0)); + + m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu + reorderAndSplitMatrix(moveThreadBlockSize); + computeDiagonal(factorizationBlockSize); + + // ensure that main stream only continues after this stream is completed + OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get())); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0)); } template void -GpuDILU::computeDiagAndMoveReorderedData(int moveThreadBlockSize, int factorizationBlockSize) +GpuDILU::reorderAndSplitMatrix(int moveThreadBlockSize) { if (m_splitMatrix) { detail::copyMatDataToReorderedSplit( @@ -290,7 +332,12 @@ GpuDILU::computeDiagAndMoveReorderedData(int moveThreadBlockSize, in m_gpuMatrixReordered->N(), moveThreadBlockSize); } +} +template +void +GpuDILU::computeDiagonal(int factorizationBlockSize) +{ int levelStartIdx = 0; for (int level = 0; level < m_levelSets.size(); ++level) { const int numOfRowsInLevel = m_levelSets[level].size(); diff --git a/opm/simulators/linalg/gpuistl/GpuDILU.hpp b/opm/simulators/linalg/gpuistl/GpuDILU.hpp index 241eef3e9..a023c05d4 100644 --- a/opm/simulators/linalg/gpuistl/GpuDILU.hpp +++ b/opm/simulators/linalg/gpuistl/GpuDILU.hpp @@ -24,8 +24,10 @@ #include #include #include +#include #include - +#include +#include namespace Opm::gpuistl @@ -82,8 +84,11 @@ public: //! \brief Updates the matrix data. void update() final; + //! \brief perform matrix splitting and reordering + void reorderAndSplitMatrix(int moveThreadBlockSize); + //! \brief Compute the diagonal of the DILU, and update the data of the reordered matrix - void computeDiagAndMoveReorderedData(int moveThreadBlockSize, int factorizationThreadBlockSize); + void computeDiagonal(int factorizationThreadBlockSize); //! \brief function that will experimentally tune the thread block sizes of the important cuda kernels void tuneThreadBlockSizes(); @@ -153,6 +158,16 @@ private: int m_lowerSolveThreadBlockSize = -1; int m_moveThreadBlockSize = -1; int m_DILUFactorizationThreadBlockSize = -1; + + // Graphs for Apply + std::map, GPUGraph> m_apply_graphs; + std::map, GPUGraphExec> m_executableGraphs; + + // Stream for the DILU operations on the GPU + GPUStream m_stream{}; + // Events for synchronization with main stream + GPUEvent m_before{}; + GPUEvent m_after{}; }; } // end namespace Opm::gpuistl diff --git a/opm/simulators/linalg/gpuistl/GpuVector.cpp b/opm/simulators/linalg/gpuistl/GpuVector.cpp index ddafc1dbf..af39e9941 100644 --- a/opm/simulators/linalg/gpuistl/GpuVector.cpp +++ b/opm/simulators/linalg/gpuistl/GpuVector.cpp @@ -266,6 +266,19 @@ GpuVector::copyFromHost(const T* dataPointer, size_t numberOfElements) OPM_GPU_SAFE_CALL(cudaMemcpy(data(), dataPointer, numberOfElements * sizeof(T), cudaMemcpyHostToDevice)); } +template +void +GpuVector::copyFromHost(const T* dataPointer, size_t numberOfElements, cudaStream_t stream) +{ + if (numberOfElements > dim()) { + OPM_THROW(std::runtime_error, + fmt::format("Requesting to copy too many elements. Vector has {} elements, while {} was requested.", + dim(), + numberOfElements)); + } + OPM_GPU_SAFE_CALL(cudaMemcpyAsync(data(), dataPointer, numberOfElements * sizeof(T), cudaMemcpyHostToDevice, stream)); +} + template void GpuVector::copyToHost(T* dataPointer, size_t numberOfElements) const diff --git a/opm/simulators/linalg/gpuistl/GpuVector.hpp b/opm/simulators/linalg/gpuistl/GpuVector.hpp index e369d1d9d..66a1ebd65 100644 --- a/opm/simulators/linalg/gpuistl/GpuVector.hpp +++ b/opm/simulators/linalg/gpuistl/GpuVector.hpp @@ -203,6 +203,7 @@ public: * @note assumes that this vector has numberOfElements elements */ void copyFromHost(const T* dataPointer, size_t numberOfElements); + void copyFromHost(const T* dataPointer, size_t numberOfElements, cudaStream_t stream); /** * @brief copyFromHost copies numberOfElements to the CPU memory dataPointer diff --git a/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp b/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp index 686b8f09f..42b535b76 100644 --- a/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp +++ b/opm/simulators/linalg/gpuistl/OpmGpuILU0.cpp @@ -37,6 +37,7 @@ #include #include #include + namespace Opm::gpuistl { @@ -58,7 +59,6 @@ OpmGpuILU0::OpmGpuILU0(const M& A, bool splitMatrix, bool tuneKernel , m_tuneThreadBlockSizes(tuneKernels) , m_mixedPrecisionScheme(makeMatrixStorageMPScheme(mixedPrecisionScheme)) { - // TODO: Should in some way verify that this matrix is symmetric, only do it debug mode? // Some sanity check OPM_ERROR_IF(A.N() != m_gpuMatrix.N(), @@ -98,7 +98,8 @@ OpmGpuILU0::OpmGpuILU0(const M& A, bool splitMatrix, bool tuneKernel } } - LUFactorizeAndMoveData(m_moveThreadBlockSize, m_ILU0FactorizationThreadBlockSize); + reorderAndSplitMatrix(m_moveThreadBlockSize); + LUFactorizeMatrix(m_ILU0FactorizationThreadBlockSize); if (m_tuneThreadBlockSizes) { tuneThreadBlockSizes(); @@ -117,7 +118,29 @@ OpmGpuILU0::apply(X& v, const Y& d) { OPM_TIMEBLOCK(prec_apply); { - apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); + // ensure that this stream only starts doing work when main stream is completed up to this point + OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0)); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0)); + + const auto ptrs = std::make_pair(v.data(), d.data()); + + auto it = m_apply_graphs.find(ptrs); + + if (it == m_apply_graphs.end()) { + OPM_GPU_SAFE_CALL(cudaStreamBeginCapture(m_stream.get(), cudaStreamCaptureModeGlobal)); + + // The apply functions contains lots of small function calls which call a kernel each + apply(v, d, m_lowerSolveThreadBlockSize, m_upperSolveThreadBlockSize); + + OPM_GPU_SAFE_CALL(cudaStreamEndCapture(m_stream.get(), &m_apply_graphs[ptrs].get())); + OPM_GPU_SAFE_CALL(cudaGraphInstantiate(&m_executableGraphs[ptrs].get(), m_apply_graphs[ptrs].get(), nullptr, nullptr, 0)); + } + OPM_GPU_SAFE_CALL(cudaGraphLaunch(m_executableGraphs[ptrs].get(), 0)); + + + // ensure that main stream only continues after this stream is completed + OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get())); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0)); } } @@ -142,7 +165,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } else{ detail::ILU0::solveLowerLevelSetSplit( @@ -154,7 +178,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } } else { detail::ILU0::solveLowerLevelSet(m_gpuReorderedLU->getNonZeroValues().data(), @@ -165,7 +190,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, d.data(), v.data(), - lowerSolveThreadBlockSize); + lowerSolveThreadBlockSize, + m_stream.get()); } levelStartIdx += numOfRowsInLevel; } @@ -185,7 +211,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, m_gpuMatrixReorderedDiagFloat.value().data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } else if (m_mixedPrecisionScheme == MatrixStorageMPScheme::DOUBLE_DIAG_FLOAT_OFFDIAG) { detail::ILU0::solveUpperLevelSetSplit( @@ -197,7 +224,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, m_gpuMatrixReorderedDiag.value().data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } else{ detail::ILU0::solveUpperLevelSetSplit( @@ -209,7 +237,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i numOfRowsInLevel, m_gpuMatrixReorderedDiag.value().data(), v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } } else { detail::ILU0::solveUpperLevelSet(m_gpuReorderedLU->getNonZeroValues().data(), @@ -219,7 +248,8 @@ OpmGpuILU0::apply(X& v, const Y& d, int lowerSolveThreadBlockSize, i levelStartIdx, numOfRowsInLevel, v.data(), - upperSolveThreadBlockSize); + upperSolveThreadBlockSize, + m_stream.get()); } } } @@ -241,10 +271,9 @@ template void OpmGpuILU0::update() { - OPM_TIMEBLOCK(prec_update); - { - update(m_moveThreadBlockSize, m_ILU0FactorizationThreadBlockSize); - } + m_gpuMatrix.updateNonzeroValues(m_cpuMatrix); // send updated matrix to the gpu + reorderAndSplitMatrix(m_moveThreadBlockSize); + LUFactorizeMatrix(m_ILU0FactorizationThreadBlockSize); } template @@ -253,13 +282,23 @@ OpmGpuILU0::update(int moveThreadBlockSize, int factorizationThreadB { OPM_TIMEBLOCK(prec_update); { + // ensure that this stream only starts doing work when main stream is completed up to this point + OPM_GPU_SAFE_CALL(cudaEventRecord(m_before.get(), 0)); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(m_stream.get(), m_before.get(), 0)); + m_gpuMatrix.updateNonzeroValues(m_cpuMatrix, true); // send updated matrix to the gpu - LUFactorizeAndMoveData(moveThreadBlockSize, factorizationThreadBlockSize); + reorderAndSplitMatrix(moveThreadBlockSize); + LUFactorizeMatrix(factorizationThreadBlockSize); + + // ensure that main stream only continues after this stream is completed + OPM_GPU_SAFE_CALL(cudaEventRecord(m_after.get(), m_stream.get())); + OPM_GPU_SAFE_CALL(cudaStreamWaitEvent(0, m_after.get(), 0)); } } + template void -OpmGpuILU0::LUFactorizeAndMoveData(int moveThreadBlockSize, int factorizationThreadBlockSize) +OpmGpuILU0::reorderAndSplitMatrix(int moveThreadBlockSize) { if (m_splitMatrix) { detail::copyMatDataToReorderedSplit( @@ -283,6 +322,12 @@ OpmGpuILU0::LUFactorizeAndMoveData(int moveThreadBlockSize, int fact m_gpuReorderedLU->N(), moveThreadBlockSize); } +} + +template +void +OpmGpuILU0::LUFactorizeMatrix(int factorizationThreadBlockSize) +{ int levelStartIdx = 0; for (int level = 0; level < m_levelSets.size(); ++level) { const int numOfRowsInLevel = m_levelSets[level].size(); diff --git a/opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp b/opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp index c163e0483..e4cc3f059 100644 --- a/opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp +++ b/opm/simulators/linalg/gpuistl/OpmGpuILU0.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -84,8 +85,11 @@ public: //! \brief Updates the matrix data. void update() final; + //! \brief perform matrix splitting and reordering + void reorderAndSplitMatrix(int moveThreadBlockSize); + //! \brief Compute LU factorization, and update the data of the reordered matrix - void LUFactorizeAndMoveData(int moveThreadBlockSize, int factorizationThreadBlockSize); + void LUFactorizeMatrix(int factorizationThreadBlockSize); //! \brief function that will experimentally tune the thread block sizes of the important cuda kernels void tuneThreadBlockSizes(); @@ -152,6 +156,16 @@ private: int m_lowerSolveThreadBlockSize = -1; int m_moveThreadBlockSize = -1; int m_ILU0FactorizationThreadBlockSize = -1; + + // Graphs for Apply + std::map, GPUGraph> m_apply_graphs; + std::map, GPUGraphExec> m_executableGraphs; + + // Stream for the DILU operations on the GPU + GPUStream m_stream{}; + // Events for synchronization with main stream + GPUEvent m_before{}; + GPUEvent m_after{}; }; } // end namespace Opm::gpuistl diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu index b8f9f48ca..8af6765df 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.cu @@ -85,10 +85,12 @@ namespace // TODO: removce the first condition in the for loop for (int block = nnzIdx; block < nnzIdxLim; ++block) { const int col = colIndices[block]; - mmvMixedGeneral(&mat[block * blocksize * blocksize], &v[col * blocksize], rhs); + mmvMixedGeneral( + &mat[block * blocksize * blocksize], &v[col * blocksize], rhs); } - mvMixedGeneral(&dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]); + mvMixedGeneral( + &dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]); } } @@ -137,10 +139,12 @@ namespace LinearSolverScalar rhs[blocksize] = {0}; for (int block = nnzIdx; block < nnzIdxLim; ++block) { const int col = colIndices[block]; - umvMixedGeneral(&mat[block * blocksize * blocksize], &v[col * blocksize], rhs); + umvMixedGeneral( + &mat[block * blocksize * blocksize], &v[col * blocksize], rhs); } - mmvMixedGeneral(&dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]); + mmvMixedGeneral( + &dInv[reorderedRowIdx * blocksize * blocksize], rhs, &v[naturalRowIdx * blocksize]); } } @@ -211,8 +215,8 @@ namespace } } - // TODO: rewrite such that during the factorization there is a dInv of InputScalar type that stores intermediate results - // TOOD: The important part is to only cast after that is fully computed + // TODO: rewrite such that during the factorization there is a dInv of InputScalar type that stores intermediate + // results TOOD: The important part is to only cast after that is fully computed template __global__ void cuComputeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat, int* lowerRowIndices, @@ -239,7 +243,8 @@ namespace InputScalar dInvTmp[blocksize * blocksize]; for (int i = 0; i < blocksize; ++i) { for (int j = 0; j < blocksize; ++j) { - dInvTmp[i * blocksize + j] = srcDiagonal[reorderedRowIdx * blocksize * blocksize + i * blocksize + j]; + dInvTmp[i * blocksize + j] + = srcDiagonal[reorderedRowIdx * blocksize * blocksize + i * blocksize + j]; } } @@ -257,21 +262,26 @@ namespace if constexpr (detail::storeOffDiagonalAsFloat(mixedPrecisionScheme)) { // TODO: think long and hard about whether this performs only the wanted memory transfers - moveBlock(&srcReorderedLowerMat[block * blocksize * blocksize], &dstLowerMat[block * blocksize * blocksize]); - moveBlock(&srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], &dstUpperMat[symOppositeBlock * blocksize * blocksize]); + moveBlock( + &srcReorderedLowerMat[block * blocksize * blocksize], + &dstLowerMat[block * blocksize * blocksize]); + moveBlock( + &srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], + &dstUpperMat[symOppositeBlock * blocksize * blocksize]); } mmx2Subtraction(&srcReorderedLowerMat[block * blocksize * blocksize], - &dInv[col * blocksize * blocksize], - &srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], - dInvTmp); + &dInv[col * blocksize * blocksize], + &srcReorderedUpperMat[symOppositeBlock * blocksize * blocksize], + dInvTmp); } invBlockInPlace(dInvTmp); moveBlock(dInvTmp, &dInv[reorderedRowIdx * blocksize * blocksize]); if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)) { - moveBlock(dInvTmp, &dstDiag[reorderedRowIdx * blocksize * blocksize]); // important! + moveBlock( + dInvTmp, &dstDiag[reorderedRowIdx * blocksize * blocksize]); // important! } } } @@ -289,12 +299,13 @@ solveLowerLevelSet(T* reorderedMat, const T* dInv, const T* d, T* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveLowerLevelSet, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveLowerLevelSet<<>>( + cuSolveLowerLevelSet<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v); } @@ -310,13 +321,15 @@ solveLowerLevelSetSplit(MatrixScalar* reorderedMat, const DiagonalScalar* dInv, const LinearSolverScalar* d, LinearSolverScalar* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuSolveLowerLevelSetSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveLowerLevelSetSplit<<>>( - reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v); + cuSolveLowerLevelSetSplit + <<>>( + reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, d, v); } // perform the upper solve for all rows in the same level set template @@ -329,12 +342,13 @@ solveUpperLevelSet(T* reorderedMat, int rowsInLevelSet, const T* dInv, T* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveUpperLevelSet, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveUpperLevelSet<<>>( + cuSolveUpperLevelSet<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); } @@ -348,13 +362,15 @@ solveUpperLevelSetSplit(MatrixScalar* reorderedMat, int rowsInLevelSet, const DiagonalScalar* dInv, LinearSolverScalar* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuSolveUpperLevelSetSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveUpperLevelSetSplit<<>>( - reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); + cuSolveUpperLevelSetSplit + <<>>( + reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); } template @@ -409,43 +425,134 @@ computeDiluDiagonalSplit(const InputScalar* srcReorderedLowerMat, int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuComputeDiluDiagonalSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuComputeDiluDiagonalSplit<<>>(srcReorderedLowerMat, - lowerRowIndices, - lowerColIndices, - srcReorderedUpperMat, - upperRowIndices, - upperColIndices, - srcDiagonal, - reorderedToNatural, - naturalToReordered, - startIdx, - rowsInLevelSet, - dInv, - dstDiag, - dstLowerMat, - dstUpperMat); + cuComputeDiluDiagonalSplit + <<>>(srcReorderedLowerMat, + lowerRowIndices, + lowerColIndices, + srcReorderedUpperMat, + upperRowIndices, + upperColIndices, + srcDiagonal, + reorderedToNatural, + naturalToReordered, + startIdx, + rowsInLevelSet, + dInv, + dstDiag, + dstLowerMat, + dstUpperMat); } else { OPM_THROW(std::invalid_argument, "Inverting diagonal is not implemented for blocksizes > 3"); } } -// TODO: format #define INSTANTIATE_KERNEL_WRAPPERS(T, blocksize) \ template void computeDiluDiagonal(T*, int*, int*, int*, int*, const int, int, T*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, float*, float*, float*, int); \ - template void computeDiluDiagonalSplit( \ - const T*, int*, int*, const T*, int*, int*, const T*, int*, int*, const int, int, T*, double*, double*, double*, int); \ - template void solveUpperLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int); \ - template void solveLowerLevelSet(T*, int*, int*, int*, int, int, const T*, const T*, T*, int); + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + double*, \ + double*, \ + double*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + float*, \ + float*, \ + float*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + float*, \ + float*, \ + float*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + double*, \ + double*, \ + double*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + float*, \ + float*, \ + float*, \ + int); \ + template void computeDiluDiagonalSplit( \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const T*, \ + int*, \ + int*, \ + const int, \ + int, \ + T*, \ + double*, \ + double*, \ + double*, \ + int); \ + template void solveUpperLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int, cudaStream_t); \ + template void solveLowerLevelSet( \ + T*, int*, int*, int*, int, int, const T*, const T*, T*, int, cudaStream_t); INSTANTIATE_KERNEL_WRAPPERS(float, 1); INSTANTIATE_KERNEL_WRAPPERS(float, 2); @@ -460,19 +567,29 @@ INSTANTIATE_KERNEL_WRAPPERS(double, 4); INSTANTIATE_KERNEL_WRAPPERS(double, 5); INSTANTIATE_KERNEL_WRAPPERS(double, 6); -#define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar) \ - template void solveUpperLevelSetSplit( \ - MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, LinearSolverScalar*, int); \ - template void solveLowerLevelSetSplit( \ - MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, const LinearSolverScalar*, LinearSolverScalar*, int); +#define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, LinearSolverScalar, MatrixScalar, DiagonalScalar) \ + template void solveUpperLevelSetSplit( \ + MatrixScalar*, int*, int*, int*, int, int, const DiagonalScalar*, LinearSolverScalar*, int, cudaStream_t); \ + template void solveLowerLevelSetSplit( \ + MatrixScalar*, \ + int*, \ + int*, \ + int*, \ + int, \ + int, \ + const DiagonalScalar*, \ + const LinearSolverScalar*, \ + LinearSolverScalar*, \ + int, \ + cudaStream_t); // TODO: be smarter about this... Surely this instantiates many more combinations that are actually needed -#define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT_ALL(blocksize) \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, float, float, float); \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, double, float); \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, float, float); \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, float, float, double); \ - INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, double, double); \ +#define INSTANTIATE_SOLVE_LEVEL_SET_SPLIT_ALL(blocksize) \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, float, float, float); \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, double, float); \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, float, float); \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, float, float, double); \ + INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, double, double); \ INSTANTIATE_SOLVE_LEVEL_SET_SPLIT(blocksize, double, float, double); INSTANTIATE_SOLVE_LEVEL_SET_SPLIT_ALL(1); diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp index b7141834f..77a9cae69 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/DILUKernels.hpp @@ -54,7 +54,8 @@ void solveLowerLevelSet(T* reorderedMat, const T* dInv, const T* d, T* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform a lower solve on certain rows in a matrix that can safely be computed in parallel @@ -82,7 +83,8 @@ void solveLowerLevelSetSplit(MatrixScalar* reorderedUpperMat, const DiagonalScalar* dInv, const LinearSolverScalar* d, LinearSolverScalar* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel @@ -108,7 +110,8 @@ void solveUpperLevelSet(T* reorderedMat, int rowsInLevelSet, const T* dInv, T* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel @@ -134,7 +137,8 @@ void solveUpperLevelSetSplit(MatrixScalar* reorderedUpperMat, int rowsInLevelSet, const DiagonalScalar* dInv, LinearSolverScalar* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Computes the ILU0 of the diagonal elements of the reordered matrix and stores it in a reordered vector diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu index 46183efd2..e66c42f8a 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.cu @@ -215,9 +215,9 @@ namespace // as of now, if we are using mixed precision, then we are always storing the off-diagonals as floats, // and sometimes also the diagonal. if constexpr (detail::usingMixedPrecision(mixedPrecisionScheme)) { - // if we are want to store the entire matrix as a float then we must also move the diagonal block from double to float - // if not then we just use the double diagonal that is already now stored in srcDiagonal - if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)){ + // if we are want to store the entire matrix as a float then we must also move the diagonal block from + // double to float if not then we just use the double diagonal that is already now stored in srcDiagonal + if constexpr (detail::storeDiagonalAsFloat(mixedPrecisionScheme)) { moveBlock(&srcDiagonal[reorderedIdx * scalarsInBlock], &dstDiagonal[reorderedIdx * scalarsInBlock]); } @@ -362,12 +362,13 @@ solveLowerLevelSet(T* reorderedMat, int rowsInLevelSet, const T* d, T* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveLowerLevelSet, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveLowerLevelSet<<>>( + cuSolveLowerLevelSet<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, d, v); } // perform the upper solve for all rows in the same level set @@ -380,12 +381,13 @@ solveUpperLevelSet(T* reorderedMat, int startIdx, int rowsInLevelSet, T* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize(cuSolveUpperLevelSet, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveUpperLevelSet<<>>( + cuSolveUpperLevelSet<<>>( reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, v); } @@ -399,13 +401,15 @@ solveLowerLevelSetSplit(MatrixScalar* reorderedMat, int rowsInLevelSet, const LinearSolverScalar* d, LinearSolverScalar* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuSolveLowerLevelSetSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveLowerLevelSetSplit<<>>( - reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, d, v); + cuSolveLowerLevelSetSplit + <<>>( + reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, d, v); } // perform the upper solve for all rows in the same level set template @@ -418,13 +422,15 @@ solveUpperLevelSetSplit(MatrixScalar* reorderedMat, int rowsInLevelSet, const DiagonalScalar* dInv, LinearSolverScalar* v, - int thrBlockSize) + int thrBlockSize, + cudaStream_t stream) { int threadBlockSize = ::Opm::gpuistl::detail::getCudaRecomendedThreadBlockSize( cuSolveUpperLevelSetSplit, thrBlockSize); int nThreadBlocks = ::Opm::gpuistl::detail::getNumberOfBlocks(rowsInLevelSet, threadBlockSize); - cuSolveUpperLevelSetSplit<<>>( - reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); + cuSolveUpperLevelSetSplit + <<>>( + reorderedMat, rowIndices, colIndices, indexConversion, startIdx, rowsInLevelSet, dInv, v); } template @@ -484,28 +490,28 @@ LUFactorizationSplit(InputScalar* srcReorderedLowerMat, } #define INSTANTIATE_KERNEL_WRAPPERS(T, blocksize) \ - template void solveUpperLevelSet(T*, int*, int*, int*, int, int, T*, int); \ - template void solveLowerLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int); \ + template void solveUpperLevelSet(T*, int*, int*, int*, int, int, T*, int, cudaStream_t); \ + template void solveLowerLevelSet(T*, int*, int*, int*, int, int, const T*, T*, int, cudaStream_t); \ template void LUFactorization(T*, int*, int*, int*, int*, size_t, int, int); \ - template void LUFactorizationSplit( \ + template void LUFactorizationSplit( \ T*, int*, int*, T*, int*, int*, T*, float*, float*, float*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ + template void LUFactorizationSplit( \ T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ + template void LUFactorizationSplit( \ T*, int*, int*, T*, int*, int*, T*, float*, float*, float*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ - T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ + template void LUFactorizationSplit( \ + T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); \ + template void LUFactorizationSplit( \ T*, int*, int*, T*, int*, int*, T*, float*, float*, float*, int*, int*, const int, int, int); \ - template void LUFactorizationSplit( \ - T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); + template void LUFactorizationSplit( \ + T*, int*, int*, T*, int*, int*, T*, double*, double*, double*, int*, int*, const int, int, int); -#define INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(T) \ - INSTANTIATE_KERNEL_WRAPPERS(T, 1); \ - INSTANTIATE_KERNEL_WRAPPERS(T, 2); \ - INSTANTIATE_KERNEL_WRAPPERS(T, 3); \ - INSTANTIATE_KERNEL_WRAPPERS(T, 4); \ - INSTANTIATE_KERNEL_WRAPPERS(T, 5); \ +#define INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(T) \ + INSTANTIATE_KERNEL_WRAPPERS(T, 1); \ + INSTANTIATE_KERNEL_WRAPPERS(T, 2); \ + INSTANTIATE_KERNEL_WRAPPERS(T, 3); \ + INSTANTIATE_KERNEL_WRAPPERS(T, 4); \ + INSTANTIATE_KERNEL_WRAPPERS(T, 5); \ INSTANTIATE_KERNEL_WRAPPERS(T, 6); INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(float) @@ -514,32 +520,32 @@ INSTANTIATE_BLOCK_SIZED_KERNEL_WRAPPERS(double) #define INSTANTIATE_MIXED_PRECISION_KERNEL_WRAPPERS(blocksize) \ /* double preconditioner */ \ template void solveLowerLevelSetSplit( \ - double*, int*, int*, int*, int, int, const double*, double*, int); \ + double*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t); \ /* float matrix, double compute preconditioner */ \ template void solveLowerLevelSetSplit( \ - float*, int*, int*, int*, int, int, const double*, double*, int); \ + float*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t); \ /* float preconditioner */ \ template void solveLowerLevelSetSplit( \ - float*, int*, int*, int*, int, int, const float*, float*, int); \ + float*, int*, int*, int*, int, int, const float*, float*, int, cudaStream_t); \ \ /* double preconditioner */ \ - template void solveUpperLevelSetSplit( \ - double*, int*, int*, int*, int, int, const double*, double*, int); \ + template void solveUpperLevelSetSplit( \ + double*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t); \ /* float matrix, double compute preconditioner */ \ - template void solveUpperLevelSetSplit( \ - float*, int*, int*, int*, int, int, const double*, double*, int); \ + template void solveUpperLevelSetSplit( \ + float*, int*, int*, int*, int, int, const double*, double*, int, cudaStream_t); \ /* float preconditioner */ \ - template void solveUpperLevelSetSplit( \ - float*, int*, int*, int*, int, int, const double*, float*, int); \ + template void solveUpperLevelSetSplit( \ + float*, int*, int*, int*, int, int, const double*, float*, int, cudaStream_t); \ /* double preconditioner */ \ - template void solveUpperLevelSetSplit( \ - double*, int*, int*, int*, int, int, const float*, double*, int); \ + template void solveUpperLevelSetSplit( \ + double*, int*, int*, int*, int, int, const float*, double*, int, cudaStream_t); \ /* float matrix, double compute preconditioner */ \ - template void solveUpperLevelSetSplit( \ - float*, int*, int*, int*, int, int, const float*, double*, int); \ + template void solveUpperLevelSetSplit( \ + float*, int*, int*, int*, int, int, const float*, double*, int, cudaStream_t); \ /* float preconditioner */ \ - template void solveUpperLevelSetSplit( \ - float*, int*, int*, int*, int, int, const float*, float*, int); + template void solveUpperLevelSetSplit( \ + float*, int*, int*, int*, int, int, const float*, float*, int, cudaStream_t); INSTANTIATE_MIXED_PRECISION_KERNEL_WRAPPERS(1); diff --git a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp index cdd78ec91..04390ff77 100644 --- a/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp +++ b/opm/simulators/linalg/gpuistl/detail/preconditionerKernels/ILU0Kernels.hpp @@ -46,7 +46,8 @@ void solveUpperLevelSet(T* reorderedMat, int startIdx, int rowsInLevelSet, T* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform a lower solve on certain rows in a matrix that can safely be computed in parallel @@ -72,7 +73,8 @@ void solveLowerLevelSet(T* reorderedMat, int rowsInLevelSet, const T* d, T* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform an upper solve on certain rows in a matrix that can safely be computed in parallel @@ -99,7 +101,8 @@ void solveUpperLevelSetSplit(MatrixScalar* reorderedMat, int rowsInLevelSet, const DiagonalScalar* dInv, LinearSolverScalar* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Perform an lower solve on certain rows in a matrix that can safely be computed in parallel @@ -127,7 +130,8 @@ void solveLowerLevelSetSplit(MatrixScalar* reorderedLowerMat, int rowsInLevelSet, const LinearSolverScalar* d, LinearSolverScalar* v, - int threadBlockSize); + int threadBlockSize, + cudaStream_t stream); /** * @brief Computes the ILU Factorization of the input bcsr matrix, which is stored in a reordered way. The diagonal From ec1359c9159db84d1c8066d1b84e7f3a8595e50d Mon Sep 17 00:00:00 2001 From: Halvor M Nilsen Date: Tue, 4 Feb 2025 12:37:39 +0100 Subject: [PATCH 14/36] Change primary variables to be consistent when updated from wellstate --- opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp | 2 ++ opm/simulators/wells/StandardWellPrimaryVariables.cpp | 3 +++ 2 files changed, 5 insertions(+) diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index 1c2a99507..c1722cb40 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -152,6 +152,7 @@ update(const WellState& well_state, } } } + init(); } template @@ -208,6 +209,7 @@ updateNewton(const BVectorWell& dwells, if (stop_or_zero_rate_target) { value_[0][WQTotal] = 0.; } + init(); } template diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp index 17410e2f3..0c43bf9e9 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.cpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.cpp @@ -223,6 +223,7 @@ update(const WellState& well_state, // BHP value_[Bhp] = ws.bhp; + init(); } template @@ -301,6 +302,7 @@ updateNewton(const BVectorWell& dwells, // so that bhp constaint can be an active control when needed. constexpr Scalar bhp_lower_limit = 1. * unit::barsa - 1. * unit::Pascal; value_[Bhp] = std::max(value_[Bhp] - dx1_limited, bhp_lower_limit); + init(); } template @@ -320,6 +322,7 @@ updateNewtonPolyMW(const BVectorWell& dwells) value_[pskin_index] -= relaxation_factor * dx_pskin; } } + init(); } template From e7a2d954fe4bd0bbc66d31c38e586c6ebb10c6b0 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Thu, 13 Feb 2025 11:08:14 +0100 Subject: [PATCH 15/36] 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); From de1e2480f883640ec291e7b912c455afcb9376d7 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 7 Jan 2025 10:40:18 +0100 Subject: [PATCH 16/36] replace typedefs with type aliases --- opm/simulators/wells/BlackoilWellModel.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index a13874990..c5b665ff6 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -118,13 +118,13 @@ template class WellContributions; // TODO: where we should put these types, WellInterface or Well Model? // or there is some other strategy, like TypeTag - typedef Dune::FieldVector VectorBlockType; - typedef Dune::BlockVector BVector; + using VectorBlockType = Dune::FieldVector; + using BVector = Dune::BlockVector; - typedef BlackOilPolymerModule PolymerModule; - typedef BlackOilMICPModule MICPModule; + using PolymerModule = BlackOilPolymerModule; + using MICPModule = BlackOilMICPModule; - // For the conversion between the surface volume rate and resrevoir voidage rate + // For the conversion between the surface volume rate and reservoir voidage rate using RateConverterType = RateConverter:: SurfaceToReservoirVoidage >; From 4bfadcd80d75ccace784b96ed9525e2c80cdcd12 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 4 Feb 2025 16:31:52 +0100 Subject: [PATCH 17/36] add CompositionalContainer, a container for compositional data output start by moving data members into it --- CMakeLists_files.cmake | 2 + .../flow/CompositionalContainer.cpp | 29 +++++++++++ .../flow/CompositionalContainer.hpp | 52 +++++++++++++++++++ .../flow/GenericOutputBlackoilModule.cpp | 12 ++--- .../flow/GenericOutputBlackoilModule.hpp | 6 +-- .../flow/OutputCompositionalModule.hpp | 16 +++--- 6 files changed, 101 insertions(+), 16 deletions(-) create mode 100644 opm/simulators/flow/CompositionalContainer.cpp create mode 100644 opm/simulators/flow/CompositionalContainer.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 51372468d..f4aa0dec2 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -90,6 +90,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/flow/BlackoilModelParameters.cpp opm/simulators/flow/BlackoilModelConvergenceMonitor.cpp opm/simulators/flow/CollectDataOnIORank.cpp + opm/simulators/flow/CompositionalContainer.cpp opm/simulators/flow/ConvergenceOutputConfiguration.cpp opm/simulators/flow/EclGenericWriter.cpp opm/simulators/flow/ExtboContainer.cpp @@ -821,6 +822,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/BlackoilModelProperties.hpp opm/simulators/flow/CollectDataOnIORank.hpp opm/simulators/flow/CollectDataOnIORank_impl.hpp + opm/simulators/flow/CompositionalContainer.hpp opm/simulators/flow/ConvergenceOutputConfiguration.hpp opm/simulators/flow/countGlobalCells.hpp opm/simulators/flow/CpGridVanguard.hpp diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp new file mode 100644 index 000000000..ff82dd9de --- /dev/null +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -0,0 +1,29 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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. +*/ + +#include +#include + + +namespace Opm { + +} // namespace Opm diff --git a/opm/simulators/flow/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp new file mode 100644 index 000000000..a9c885107 --- /dev/null +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -0,0 +1,52 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 + * \copydoc Opm::OutputBlackOilModule + */ +#ifndef OPM_COMPOSITIONAL_CONTAINER_HPP +#define OPM_COMPOSITIONAL_CONTAINER_HPP + +#include +#include + +namespace Opm { + +template +class CompositionalContainer +{ + using Scalar = typename FluidSystem::Scalar; + using ScalarBuffer = std::vector; + + static constexpr int numComponents = FluidSystem::numComponents; + static constexpr int numPhases = FluidSystem::numPhases; + +public: + // total mole fractions for each component + std::array moleFractions_; + // mole fractions for each component in each phase + std::array, numPhases> phaseMoleFractions_; +}; + +} // namespace Opm + +#endif // OPM_COMPOSITIONAL_CONTAINER_HPP diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 46f556241..4608cb071 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -536,21 +536,21 @@ assignToSolution(data::Solution& sol) // ZMF for (int i = 0; i < numComponents; ++i) { const auto name = fmt::format("ZMF{}", i + 1); // Generate ZMF1, ZMF2, ... - compositionalEntries.emplace_back(name, UnitSystem::measure::identity, moleFractions_[i]); + compositionalEntries.emplace_back(name, UnitSystem::measure::identity, compC_.moleFractions_[i]); } // XMF for (int i = 0; i < numComponents; ++i) { const auto name = fmt::format("XMF{}", i + 1); // Generate XMF1, XMF2, ... compositionalEntries.emplace_back(name, UnitSystem::measure::identity, - phaseMoleFractions_[oilPhaseIdx][i]); + compC_.phaseMoleFractions_[oilPhaseIdx][i]); } // YMF for (int i = 0; i < numComponents; ++i) { const auto name = fmt::format("YMF{}", i + 1); // Generate YMF1, YMF2, ... compositionalEntries.emplace_back(name, UnitSystem::measure::identity, - phaseMoleFractions_[gasPhaseIdx][i]); + compC_.phaseMoleFractions_[gasPhaseIdx][i]); } } @@ -1297,21 +1297,21 @@ doAllocBuffers(const unsigned bufferSize, if (rstKeywords["ZMF"] > 0) { rstKeywords["ZMF"] = 0; for (int i = 0; i < numComponents; ++i) { - moleFractions_[i].resize(bufferSize, 0.0); + compC_.moleFractions_[i].resize(bufferSize, 0.0); } } if (rstKeywords["XMF"] > 0 && FluidSystem::phaseIsActive(oilPhaseIdx)) { rstKeywords["XMF"] = 0; for (int i = 0; i < numComponents; ++i) { - phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0); + compC_.phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0); } } if (rstKeywords["YMF"] > 0 && FluidSystem::phaseIsActive(gasPhaseIdx)) { rstKeywords["YMF"] = 0; for (int i = 0; i < numComponents; ++i) { - phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0); + compC_.phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0); } } } diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index c2d45aafa..0a7b0a11b 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -32,6 +32,7 @@ #include #include +#include #include #include #include @@ -468,10 +469,7 @@ protected: std::array viscosity_; std::array relativePermeability_; - // total mole fractions for each component - std::array moleFractions_; - // mole fractions for each component in each phase - std::array, numPhases> phaseMoleFractions_; + CompositionalContainer compC_; std::vector freeTracerConcentrations_; std::vector solTracerConcentrations_; diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index d6102528d..5f369f47d 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -32,12 +32,16 @@ #include #include +#include #include #include #include #include + +#include +#include #include #include @@ -188,19 +192,19 @@ public: } for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - if (this->moleFractions_[compIdx].empty()) continue; + if (this->compC_.moleFractions_[compIdx].empty()) continue; - this->moleFractions_[compIdx][globalDofIdx] = getValue(fs.moleFraction(compIdx)); + this->compC_.moleFractions_[compIdx][globalDofIdx] = getValue(fs.moleFraction(compIdx)); } // XMF and YMF for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { if (FluidSystem::phaseIsActive(oilPhaseIdx)) { - if (this->phaseMoleFractions_[oilPhaseIdx][compIdx].empty()) continue; - this->phaseMoleFractions_[oilPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(oilPhaseIdx, compIdx)); + if (this->compC_.phaseMoleFractions_[oilPhaseIdx][compIdx].empty()) continue; + this->compC_.phaseMoleFractions_[oilPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(oilPhaseIdx, compIdx)); } if (FluidSystem::phaseIsActive(gasPhaseIdx)) { - if (this->phaseMoleFractions_[gasPhaseIdx][compIdx].empty()) continue; - this->phaseMoleFractions_[gasPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(gasPhaseIdx, compIdx)); + if (this->compC_.phaseMoleFractions_[gasPhaseIdx][compIdx].empty()) continue; + this->compC_.phaseMoleFractions_[gasPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(gasPhaseIdx, compIdx)); } } From 686f110a3489f7837873cc377e3134c844ad4259 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 11 Feb 2025 14:58:09 +0100 Subject: [PATCH 18/36] move allocation of compositional data into CompositionalContainer --- .../flow/CompositionalContainer.cpp | 56 +++++++++++++++++++ .../flow/CompositionalContainer.hpp | 10 ++++ .../flow/GenericOutputBlackoilModule.cpp | 21 +------ 3 files changed, 67 insertions(+), 20 deletions(-) diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp index ff82dd9de..e1f196e88 100644 --- a/opm/simulators/flow/CompositionalContainer.cpp +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -23,7 +23,63 @@ #include #include +#include +#include +#include namespace Opm { +template +void CompositionalContainer:: +allocate(const unsigned bufferSize, + std::map& rstKeywords) +{ + if (auto& zmf = rstKeywords["ZMF"]; zmf > 0) { + zmf = 0; + this->allocated_ = true; + for (int i = 0; i < numComponents; ++i) { + moleFractions_[i].resize(bufferSize, 0.0); + } + } + + if (auto& xmf = rstKeywords["XMF"]; xmf > 0 && FluidSystem::phaseIsActive(oilPhaseIdx)) { + this->allocated_ = true; + xmf = 0; + for (int i = 0; i < numComponents; ++i) { + phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0); + } + } + + if (auto& ymf = rstKeywords["YMF"]; ymf > 0 && FluidSystem::phaseIsActive(gasPhaseIdx)) { + this->allocated_ = true; + ymf = 0; + for (int i = 0; i < numComponents; ++i) { + phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0); + } + } +} + +template using FS = BlackOilFluidSystem; + +#define INSTANTIATE_TYPE(T) \ + template class CompositionalContainer>; + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +#define INSTANTIATE_COMP(NUM) \ + template using FS##NUM = GenericOilGasFluidSystem; \ + template class CompositionalContainer>; + +INSTANTIATE_COMP(0) +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/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp index a9c885107..07325bce9 100644 --- a/opm/simulators/flow/CompositionalContainer.hpp +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -27,6 +27,8 @@ #define OPM_COMPOSITIONAL_CONTAINER_HPP #include +#include +#include #include namespace Opm { @@ -38,9 +40,17 @@ class CompositionalContainer using ScalarBuffer = std::vector; static constexpr int numComponents = FluidSystem::numComponents; + static constexpr int numPhases = FluidSystem::numPhases; + static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx; + static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx; + static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx; public: + void allocate(const unsigned bufferSize, + std::map& rstKeywords); + + bool allocated_ = false; // total mole fractions for each component std::array moleFractions_; // mole fractions for each component in each phase diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 4608cb071..bd63f1c39 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -1294,26 +1294,7 @@ doAllocBuffers(const unsigned bufferSize, } if (this->isCompositional_) { - if (rstKeywords["ZMF"] > 0) { - rstKeywords["ZMF"] = 0; - for (int i = 0; i < numComponents; ++i) { - compC_.moleFractions_[i].resize(bufferSize, 0.0); - } - } - - if (rstKeywords["XMF"] > 0 && FluidSystem::phaseIsActive(oilPhaseIdx)) { - rstKeywords["XMF"] = 0; - for (int i = 0; i < numComponents; ++i) { - compC_.phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0); - } - } - - if (rstKeywords["YMF"] > 0 && FluidSystem::phaseIsActive(gasPhaseIdx)) { - rstKeywords["YMF"] = 0; - for (int i = 0; i < numComponents; ++i) { - compC_.phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0); - } - } + this->compC_.allocate(bufferSize, rstKeywords); } From f4b1c093ce43f6b8af030067ca2973b5ba7f79e5 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 11 Feb 2025 15:08:58 +0100 Subject: [PATCH 19/36] move output of restart data to CompositionalContainer --- .../flow/CompositionalContainer.cpp | 62 +++++++++++++++++++ .../flow/CompositionalContainer.hpp | 4 ++ .../flow/GenericOutputBlackoilModule.cpp | 30 +-------- 3 files changed, 67 insertions(+), 29 deletions(-) diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp index e1f196e88..a7343c921 100644 --- a/opm/simulators/flow/CompositionalContainer.cpp +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -27,6 +27,13 @@ #include #include +#include + +#include +#include + +#include + namespace Opm { template @@ -59,6 +66,61 @@ allocate(const unsigned bufferSize, } } +template +void CompositionalContainer:: +outputRestart(data::Solution& sol) +{ + using DataEntry = + std::tuple&>; + + auto doInsert = [&sol](DataEntry& entry, + const data::TargetType target) + { + if (std::get<2>(entry).empty()) { + return; + } + + sol.insert(std::get(entry), + std::get(entry), + std::move(std::get<2>(entry)), + target); + }; + + auto entries = std::vector{}; + + // ZMF + if (!moleFractions_[0].empty()) { + for (int i = 0; i < numComponents; ++i) { + const auto name = fmt::format("ZMF{}", i + 1); // Generate ZMF1, ZMF2, ... + entries.emplace_back(name, UnitSystem::measure::identity, moleFractions_[i]); + } + } + + // XMF + if (!phaseMoleFractions_[oilPhaseIdx][0].empty()) { + for (int i = 0; i < numComponents; ++i) { + const auto name = fmt::format("XMF{}", i + 1); // Generate XMF1, XMF2, ... + entries.emplace_back(name, UnitSystem::measure::identity, + phaseMoleFractions_[oilPhaseIdx][i]); + } + } + + // YMF + if (!phaseMoleFractions_[gasPhaseIdx][0].empty()) { + for (int i = 0; i < numComponents; ++i) { + const auto name = fmt::format("YMF{}", i + 1); // Generate YMF1, YMF2, ... + entries.emplace_back(name, UnitSystem::measure::identity, + phaseMoleFractions_[gasPhaseIdx][i]); + } + } + + std::for_each(entries.begin(), entries.end(), + [&doInsert](auto& array) + { doInsert(array, data::TargetType::RESTART_SOLUTION); }); + + this->allocated_ = false; +} + template using FS = BlackOilFluidSystem; #define INSTANTIATE_TYPE(T) \ diff --git a/opm/simulators/flow/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp index 07325bce9..feca23ea6 100644 --- a/opm/simulators/flow/CompositionalContainer.hpp +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -33,6 +33,8 @@ namespace Opm { +namespace data { class Solution; } + template class CompositionalContainer { @@ -50,6 +52,8 @@ public: void allocate(const unsigned bufferSize, std::map& rstKeywords); + void outputRestart(data::Solution& sol); + bool allocated_ = false; // total mole fractions for each component std::array moleFractions_; diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index bd63f1c39..405d70887 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -527,36 +527,8 @@ assignToSolution(data::Solution& sol) DataEntry{"TMULT_RC", UnitSystem::measure::identity, rockCompTransMultiplier_}, }; - // basically, for compositional, we can not use std::array for this. We need to generate the ZMF1, ZMF2, and so on - // and also, we need to map these values. - // TODO: the following should go to a function if (this->isCompositional_) { - auto compositionalEntries = std::vector{}; - { - // ZMF - for (int i = 0; i < numComponents; ++i) { - const auto name = fmt::format("ZMF{}", i + 1); // Generate ZMF1, ZMF2, ... - compositionalEntries.emplace_back(name, UnitSystem::measure::identity, compC_.moleFractions_[i]); - } - - // XMF - for (int i = 0; i < numComponents; ++i) { - const auto name = fmt::format("XMF{}", i + 1); // Generate XMF1, XMF2, ... - compositionalEntries.emplace_back(name, UnitSystem::measure::identity, - compC_.phaseMoleFractions_[oilPhaseIdx][i]); - } - - // YMF - for (int i = 0; i < numComponents; ++i) { - const auto name = fmt::format("YMF{}", i + 1); // Generate YMF1, YMF2, ... - compositionalEntries.emplace_back(name, UnitSystem::measure::identity, - compC_.phaseMoleFractions_[gasPhaseIdx][i]); - } - } - - for (auto& array: compositionalEntries) { - doInsert(array, data::TargetType::RESTART_SOLUTION); - } + this->compC_.outputRestart(sol); } for (auto& array : baseSolutionVector) { From 0496f4bb8ac29c24ac1bea8eac008c007d507578 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 12 Feb 2025 07:51:22 +0100 Subject: [PATCH 20/36] move assignment of mole fractions into CompositionalContainer --- opm/simulators/flow/CompositionalContainer.cpp | 14 ++++++++++++++ opm/simulators/flow/CompositionalContainer.hpp | 9 +++++++++ opm/simulators/flow/OutputCompositionalModule.hpp | 8 ++++---- 3 files changed, 27 insertions(+), 4 deletions(-) diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp index a7343c921..b6387e91d 100644 --- a/opm/simulators/flow/CompositionalContainer.cpp +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -66,6 +66,20 @@ allocate(const unsigned bufferSize, } } +template +void CompositionalContainer:: +assignMoleFractions(const unsigned globalDofIdx, + const AssignFunction& fractions) +{ + if (moleFractions_.empty()) { + return; + } + + std::for_each(moleFractions_.begin(), moleFractions_.end(), + [&fractions, globalDofIdx, c = 0](auto& comp) mutable + { comp[globalDofIdx] = fractions(c++); }); +} + template void CompositionalContainer:: outputRestart(data::Solution& sol) diff --git a/opm/simulators/flow/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp index feca23ea6..f3aef1e42 100644 --- a/opm/simulators/flow/CompositionalContainer.hpp +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -27,6 +27,7 @@ #define OPM_COMPOSITIONAL_CONTAINER_HPP #include +#include #include #include #include @@ -52,8 +53,16 @@ public: void allocate(const unsigned bufferSize, std::map& rstKeywords); + using AssignFunction = std::function; + + void assignMoleFractions(const unsigned globalDofIdx, + const AssignFunction& fractions); + void outputRestart(data::Solution& sol); + bool allocated() const + { return allocated_; } + bool allocated_ = false; // total mole fractions for each component std::array moleFractions_; diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index 5f369f47d..45729f4f2 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -191,10 +191,10 @@ public: Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]); } - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - if (this->compC_.moleFractions_[compIdx].empty()) continue; - - this->compC_.moleFractions_[compIdx][globalDofIdx] = getValue(fs.moleFraction(compIdx)); + if (this->compC_.allocated()) { + this->compC_.assignMoleFractions(globalDofIdx, + [&fs](const unsigned compIdx) + { return getValue(fs.moleFraction(compIdx)); }); } // XMF and YMF for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { From ae23b39d40cf0be7b285174d6078df55fe986d8c Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 12 Feb 2025 07:51:22 +0100 Subject: [PATCH 21/36] move assignment of oil fractions into CompositionalContainer --- opm/simulators/flow/CompositionalContainer.cpp | 14 ++++++++++++++ opm/simulators/flow/CompositionalContainer.hpp | 3 +++ opm/simulators/flow/OutputCompositionalModule.hpp | 10 ++++++---- 3 files changed, 23 insertions(+), 4 deletions(-) diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp index b6387e91d..f9bd6fc40 100644 --- a/opm/simulators/flow/CompositionalContainer.cpp +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -80,6 +80,20 @@ assignMoleFractions(const unsigned globalDofIdx, { comp[globalDofIdx] = fractions(c++); }); } +template +void CompositionalContainer:: +assignOilFractions(const unsigned globalDofIdx, + const AssignFunction& fractions) +{ + if (phaseMoleFractions_[oilPhaseIdx][0].empty()) { + return; + } + std::for_each(phaseMoleFractions_[oilPhaseIdx].begin(), + phaseMoleFractions_[oilPhaseIdx].end(), + [globalDofIdx, &fractions, c = 0](auto& comp) mutable + { comp[globalDofIdx] = fractions(c++); }); +} + template void CompositionalContainer:: outputRestart(data::Solution& sol) diff --git a/opm/simulators/flow/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp index f3aef1e42..f7e4aca84 100644 --- a/opm/simulators/flow/CompositionalContainer.hpp +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -58,6 +58,9 @@ public: void assignMoleFractions(const unsigned globalDofIdx, const AssignFunction& fractions); + void assignOilFractions(const unsigned globalDofIdx, + const AssignFunction& fractions); + void outputRestart(data::Solution& sol); bool allocated() const diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index 45729f4f2..db8dcec9b 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -195,13 +195,15 @@ public: this->compC_.assignMoleFractions(globalDofIdx, [&fs](const unsigned compIdx) { return getValue(fs.moleFraction(compIdx)); }); + + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { + this->compC_.assignOilFractions(globalDofIdx, + [&fs](const unsigned compIdx) + { return getValue(fs.moleFraction(oilPhaseIdx, compIdx)); }); + } } // XMF and YMF for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - if (FluidSystem::phaseIsActive(oilPhaseIdx)) { - if (this->compC_.phaseMoleFractions_[oilPhaseIdx][compIdx].empty()) continue; - this->compC_.phaseMoleFractions_[oilPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(oilPhaseIdx, compIdx)); - } if (FluidSystem::phaseIsActive(gasPhaseIdx)) { if (this->compC_.phaseMoleFractions_[gasPhaseIdx][compIdx].empty()) continue; this->compC_.phaseMoleFractions_[gasPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(gasPhaseIdx, compIdx)); From 95e37356780e888e6ea347204b7b4aaaa9ea2bfe Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 12 Feb 2025 07:51:22 +0100 Subject: [PATCH 22/36] move assignment of gas fractions into CompositionalContainer --- opm/simulators/flow/CompositionalContainer.cpp | 16 ++++++++++++++++ opm/simulators/flow/CompositionalContainer.hpp | 3 +++ .../flow/OutputCompositionalModule.hpp | 16 +++++++--------- 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp index f9bd6fc40..8f5034b0c 100644 --- a/opm/simulators/flow/CompositionalContainer.cpp +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -66,6 +66,21 @@ allocate(const unsigned bufferSize, } } +template +void CompositionalContainer:: +assignGasFractions(const unsigned globalDofIdx, + const AssignFunction& fractions) +{ + if (phaseMoleFractions_[gasPhaseIdx][0].empty()) { + return; + } + + std::for_each(phaseMoleFractions_[gasPhaseIdx].begin(), + phaseMoleFractions_[gasPhaseIdx].end(), + [globalDofIdx, &fractions, c = 0](auto& comp) mutable + { comp[globalDofIdx] = fractions(c++); }); +} + template void CompositionalContainer:: assignMoleFractions(const unsigned globalDofIdx, @@ -88,6 +103,7 @@ assignOilFractions(const unsigned globalDofIdx, if (phaseMoleFractions_[oilPhaseIdx][0].empty()) { return; } + std::for_each(phaseMoleFractions_[oilPhaseIdx].begin(), phaseMoleFractions_[oilPhaseIdx].end(), [globalDofIdx, &fractions, c = 0](auto& comp) mutable diff --git a/opm/simulators/flow/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp index f7e4aca84..9afb7eee8 100644 --- a/opm/simulators/flow/CompositionalContainer.hpp +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -55,6 +55,9 @@ public: using AssignFunction = std::function; + void assignGasFractions(const unsigned globalDofIdx, + const AssignFunction& fractions); + void assignMoleFractions(const unsigned globalDofIdx, const AssignFunction& fractions); diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index db8dcec9b..a4fdd6ccb 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -57,8 +57,7 @@ #include -namespace Opm -{ +namespace Opm { // forward declaration template @@ -196,19 +195,18 @@ public: [&fs](const unsigned compIdx) { return getValue(fs.moleFraction(compIdx)); }); + if (FluidSystem::phaseIsActive(gasPhaseIdx)) { + this->compC_.assignGasFractions(globalDofIdx, + [&fs](const unsigned compIdx) + { return getValue(fs.moleFraction(gasPhaseIdx, compIdx)); }); + } + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { this->compC_.assignOilFractions(globalDofIdx, [&fs](const unsigned compIdx) { return getValue(fs.moleFraction(oilPhaseIdx, compIdx)); }); } } - // XMF and YMF - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - if (FluidSystem::phaseIsActive(gasPhaseIdx)) { - if (this->compC_.phaseMoleFractions_[gasPhaseIdx][compIdx].empty()) continue; - this->compC_.phaseMoleFractions_[gasPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(gasPhaseIdx, compIdx)); - } - } if (!this->fluidPressure_.empty()) { if (FluidSystem::phaseIsActive(oilPhaseIdx)) { From 3b819349f8d0724d40c5b9e6692359a67194bb28 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 12 Feb 2025 13:53:09 +0100 Subject: [PATCH 23/36] CompositionalContainer: make data private --- opm/simulators/flow/CompositionalContainer.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/simulators/flow/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp index 9afb7eee8..f2b2a2a2f 100644 --- a/opm/simulators/flow/CompositionalContainer.hpp +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -69,6 +69,7 @@ public: bool allocated() const { return allocated_; } +private: bool allocated_ = false; // total mole fractions for each component std::array moleFractions_; From c8209c39a79e87ac6d15c7551f8a862f37d1bf58 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 12 Feb 2025 15:07:36 +0100 Subject: [PATCH 24/36] move CompositionalContainer member into OutputCompositionalModule --- .../flow/CompositionalContainer.cpp | 22 +++++--------- .../flow/CompositionalContainer.hpp | 3 +- .../flow/GenericOutputBlackoilModule.cpp | 30 +++++-------------- .../flow/GenericOutputBlackoilModule.hpp | 11 +++---- .../flow/OutputCompositionalModule.hpp | 28 ++++++++++++----- 5 files changed, 41 insertions(+), 53 deletions(-) diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp index 8f5034b0c..28135641d 100644 --- a/opm/simulators/flow/CompositionalContainer.cpp +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -23,8 +23,6 @@ #include #include -#include -#include #include #include @@ -42,8 +40,8 @@ allocate(const unsigned bufferSize, std::map& rstKeywords) { if (auto& zmf = rstKeywords["ZMF"]; zmf > 0) { - zmf = 0; this->allocated_ = true; + zmf = 0; for (int i = 0; i < numComponents; ++i) { moleFractions_[i].resize(bufferSize, 0.0); } @@ -112,7 +110,8 @@ assignOilFractions(const unsigned globalDofIdx, template void CompositionalContainer:: -outputRestart(data::Solution& sol) +outputRestart(data::Solution& sol, + ScalarBuffer& oil_saturation) { using DataEntry = std::tuple&>; @@ -158,6 +157,10 @@ outputRestart(data::Solution& sol) } } + if (!oil_saturation.empty()) { + entries.emplace_back("SOIL", UnitSystem::measure::identity, oil_saturation); + } + std::for_each(entries.begin(), entries.end(), [&doInsert](auto& array) { doInsert(array, data::TargetType::RESTART_SOLUTION); }); @@ -165,17 +168,6 @@ outputRestart(data::Solution& sol) this->allocated_ = false; } -template using FS = BlackOilFluidSystem; - -#define INSTANTIATE_TYPE(T) \ - template class CompositionalContainer>; - -INSTANTIATE_TYPE(double) - -#if FLOW_INSTANTIATE_FLOAT -INSTANTIATE_TYPE(float) -#endif - #define INSTANTIATE_COMP(NUM) \ template using FS##NUM = GenericOilGasFluidSystem; \ template class CompositionalContainer>; diff --git a/opm/simulators/flow/CompositionalContainer.hpp b/opm/simulators/flow/CompositionalContainer.hpp index f2b2a2a2f..67a9ca4f9 100644 --- a/opm/simulators/flow/CompositionalContainer.hpp +++ b/opm/simulators/flow/CompositionalContainer.hpp @@ -64,7 +64,8 @@ public: void assignOilFractions(const unsigned globalDofIdx, const AssignFunction& fractions); - void outputRestart(data::Solution& sol); + void outputRestart(data::Solution& sol, + ScalarBuffer& oil_saturation); bool allocated() const { return allocated_; } diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 405d70887..1d4a19fb3 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -104,8 +104,7 @@ GenericOutputBlackoilModule(const EclipseState& eclState, bool enableBrine, bool enableSaltPrecipitation, bool enableExtbo, - bool enableMICP, - bool isCompositional) + bool enableMICP) : eclState_(eclState) , schedule_(schedule) , summaryState_(summaryState) @@ -124,7 +123,6 @@ GenericOutputBlackoilModule(const EclipseState& eclState, , enableSaltPrecipitation_(enableSaltPrecipitation) , enableExtbo_(enableExtbo) , enableMICP_(enableMICP) - , isCompositional_(isCompositional) , local_data_valid_(false) { const auto& fp = eclState_.fieldProps(); @@ -527,10 +525,6 @@ assignToSolution(data::Solution& sol) DataEntry{"TMULT_RC", UnitSystem::measure::identity, rockCompTransMultiplier_}, }; - if (this->isCompositional_) { - this->compC_.outputRestart(sol); - } - for (auto& array : baseSolutionVector) { doInsert(array, data::TargetType::RESTART_SOLUTION); } @@ -574,15 +568,6 @@ assignToSolution(data::Solution& sol) data::TargetType::RESTART_SOLUTION); } - if (this->isCompositional_ && FluidSystem::phaseIsActive(oilPhaseIdx) && - ! this->saturation_[oilPhaseIdx].empty()) - { - sol.insert("SOIL", UnitSystem::measure::identity, - std::move(this->saturation_[oilPhaseIdx]), - data::TargetType::RESTART_SOLUTION); - } - - if ((eclState_.runspec().co2Storage() || eclState_.runspec().h2Storage()) && !rsw_.empty()) { auto mfrac = std::vector(this->rsw_.size(), 0.0); @@ -806,10 +791,14 @@ doAllocBuffers(const unsigned bufferSize, const bool enableWettingHysteresis, const unsigned numTracers, const std::vector& enableSolTracers, - const unsigned numOutputNnc) + const unsigned numOutputNnc, + std::map rstKeywords) { + if (rstKeywords.empty()) { + rstKeywords = schedule_.rst_keywords(reportStepNum); + } + // Output RESTART_OPM_EXTENDED only when explicitly requested by user. - std::map rstKeywords = schedule_.rst_keywords(reportStepNum); for (auto& [keyword, should_write] : rstKeywords) { if (this->isOutputCreationDirective_(keyword)) { // 'BASIC', 'FREQ' and similar. Don't attempt to create @@ -1265,11 +1254,6 @@ doAllocBuffers(const unsigned bufferSize, overburdenPressure_.resize(bufferSize, 0.0); } - if (this->isCompositional_) { - this->compC_.allocate(bufferSize, rstKeywords); - } - - //Warn for any unhandled keyword if (log) { for (auto& keyValue: rstKeywords) { diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index 0a7b0a11b..5c570803e 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -32,7 +32,6 @@ #include #include -#include #include #include #include @@ -320,8 +319,7 @@ protected: bool enableBrine, bool enableSaltPrecipitation, bool enableExtbo, - bool enableMICP, - bool isCompositional = false); + bool enableMICP); void doAllocBuffers(unsigned bufferSize, unsigned reportStepNum, @@ -330,11 +328,12 @@ protected: const bool isRestart, const bool vapparsActive = false, const bool enablePCHysteresis = false, - const bool enableNonWettingHysteresis =false, + const bool enableNonWettingHysteresis = false, const bool enableWettingHysteresis = false, unsigned numTracers = 0, const std::vector& enableSolTracers = {}, - unsigned numOutputNnc = 0); + unsigned numOutputNnc = 0, + std::map rstKeywords = {}); void makeRegionSum(Inplace& inplace, const std::string& region_name, @@ -391,7 +390,6 @@ protected: bool enableSaltPrecipitation_{false}; bool enableExtbo_{false}; bool enableMICP_{false}; - bool isCompositional_{false}; bool forceDisableFipOutput_{false}; bool forceDisableFipresvOutput_{false}; @@ -469,7 +467,6 @@ protected: std::array viscosity_; std::array relativePermeability_; - CompositionalContainer compC_; std::vector freeTracerConcentrations_; std::vector solTracerConcentrations_; diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index a4fdd6ccb..34b161164 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -45,6 +45,7 @@ #include #include +#include #include #include @@ -105,8 +106,7 @@ public: getPropValue(), getPropValue(), getPropValue(), - getPropValue(), - true) + getPropValue()) , simulator_(simulator) { for (auto& region_pair : this->regions_) { @@ -158,11 +158,24 @@ public: return; } - this->doAllocBuffers(bufferSize, - reportStepNum, - substep, - log, - isRestart); + auto rstKeywords = this->schedule_.rst_keywords(reportStepNum); + this->compC_.allocate(bufferSize, rstKeywords); + + this->doAllocBuffers(bufferSize, reportStepNum, substep, log, isRestart, + /* vapparsActive =*/ false, + /* enablePCHysteresis = */ false, + /* enableNonWettingHysteresis =*/ false, + /* enableWettingHysteresis =*/ false, + /* numTracers = */ 0, + /* enableSoltracers =*/ {}, + /* numOutputNnc =*/ 0, + std::move(rstKeywords)); + } + + void assignToSolution(data::Solution& sol) + { + this->compC_.outputRestart(sol, this->saturation_[oilPhaseIdx]); + BaseType::assignToSolution(sol); } /*! @@ -336,6 +349,7 @@ private: } const Simulator& simulator_; + CompositionalContainer compC_; }; } // namespace Opm From 41fbfb3d22ed5ef506d41a99a979d73ec6519a72 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 13 Feb 2025 13:13:59 +0100 Subject: [PATCH 25/36] fix: instantiation of CompositionalContainer --- opm/simulators/flow/CompositionalContainer.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/opm/simulators/flow/CompositionalContainer.cpp b/opm/simulators/flow/CompositionalContainer.cpp index 28135641d..0a51ac903 100644 --- a/opm/simulators/flow/CompositionalContainer.cpp +++ b/opm/simulators/flow/CompositionalContainer.cpp @@ -23,7 +23,7 @@ #include #include -#include +#include #include @@ -168,11 +168,19 @@ outputRestart(data::Solution& sol, this->allocated_ = false; } -#define INSTANTIATE_COMP(NUM) \ - template using FS##NUM = GenericOilGasFluidSystem; \ +#define INSTANTIATE_COMP_THREEPHASE(NUM) \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ template class CompositionalContainer>; -INSTANTIATE_COMP(0) +#define INSTANTIATE_COMP_TWOPHASE(NUM) \ + template using GFS##NUM = GenericOilGasWaterFluidSystem; \ + template class CompositionalContainer>; + +#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) From 13fb52bf9e40175d6e8e3de8c8ab070624a609e8 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 4 Feb 2025 16:31:52 +0100 Subject: [PATCH 26/36] add TracerContainer, a container for tracer data output start by moving data members into it --- CMakeLists_files.cmake | 2 + .../flow/GenericOutputBlackoilModule.cpp | 22 ++++----- .../flow/GenericOutputBlackoilModule.hpp | 6 +-- opm/simulators/flow/OutputBlackoilModule.hpp | 12 ++--- opm/simulators/flow/TracerContainer.cpp | 28 +++++++++++ opm/simulators/flow/TracerContainer.hpp | 46 +++++++++++++++++++ 6 files changed, 96 insertions(+), 20 deletions(-) create mode 100644 opm/simulators/flow/TracerContainer.cpp create mode 100644 opm/simulators/flow/TracerContainer.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f4aa0dec2..28d3fefa9 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -120,6 +120,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/flow/SimulatorReportBanners.cpp opm/simulators/flow/SimulatorSerializer.cpp opm/simulators/flow/SolutionContainers.cpp + opm/simulators/flow/TracerContainer.cpp opm/simulators/flow/Transmissibility.cpp opm/simulators/flow/ValidationFunctions.cpp opm/simulators/flow/equil/EquilibrationHelpers.cpp @@ -880,6 +881,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/SolutionContainers.hpp opm/simulators/flow/SubDomain.hpp opm/simulators/flow/TTagFlowProblemTPFA.hpp + opm/simulators/flow/TracerContainer.hpp opm/simulators/flow/TracerModel.hpp opm/simulators/flow/Transmissibility.hpp opm/simulators/flow/Transmissibility_impl.hpp diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 9c41bb916..8e9f8590a 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -631,40 +631,40 @@ assignToSolution(data::Solution& sol) this->fipC_.outputRestart(sol); // Tracers - if (! this->freeTracerConcentrations_.empty()) { + if (! this->tracerC_.freeConcentrations_.empty()) { const auto& tracers = this->eclState_.tracer(); for (auto tracerIdx = 0*tracers.size(); tracerIdx < tracers.size(); ++tracerIdx) { sol.insert(tracers[tracerIdx].fname(), UnitSystem::measure::identity, - std::move(freeTracerConcentrations_[tracerIdx]), + std::move(this->tracerC_.freeConcentrations_[tracerIdx]), data::TargetType::RESTART_TRACER_SOLUTION); } // Put freeTracerConcentrations container into a valid state. Otherwise // we'll move from vectors that have already been moved from if we // get here and it's not a restart step. - this->freeTracerConcentrations_.clear(); + this->tracerC_.freeConcentrations_.clear(); } - if (! this->solTracerConcentrations_.empty()) { + if (! this->tracerC_.solConcentrations_.empty()) { const auto& tracers = this->eclState_.tracer(); for (auto tracerIdx = 0*tracers.size(); tracerIdx < tracers.size(); ++tracerIdx) { - if (solTracerConcentrations_[tracerIdx].empty()) + if (this->tracerC_.solConcentrations_[tracerIdx].empty()) continue; sol.insert(tracers[tracerIdx].sname(), UnitSystem::measure::identity, - std::move(solTracerConcentrations_[tracerIdx]), + std::move(this->tracerC_.solConcentrations_[tracerIdx]), data::TargetType::RESTART_TRACER_SOLUTION); } // Put solTracerConcentrations container into a valid state. Otherwise // we'll move from vectors that have already been moved from if we // get here and it's not a restart step. - this->solTracerConcentrations_.clear(); + this->tracerC_.solConcentrations_.clear(); } } @@ -1221,16 +1221,16 @@ doAllocBuffers(const unsigned bufferSize, // tracers if (numTracers > 0) { - freeTracerConcentrations_.resize(numTracers); + tracerC_.freeConcentrations_.resize(numTracers); for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) { - freeTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0); + tracerC_.freeConcentrations_[tracerIdx].resize(bufferSize, 0.0); } - solTracerConcentrations_.resize(numTracers); + tracerC_.solConcentrations_.resize(numTracers); for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) { if (enableSolTracers[tracerIdx]) - solTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0); + tracerC_.solConcentrations_[tracerIdx].resize(bufferSize, 0.0); } } diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index 5c570803e..059f5a093 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -40,6 +40,7 @@ #include #include #include +#include #include @@ -330,7 +331,7 @@ protected: const bool enablePCHysteresis = false, const bool enableNonWettingHysteresis = false, const bool enableWettingHysteresis = false, - unsigned numTracers = 0, + const unsigned numTracers = 0, const std::vector& enableSolTracers = {}, unsigned numOutputNnc = 0, std::map rstKeywords = {}); @@ -467,8 +468,7 @@ protected: std::array viscosity_; std::array relativePermeability_; - std::vector freeTracerConcentrations_; - std::vector solTracerConcentrations_; + TracerContainer tracerC_; std::array residual_; diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 650cb8f3e..71d44ea29 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -647,21 +647,21 @@ public: // tracers const auto& tracerModel = simulator_.problem().tracerModel(); - if (! this->freeTracerConcentrations_.empty()) { + if (! this->tracerC_.freeConcentrations_.empty()) { for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { - if (this->freeTracerConcentrations_[tracerIdx].empty()) { + if (this->tracerC_.freeConcentrations_[tracerIdx].empty()) { continue; } - this->freeTracerConcentrations_[tracerIdx][globalDofIdx] = + this->tracerC_.freeConcentrations_[tracerIdx][globalDofIdx] = tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx); } } - if (! this->solTracerConcentrations_.empty()) { + if (! this->tracerC_.solConcentrations_.empty()) { for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { - if (this->solTracerConcentrations_[tracerIdx].empty()) { + if (this->tracerC_.solConcentrations_[tracerIdx].empty()) { continue; } - this->solTracerConcentrations_[tracerIdx][globalDofIdx] = + this->tracerC_.solConcentrations_[tracerIdx][globalDofIdx] = tracerModel.solTracerConcentration(tracerIdx, globalDofIdx); } diff --git a/opm/simulators/flow/TracerContainer.cpp b/opm/simulators/flow/TracerContainer.cpp new file mode 100644 index 000000000..6ceacbee1 --- /dev/null +++ b/opm/simulators/flow/TracerContainer.cpp @@ -0,0 +1,28 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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. +*/ + +#include +#include + +namespace Opm { + +} // namespace Opm diff --git a/opm/simulators/flow/TracerContainer.hpp b/opm/simulators/flow/TracerContainer.hpp new file mode 100644 index 000000000..17f40f501 --- /dev/null +++ b/opm/simulators/flow/TracerContainer.hpp @@ -0,0 +1,46 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 + * \copydoc Opm::OutputBlackOilModule + */ +#ifndef OPM_TRACER_CONTAINER_HPP +#define OPM_TRACER_CONTAINER_HPP + +#include + +namespace Opm { + +template +class TracerContainer +{ + using Scalar = typename FluidSystem::Scalar; + using ScalarBuffer = std::vector; + +public: + std::vector freeConcentrations_; + std::vector solConcentrations_; +}; + +} // namespace Opm + +#endif // OPM_TRACER_CONTAINER_HPP From 7f4e3624f7ad9fa3bdf912dbd0ea06bff3d6d538 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 14 Feb 2025 08:58:44 +0100 Subject: [PATCH 27/36] move allocation of tracer data into TracerContainer --- .../flow/GenericOutputBlackoilModule.cpp | 17 +---- .../flow/GenericOutputBlackoilModule.hpp | 2 - opm/simulators/flow/OutputBlackoilModule.hpp | 2 - .../flow/OutputCompositionalModule.hpp | 2 - opm/simulators/flow/TracerContainer.cpp | 64 +++++++++++++++++++ opm/simulators/flow/TracerContainer.hpp | 15 ++++- 6 files changed, 79 insertions(+), 23 deletions(-) diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 8e9f8590a..7eaedb64d 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -123,6 +123,7 @@ GenericOutputBlackoilModule(const EclipseState& eclState, , enableSaltPrecipitation_(enableSaltPrecipitation) , enableExtbo_(enableExtbo) , enableMICP_(enableMICP) + , tracerC_(eclState_) , local_data_valid_(false) { const auto& fp = eclState_.fieldProps(); @@ -789,8 +790,6 @@ doAllocBuffers(const unsigned bufferSize, const bool enablePCHysteresis, const bool enableNonWettingHysteresis, const bool enableWettingHysteresis, - const unsigned numTracers, - const std::vector& enableSolTracers, const unsigned numOutputNnc, std::map rstKeywords) { @@ -1220,19 +1219,7 @@ doAllocBuffers(const unsigned bufferSize, } // tracers - if (numTracers > 0) { - tracerC_.freeConcentrations_.resize(numTracers); - for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) - { - tracerC_.freeConcentrations_[tracerIdx].resize(bufferSize, 0.0); - } - tracerC_.solConcentrations_.resize(numTracers); - for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx) - { - if (enableSolTracers[tracerIdx]) - tracerC_.solConcentrations_[tracerIdx].resize(bufferSize, 0.0); - } - } + this->tracerC_.allocate(bufferSize); if (rstKeywords["RESIDUAL"] > 0) { rstKeywords["RESIDUAL"] = 0; diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index 059f5a093..f8ac70fbf 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -331,8 +331,6 @@ protected: const bool enablePCHysteresis = false, const bool enableNonWettingHysteresis = false, const bool enableWettingHysteresis = false, - const unsigned numTracers = 0, - const std::vector& enableSolTracers = {}, unsigned numOutputNnc = 0, std::map rstKeywords = {}); diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 71d44ea29..edcae0a96 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -188,8 +188,6 @@ public: problem.materialLawManager()->enablePCHysteresis(), problem.materialLawManager()->enableNonWettingHysteresis(), problem.materialLawManager()->enableWettingHysteresis(), - problem.tracerModel().numTracers(), - problem.tracerModel().enableSolTracers(), problem.eclWriter()->getOutputNnc().size()); } diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index 34b161164..d0dd19727 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -166,8 +166,6 @@ public: /* enablePCHysteresis = */ false, /* enableNonWettingHysteresis =*/ false, /* enableWettingHysteresis =*/ false, - /* numTracers = */ 0, - /* enableSoltracers =*/ {}, /* numOutputNnc =*/ 0, std::move(rstKeywords)); } diff --git a/opm/simulators/flow/TracerContainer.cpp b/opm/simulators/flow/TracerContainer.cpp index 6ceacbee1..9bf80e8c6 100644 --- a/opm/simulators/flow/TracerContainer.cpp +++ b/opm/simulators/flow/TracerContainer.cpp @@ -23,6 +23,70 @@ #include #include +#include + +#include +#include +#include + +#include + namespace Opm { +template +void TracerContainer:: +allocate(const unsigned bufferSize) +{ + const auto& tracers = eclState_.tracer(); + if (!tracers.empty()) { + allocated_ = true; + freeConcentrations_.resize(tracers.size()); + solConcentrations_.resize(tracers.size()); + std::for_each(tracers.begin(), tracers.end(), + [idx = 0, bufferSize, this](const auto& tracer) mutable + { + freeConcentrations_[idx].resize(bufferSize, 0.0); + if (((tracer.phase == Phase::GAS && FluidSystem::enableDissolvedGas()) || + (tracer.phase == Phase::OIL && FluidSystem::enableVaporizedOil())) && + (tracer.solution_concentration.has_value() || + tracer.solution_tvdp.has_value())) + { + solConcentrations_[idx].resize(bufferSize, 0.0); + } + ++idx; + }); + } +} + +template using FS = BlackOilFluidSystem; + +#define INSTANTIATE_TYPE(T) \ + template class TracerContainer>; + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +#define INSTANTIATE_COMP_THREEPHASE(NUM) \ + template using FS##NUM = GenericOilGasWaterFluidSystem; \ + template class TracerContainer>; + +#define INSTANTIATE_COMP_TWOPHASE(NUM) \ + template using GFS##NUM = GenericOilGasWaterFluidSystem; \ + template class TracerContainer>; + +#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/TracerContainer.hpp b/opm/simulators/flow/TracerContainer.hpp index 17f40f501..f34a808f6 100644 --- a/opm/simulators/flow/TracerContainer.hpp +++ b/opm/simulators/flow/TracerContainer.hpp @@ -30,6 +30,8 @@ namespace Opm { +class EclipseState; + template class TracerContainer { @@ -37,8 +39,17 @@ class TracerContainer using ScalarBuffer = std::vector; public: - std::vector freeConcentrations_; - std::vector solConcentrations_; + TracerContainer(const EclipseState& eclState) + : eclState_(eclState) + {} + + void allocate(const unsigned bufferSize); + + const EclipseState& eclState_; + + std::vector freeConcentrations_{}; + std::vector solConcentrations_{}; + bool allocated_{false}; }; } // namespace Opm From 7386655b6589fb602a48b654a932c21c8475b4dd Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 13 Feb 2025 15:02:36 +0100 Subject: [PATCH 28/36] move output of tracer data into TracerContainer --- .../flow/GenericOutputBlackoilModule.cpp | 36 +------------------ opm/simulators/flow/TracerContainer.cpp | 32 +++++++++++++++++ opm/simulators/flow/TracerContainer.hpp | 3 ++ 3 files changed, 36 insertions(+), 35 deletions(-) diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 7eaedb64d..07b20ca51 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -632,41 +632,7 @@ assignToSolution(data::Solution& sol) this->fipC_.outputRestart(sol); // Tracers - if (! this->tracerC_.freeConcentrations_.empty()) { - const auto& tracers = this->eclState_.tracer(); - for (auto tracerIdx = 0*tracers.size(); - tracerIdx < tracers.size(); ++tracerIdx) - { - sol.insert(tracers[tracerIdx].fname(), - UnitSystem::measure::identity, - std::move(this->tracerC_.freeConcentrations_[tracerIdx]), - data::TargetType::RESTART_TRACER_SOLUTION); - } - - // Put freeTracerConcentrations container into a valid state. Otherwise - // we'll move from vectors that have already been moved from if we - // get here and it's not a restart step. - this->tracerC_.freeConcentrations_.clear(); - } - if (! this->tracerC_.solConcentrations_.empty()) { - const auto& tracers = this->eclState_.tracer(); - for (auto tracerIdx = 0*tracers.size(); - tracerIdx < tracers.size(); ++tracerIdx) - { - if (this->tracerC_.solConcentrations_[tracerIdx].empty()) - continue; - - sol.insert(tracers[tracerIdx].sname(), - UnitSystem::measure::identity, - std::move(this->tracerC_.solConcentrations_[tracerIdx]), - data::TargetType::RESTART_TRACER_SOLUTION); - } - - // Put solTracerConcentrations container into a valid state. Otherwise - // we'll move from vectors that have already been moved from if we - // get here and it's not a restart step. - this->tracerC_.solConcentrations_.clear(); - } + this->tracerC_.outputRestart(sol); } template diff --git a/opm/simulators/flow/TracerContainer.cpp b/opm/simulators/flow/TracerContainer.cpp index 9bf80e8c6..7383547be 100644 --- a/opm/simulators/flow/TracerContainer.cpp +++ b/opm/simulators/flow/TracerContainer.cpp @@ -29,7 +29,10 @@ #include #include +#include + #include +#include namespace Opm { @@ -58,6 +61,35 @@ allocate(const unsigned bufferSize) } } +template +void TracerContainer:: +outputRestart(data::Solution& sol) +{ + if (!this->allocated_) { + return; + } + + const auto& tracers = this->eclState_.tracer(); + std::for_each(tracers.begin(), tracers.end(), + [idx = 0, &sol, this](const auto& tracer) mutable + { + sol.insert(tracer.fname(), + UnitSystem::measure::identity, + std::move(freeConcentrations_[idx]), + data::TargetType::RESTART_TRACER_SOLUTION); + + if (!solConcentrations_[idx].empty()) { + sol.insert(tracer.sname(), + UnitSystem::measure::identity, + std::move(solConcentrations_[idx]), + data::TargetType::RESTART_TRACER_SOLUTION); + } + ++idx; + }); + + this->allocated_ = false; +} + template using FS = BlackOilFluidSystem; #define INSTANTIATE_TYPE(T) \ diff --git a/opm/simulators/flow/TracerContainer.hpp b/opm/simulators/flow/TracerContainer.hpp index f34a808f6..5b0a7cef9 100644 --- a/opm/simulators/flow/TracerContainer.hpp +++ b/opm/simulators/flow/TracerContainer.hpp @@ -30,6 +30,7 @@ namespace Opm { +namespace data { class Solution; } class EclipseState; template @@ -45,6 +46,8 @@ public: void allocate(const unsigned bufferSize); + void outputRestart(data::Solution& sol); + const EclipseState& eclState_; std::vector freeConcentrations_{}; From 00fe5607203fdb59e8646471ccd0e2cc08e6654c Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 14 Feb 2025 09:14:56 +0100 Subject: [PATCH 29/36] move assignment of free tracer concentration into TracerContainer --- opm/simulators/flow/OutputBlackoilModule.hpp | 14 +++++--------- opm/simulators/flow/TracerContainer.cpp | 15 +++++++++++++++ opm/simulators/flow/TracerContainer.hpp | 6 ++++++ 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index edcae0a96..a2b702a3c 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -645,15 +645,11 @@ public: // tracers const auto& tracerModel = simulator_.problem().tracerModel(); - if (! this->tracerC_.freeConcentrations_.empty()) { - for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { - if (this->tracerC_.freeConcentrations_[tracerIdx].empty()) { - continue; - } - this->tracerC_.freeConcentrations_[tracerIdx][globalDofIdx] = - tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx); - } - } + this->tracerC_.assignFreeConcentrations(globalDofIdx, + [globalDofIdx, &tracerModel](const unsigned tracerIdx) + { return tracerModel.freeTracerConcentration(tracerIdx, + globalDofIdx); }); + if (! this->tracerC_.solConcentrations_.empty()) { for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { if (this->tracerC_.solConcentrations_[tracerIdx].empty()) { diff --git a/opm/simulators/flow/TracerContainer.cpp b/opm/simulators/flow/TracerContainer.cpp index 7383547be..9958c8a44 100644 --- a/opm/simulators/flow/TracerContainer.cpp +++ b/opm/simulators/flow/TracerContainer.cpp @@ -61,6 +61,21 @@ allocate(const unsigned bufferSize) } } +template +void TracerContainer:: +assignFreeConcentrations(const unsigned globalDofIdx, + const AssignFunction& concentration) +{ + std::for_each(freeConcentrations_.begin(), freeConcentrations_.end(), + [globalDofIdx, idx = 0, &concentration](auto& tracer) mutable + { + if (!tracer.empty()) { + tracer[globalDofIdx] = concentration(idx); + } + ++idx; + }); +} + template void TracerContainer:: outputRestart(data::Solution& sol) diff --git a/opm/simulators/flow/TracerContainer.hpp b/opm/simulators/flow/TracerContainer.hpp index 5b0a7cef9..87f3bd35b 100644 --- a/opm/simulators/flow/TracerContainer.hpp +++ b/opm/simulators/flow/TracerContainer.hpp @@ -26,6 +26,7 @@ #ifndef OPM_TRACER_CONTAINER_HPP #define OPM_TRACER_CONTAINER_HPP +#include #include namespace Opm { @@ -46,6 +47,11 @@ public: void allocate(const unsigned bufferSize); + using AssignFunction = std::function; + + void assignFreeConcentrations(const unsigned globalDofIdx, + const AssignFunction& concentration); + void outputRestart(data::Solution& sol); const EclipseState& eclState_; From dcc4aa30dd3319190a68e0da64add6ca24e03141 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 14 Feb 2025 09:14:56 +0100 Subject: [PATCH 30/36] move assignment of sol tracer concentration into TracerContainer --- opm/simulators/flow/OutputBlackoilModule.hpp | 15 ++++----------- opm/simulators/flow/TracerContainer.cpp | 15 +++++++++++++++ opm/simulators/flow/TracerContainer.hpp | 3 +++ 3 files changed, 22 insertions(+), 11 deletions(-) diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index a2b702a3c..44287fe71 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -649,17 +649,10 @@ public: [globalDofIdx, &tracerModel](const unsigned tracerIdx) { return tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx); }); - - if (! this->tracerC_.solConcentrations_.empty()) { - for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) { - if (this->tracerC_.solConcentrations_[tracerIdx].empty()) { - continue; - } - this->tracerC_.solConcentrations_[tracerIdx][globalDofIdx] = - tracerModel.solTracerConcentration(tracerIdx, globalDofIdx); - - } - } + this->tracerC_.assignSolConcentrations(globalDofIdx, + [globalDofIdx, &tracerModel](const unsigned tracerIdx) + { return tracerModel.solTracerConcentration(tracerIdx, + globalDofIdx); }); // output residual for ( int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx ) diff --git a/opm/simulators/flow/TracerContainer.cpp b/opm/simulators/flow/TracerContainer.cpp index 9958c8a44..9f73e758c 100644 --- a/opm/simulators/flow/TracerContainer.cpp +++ b/opm/simulators/flow/TracerContainer.cpp @@ -76,6 +76,21 @@ assignFreeConcentrations(const unsigned globalDofIdx, }); } +template +void TracerContainer:: +assignSolConcentrations(const unsigned globalDofIdx, + const AssignFunction& concentration) +{ + std::for_each(solConcentrations_.begin(), solConcentrations_.end(), + [globalDofIdx, idx = 0, &concentration](auto& tracer) mutable + { + if (!tracer.empty()) { + tracer[globalDofIdx] = concentration(idx); + } + ++idx; + }); +} + template void TracerContainer:: outputRestart(data::Solution& sol) diff --git a/opm/simulators/flow/TracerContainer.hpp b/opm/simulators/flow/TracerContainer.hpp index 87f3bd35b..d9c8b1461 100644 --- a/opm/simulators/flow/TracerContainer.hpp +++ b/opm/simulators/flow/TracerContainer.hpp @@ -52,6 +52,9 @@ public: void assignFreeConcentrations(const unsigned globalDofIdx, const AssignFunction& concentration); + void assignSolConcentrations(const unsigned globalDofIdx, + const AssignFunction& concentration); + void outputRestart(data::Solution& sol); const EclipseState& eclState_; From 2e5ad4f94b0c1e330a3356686b5070ebb23ff414 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 14 Feb 2025 09:18:57 +0100 Subject: [PATCH 31/36] TracerContainer: make data private --- opm/simulators/flow/TracerContainer.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/simulators/flow/TracerContainer.hpp b/opm/simulators/flow/TracerContainer.hpp index d9c8b1461..dfd8ff86a 100644 --- a/opm/simulators/flow/TracerContainer.hpp +++ b/opm/simulators/flow/TracerContainer.hpp @@ -57,6 +57,7 @@ public: void outputRestart(data::Solution& sol); +private: const EclipseState& eclState_; std::vector freeConcentrations_{}; From 652a5a4646fc71990abc2b33de33592f325524ae Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 14 Feb 2025 11:11:49 +0100 Subject: [PATCH 32/36] GenericOutputBlackoilModule: check schedule directly for vappars avoids a bool param in doAllocBuffers --- opm/simulators/flow/GenericOutputBlackoilModule.cpp | 3 ++- opm/simulators/flow/GenericOutputBlackoilModule.hpp | 1 - opm/simulators/flow/OutputBlackoilModule.hpp | 1 - opm/simulators/flow/OutputCompositionalModule.hpp | 1 - 4 files changed, 2 insertions(+), 4 deletions(-) diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 07b20ca51..96e0d20f8 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -752,7 +752,6 @@ doAllocBuffers(const unsigned bufferSize, const bool substep, const bool log, const bool isRestart, - const bool vapparsActive, const bool enablePCHysteresis, const bool enableNonWettingHysteresis, const bool enableWettingHysteresis, @@ -953,6 +952,8 @@ doAllocBuffers(const unsigned bufferSize, this->micpC_.allocate(bufferSize); } + const bool vapparsActive = schedule_[std::max(reportStepNum, 0u)].oilvap().getType() == + OilVaporizationProperties::OilVaporization::VAPPARS; if (vapparsActive) { soMax_.resize(bufferSize, 0.0); } diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index f8ac70fbf..79beec842 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -327,7 +327,6 @@ protected: const bool substep, const bool log, const bool isRestart, - const bool vapparsActive = false, const bool enablePCHysteresis = false, const bool enableNonWettingHysteresis = false, const bool enableWettingHysteresis = false, diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 44287fe71..d45eb6f28 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -184,7 +184,6 @@ public: substep, log, isRestart, - problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)), problem.materialLawManager()->enablePCHysteresis(), problem.materialLawManager()->enableNonWettingHysteresis(), problem.materialLawManager()->enableWettingHysteresis(), diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index d0dd19727..f94289742 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -162,7 +162,6 @@ public: this->compC_.allocate(bufferSize, rstKeywords); this->doAllocBuffers(bufferSize, reportStepNum, substep, log, isRestart, - /* vapparsActive =*/ false, /* enablePCHysteresis = */ false, /* enableNonWettingHysteresis =*/ false, /* enableWettingHysteresis =*/ false, From 41fbe1c839b736c51dd0858d2bf2c45ab4fdb00f Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 14 Feb 2025 11:30:27 +0100 Subject: [PATCH 33/36] doAllocBuffers: pass hysteresis config instead of 3 bools --- .../flow/GenericOutputBlackoilModule.cpp | 13 +++++------- .../flow/GenericOutputBlackoilModule.hpp | 5 ++--- opm/simulators/flow/OutputBlackoilModule.hpp | 20 +++++++++---------- .../flow/OutputCompositionalModule.hpp | 4 +--- 4 files changed, 17 insertions(+), 25 deletions(-) diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 96e0d20f8..e7bef59a8 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -27,9 +27,9 @@ #include +#include #include #include - #include #include @@ -55,7 +55,6 @@ #include #include #include -#include #include #include #include @@ -752,9 +751,7 @@ doAllocBuffers(const unsigned bufferSize, const bool substep, const bool log, const bool isRestart, - const bool enablePCHysteresis, - const bool enableNonWettingHysteresis, - const bool enableWettingHysteresis, + const EclHysteresisConfig* hysteresisConfig, const unsigned numOutputNnc, std::map rstKeywords) { @@ -958,7 +955,7 @@ doAllocBuffers(const unsigned bufferSize, soMax_.resize(bufferSize, 0.0); } - if (enableNonWettingHysteresis) { + if (hysteresisConfig && hysteresisConfig->enableNonWettingHysteresis()) { if (FluidSystem::phaseIsActive(oilPhaseIdx)){ if (FluidSystem::phaseIsActive(waterPhaseIdx)){ soMax_.resize(bufferSize, 0.0); @@ -970,7 +967,7 @@ doAllocBuffers(const unsigned bufferSize, //TODO add support for gas-water } } - if (enableWettingHysteresis) { + if (hysteresisConfig && hysteresisConfig->enableWettingHysteresis()) { if (FluidSystem::phaseIsActive(oilPhaseIdx)){ if (FluidSystem::phaseIsActive(waterPhaseIdx)){ swMax_.resize(bufferSize, 0.0); @@ -982,7 +979,7 @@ doAllocBuffers(const unsigned bufferSize, //TODO add support for gas-water } } - if (enablePCHysteresis) { + if (hysteresisConfig && hysteresisConfig->enablePCHysteresis()) { if (FluidSystem::phaseIsActive(oilPhaseIdx)){ if (FluidSystem::phaseIsActive(waterPhaseIdx)){ swmin_.resize(bufferSize, 0.0); diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index 79beec842..b73cc4144 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -63,6 +63,7 @@ struct ForceDisableResvFluidInPlaceOutput { static constexpr bool value = false; namespace Opm { namespace data { class Solution; } +class EclHysteresisConfig; class EclipseState; class Schedule; class SummaryConfig; @@ -327,9 +328,7 @@ protected: const bool substep, const bool log, const bool isRestart, - const bool enablePCHysteresis = false, - const bool enableNonWettingHysteresis = false, - const bool enableWettingHysteresis = false, + const EclHysteresisConfig* hysteresisConfig, unsigned numOutputNnc = 0, std::map rstKeywords = {}); diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index d45eb6f28..a6fcc58f2 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -184,9 +184,7 @@ public: substep, log, isRestart, - problem.materialLawManager()->enablePCHysteresis(), - problem.materialLawManager()->enableNonWettingHysteresis(), - problem.materialLawManager()->enableWettingHysteresis(), + &problem.materialLawManager()->hysteresisConfig(), problem.eclWriter()->getOutputNnc().size()); } @@ -514,7 +512,7 @@ public: const auto& matLawManager = problem.materialLawManager(); if (matLawManager->enableHysteresis()) { - if (FluidSystem::phaseIsActive(oilPhaseIdx) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) { Scalar somax; Scalar swmax; @@ -522,7 +520,7 @@ public: matLawManager->oilWaterHysteresisParams( somax, swmax, swmin, globalDofIdx); - + if (matLawManager->enableNonWettingHysteresis()) { if (!this->soMax_.empty()) { this->soMax_[globalDofIdx] = somax; @@ -540,14 +538,14 @@ public: } } - if (FluidSystem::phaseIsActive(oilPhaseIdx) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) { Scalar sgmax; Scalar shmax; Scalar somin; matLawManager->gasOilHysteresisParams( sgmax, shmax, somin, globalDofIdx); - + if (matLawManager->enableNonWettingHysteresis()) { if (!this->sgmax_.empty()) { this->sgmax_[globalDofIdx] = sgmax; @@ -565,7 +563,7 @@ public: } } } else { - + if (!this->soMax_.empty()) this->soMax_[globalDofIdx] = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx)); @@ -1162,7 +1160,7 @@ public: if (simulator.problem().materialLawManager()->enableHysteresis()) { auto matLawManager = simulator.problem().materialLawManager(); - if (FluidSystem::phaseIsActive(oilPhaseIdx) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(waterPhaseIdx)) { Scalar somax = 2.0; Scalar swmax = -2.0; @@ -1186,7 +1184,7 @@ public: matLawManager->setOilWaterHysteresisParams( somax, swmax, swmin, elemIdx); } - if (FluidSystem::phaseIsActive(oilPhaseIdx) + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) { Scalar sgmax = 2.0; Scalar shmax = -2.0; @@ -1524,7 +1522,7 @@ private: { this->updateCO2InGas(globalDofIdx, pv, intQuants); } - + if (this->fipC_.hasCo2InWater() && (FluidSystem::phaseIsActive(waterPhaseIdx) || FluidSystem::phaseIsActive(oilPhaseIdx))) diff --git a/opm/simulators/flow/OutputCompositionalModule.hpp b/opm/simulators/flow/OutputCompositionalModule.hpp index f94289742..97606d9d8 100644 --- a/opm/simulators/flow/OutputCompositionalModule.hpp +++ b/opm/simulators/flow/OutputCompositionalModule.hpp @@ -162,9 +162,7 @@ public: this->compC_.allocate(bufferSize, rstKeywords); this->doAllocBuffers(bufferSize, reportStepNum, substep, log, isRestart, - /* enablePCHysteresis = */ false, - /* enableNonWettingHysteresis =*/ false, - /* enableWettingHysteresis =*/ false, + /* hysteresisConfig = */ nullptr, /* numOutputNnc =*/ 0, std::move(rstKeywords)); } From cc68f7246a73541e37a32d05a860d6616303d359 Mon Sep 17 00:00:00 2001 From: Vegard Kippe Date: Fri, 17 Jan 2025 14:49:40 +0100 Subject: [PATCH 34/36] Support instantaneous flow rates in extended network (WEFAC and GEFAC item 3) --- .../utils/FullySupportedFlowKeywords.cpp | 12 +++++ .../utils/PartiallySupportedFlowKeywords.cpp | 12 ----- opm/simulators/wells/GroupState.cpp | 30 +++++++++++ opm/simulators/wells/GroupState.hpp | 7 +++ opm/simulators/wells/WellGroupHelpers.cpp | 51 +++++++++++-------- opm/simulators/wells/WellGroupHelpers.hpp | 3 +- 6 files changed, 80 insertions(+), 35 deletions(-) diff --git a/opm/simulators/utils/FullySupportedFlowKeywords.cpp b/opm/simulators/utils/FullySupportedFlowKeywords.cpp index ea13034b6..b867024b9 100644 --- a/opm/simulators/utils/FullySupportedFlowKeywords.cpp +++ b/opm/simulators/utils/FullySupportedFlowKeywords.cpp @@ -32,6 +32,12 @@ const SupportedKeywordItems& fullySupported() { static const SupportedKeywordItems fully_supported_keywords_strings = { + { + "GEFAC", + { + {3,{true, is_bool_convertible {}, "GEFAC(GRPNETWK): String value must be convertible to bool."}}, // USE_GEFAC_IN_NETWORK + }, + }, { "NEXTSTEP", { @@ -44,6 +50,12 @@ fullySupported() {3,{true, allow_values {"ORAT", "WRAT", "GRAT", "LRAT", "RESV", "BHP"}, "WCONHIST(TARGET): should be set to ORAT/WRAT/GRAT/LRAT/RESV or BHP"}}, // CMODE }, }, + { + "WEFAC", + { + {3,{true, is_bool_convertible {}, "WEFAC(WELNETWK): String value must be convertible to bool."}}, // USE_WEFAC_IN_NETWORK + }, + }, }; return fully_supported_keywords_strings; diff --git a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp index 2dada5013..3e5fc99c0 100644 --- a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp +++ b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp @@ -111,12 +111,6 @@ partiallySupported() {8,{true, allow_values {"NO"}, "GECON(ENDRUN): End run not implemented"}}, }, }, - { - "GEFAC", - { - {3,{true, allow_values {"YES"}, "GEFAC(GRPNETWK): Extended Network Model efficiency NO option not implemented"}}, // TRANSFER_EXT_NET - }, - }, { "GRIDOPTS", { @@ -273,12 +267,6 @@ partiallySupported() {5,{true, allow_values {"NO"}, "WAGHYSTR(WATER_MODEL): only the NO option is supported – will STOP"}}, // WATER_MODEL }, }, - { - "WEFAC", - { - {3,{true, allow_values {"YES"}, "WEFAC(WELNETWK): only the YES option is supported"}}, // EXTENDED_NETWORK_OPT - }, - }, { "WELSPECS", { diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index f1a78796f..ab25a065e 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -41,6 +41,7 @@ GroupState GroupState::serializationTestObject() { GroupState result(3); result.m_production_rates = {{"test1", {1.0, 2.0}}}; + result.m_network_production_rates={{"test1", {1.0, 20}}}; result.production_controls = {{"test2", Group::ProductionCMode::LRAT}}; result.prod_red_rates = {{"test3", {3.0, 4.0, 5.0}}}; result.inj_red_rates = {{"test4", {6.0, 7.0}}}; @@ -61,6 +62,7 @@ template bool GroupState::operator==(const GroupState& other) const { return this->m_production_rates == other.m_production_rates && + this->m_network_production_rates == other.m_network_production_rates && this->production_controls == other.production_controls && this->prod_red_rates == other.prod_red_rates && this->inj_red_rates == other.inj_red_rates && @@ -83,6 +85,13 @@ bool GroupState::has_production_rates(const std::string& gname) const return (group_iter != this->m_production_rates.end()); } +template +bool GroupState::has_network_production_rates(const std::string& gname) const +{ + auto group_iter = this->m_network_production_rates.find(gname); + return (group_iter != this->m_network_production_rates.end()); +} + template void GroupState::update_production_rates(const std::string& gname, const std::vector& rates) @@ -93,6 +102,16 @@ void GroupState::update_production_rates(const std::string& gname, this->m_production_rates[gname] = rates; } +template +void GroupState::update_network_production_rates(const std::string& gname, + const std::vector& rates) +{ + if (rates.size() != this->num_phases) + throw std::logic_error("Wrong number of phases"); + + this->m_network_production_rates[gname] = rates; +} + template const std::vector& GroupState::production_rates(const std::string& gname) const @@ -104,6 +123,17 @@ GroupState::production_rates(const std::string& gname) const return group_iter->second; } +template +const std::vector& +GroupState::network_production_rates(const std::string& gname) const +{ + auto group_iter = this->m_network_production_rates.find(gname); + if (group_iter == this->m_network_production_rates.end()) + throw std::logic_error("No such group"); + + return group_iter->second; +} + //------------------------------------------------------------------------- template diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index 06bf1466f..4a801eedd 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -50,9 +50,13 @@ public: bool operator==(const GroupState& other) const; bool has_production_rates(const std::string& gname) const; + bool has_network_production_rates(const std::string& gname) const; void update_production_rates(const std::string& gname, const std::vector& rates); + void update_network_production_rates(const std::string& gname, + const std::vector& rates); const std::vector& production_rates(const std::string& gname) const; + const std::vector& network_production_rates(const std::string& gname) const; void update_well_group_thp(const std::string& gname, const double& thp); Scalar well_group_thp(const std::string& gname) const; @@ -130,6 +134,7 @@ public: auto forAllGroupData = [&](auto& func) { iterateContainer(m_production_rates, func); + iterateContainer(m_network_production_rates, func); iterateContainer(prod_red_rates, func); iterateContainer(inj_red_rates, func); iterateContainer(inj_resv_rates, func); @@ -187,6 +192,7 @@ public: { serializer(num_phases); serializer(m_production_rates); + serializer(m_network_production_rates); serializer(production_controls); serializer(group_thp); serializer(prod_red_rates); @@ -205,6 +211,7 @@ public: private: std::size_t num_phases{}; std::map> m_production_rates; + std::map> m_network_production_rates; std::map production_controls; std::map> prod_red_rates; std::map> inj_red_rates; diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 86e9ca791..c0998a047 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -90,14 +90,15 @@ namespace Opm { const Opm::WellState& wellState, const int reportStepIdx, const int phasePos, - const bool injector) + const bool injector, + const bool network) { Scalar rate = 0.0; for (const std::string& groupName : group.groups()) { const auto& groupTmp = schedule.getGroup(groupName, reportStepIdx); - const auto& gefac = groupTmp.getGroupEfficiencyFactor(); - rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector); + const auto& gefac = groupTmp.getGroupEfficiencyFactor(network); + rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector, network); } // only sum satelite production once @@ -126,8 +127,7 @@ namespace Opm { if (wellEcl.getStatus() == Opm::Well::Status::SHUT) continue; - const Scalar factor = wellEcl.getEfficiencyFactor() * - wellState[wellEcl.name()].efficiency_scaling_factor; + Scalar factor = wellEcl.getEfficiencyFactor(network) * wellState[wellEcl.name()].efficiency_scaling_factor; const auto& ws = wellState.well(well_index.value()); if (res_rates) { const auto& well_rates = ws.reservoir_rates; @@ -742,6 +742,17 @@ updateGroupProductionRates(const Group& group, rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false); } group_state.update_production_rates(group.name(), rates); + + const auto& network = schedule[reportStepIdx].network(); + if (network.active()) { + std::vector network_rates = rates; + if (network.needs_instantaneous_rates(schedule, reportStepIdx)) { + for (int phase = 0; phase < np; ++phase) { + network_rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false, /*network*/ true); + } + } + group_state.update_network_production_rates(group.name(), network_rates); + } } template @@ -927,24 +938,16 @@ computeNetworkPressures(const Network::ExtNetwork& network, node_inflows[node] = zero_rates; continue; } + node_inflows[node] = group_state.network_production_rates(node); - node_inflows[node] = group_state.production_rates(node); // Add the ALQ amounts to the gas rates if requested. if (network.node(node).add_gas_lift_gas()) { const auto& group = schedule.getGroup(node, report_time_step); for (const std::string& wellname : group.wells()) { const Well& well = schedule.getWell(wellname, report_time_step); + if (well.isInjector() || !well_state.isOpen(wellname)) continue; - // Here we use the efficiency unconditionally, but if WEFAC item 3 - // for the well is false (it defaults to true) then we should NOT use - // the efficiency factor. Fixing this requires not only changing the - // code here, but also: - // - Adding a member to the well for this flag, and setting it in Schedule::handleWEFAC(). - // - Making the wells' maximum flows (i.e. not time-averaged by using a efficiency factor) - // available and using those (for wells with WEFAC(3) true only) when accumulating group - // rates, but ONLY for network calculations. - const Scalar efficiency = well.getEfficiencyFactor() * - well_state.getGlobalEfficiencyScalingFactor(wellname); + const Scalar efficiency = well.getEfficiencyFactor(/*network*/ true) * well_state.getGlobalEfficiencyScalingFactor(wellname); node_inflows[node][BlackoilPhases::Vapour] += well_state.getALQ(wellname) * efficiency; } } @@ -962,13 +965,17 @@ computeNetworkPressures(const Network::ExtNetwork& network, // Add downbranch rates to upbranch. std::vector& up = node_inflows[(*upbranch).uptree_node()]; const std::vector& down = node_inflows[node]; + // @TODO@ Also support NEFAC (for nodes that do not correspond to groups) + double efficiency = 1.0; + if (schedule.hasGroup(node, report_time_step)) { + efficiency = schedule.getGroup(node, report_time_step).getGroupEfficiencyFactor(/*network*/ true); + } if (up.empty()) { - up = down; - } else { - assert (up.size() == down.size()); - for (std::size_t ii = 0; ii < up.size(); ++ii) { - up[ii] += down[ii]; - } + up = std::vector(down.size(), 0.0); + } + assert (up.size() == down.size()); + for (std::size_t ii = 0; ii < up.size(); ++ii) { + up[ii] += (efficiency*down[ii]); } } } diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 62bd15b28..876bf2420 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -56,7 +56,8 @@ public: const Opm::WellState& wellState, const int reportStepIdx, const int phasePos, - const bool injector); + const bool injector, + const bool network = false); static Scalar satelliteProduction(const ScheduleState& sched, const std::vector& groups, From 44a87f8e658c1ecc4b82d57ad55c828df5214d82 Mon Sep 17 00:00:00 2001 From: Vegard Kippe Date: Mon, 20 Jan 2025 13:24:01 +0100 Subject: [PATCH 35/36] Support NEFAC and clean up implementation --- .../utils/FullySupportedFlowKeywords.cpp | 6 ++++ .../utils/UnsupportedFlowKeywords.cpp | 1 - .../wells/BlackoilWellModelGeneric.cpp | 4 +++ opm/simulators/wells/GroupState.cpp | 20 ++++++------- opm/simulators/wells/GroupState.hpp | 12 ++++---- opm/simulators/wells/WellGroupHelpers.cpp | 30 ++++++++++++------- opm/simulators/wells/WellGroupHelpers.hpp | 6 ++++ 7 files changed, 52 insertions(+), 27 deletions(-) diff --git a/opm/simulators/utils/FullySupportedFlowKeywords.cpp b/opm/simulators/utils/FullySupportedFlowKeywords.cpp index b867024b9..669635894 100644 --- a/opm/simulators/utils/FullySupportedFlowKeywords.cpp +++ b/opm/simulators/utils/FullySupportedFlowKeywords.cpp @@ -78,6 +78,12 @@ const SupportedKeywordItems& fullySupported() { static const SupportedKeywordItems fully_supported_keywords_double = { + { + "NEFAC", + { + {2,{true, [](double x) { return x > 0 && x <= 1.0; }, "NEFAC(EFF_FACTOR: Efficiency must be in the range (0,1]"}}, // NETWORK_EFF_FACTOR + }, + }, { "WPIMULT", { diff --git a/opm/simulators/utils/UnsupportedFlowKeywords.cpp b/opm/simulators/utils/UnsupportedFlowKeywords.cpp index da61bbd8b..ce6196ce3 100644 --- a/opm/simulators/utils/UnsupportedFlowKeywords.cpp +++ b/opm/simulators/utils/UnsupportedFlowKeywords.cpp @@ -397,7 +397,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords() {"NARROW", {true, std::nullopt}}, {"NCONSUMP", {true, std::nullopt}}, {"NCOMPS", {false, std::nullopt}}, - {"NEFAC", {true, std::nullopt}}, {"NETCOMPA", {true, std::nullopt}}, {"NEXT", {false, std::nullopt}}, {"NEXTSTPL", {true, std::nullopt}}, diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index ca69376de..d1c3b48f6 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1373,6 +1373,10 @@ updateAndCommunicateGroupData(const int reportStepIdx, reportStepIdx, well_state_nupcol, this->groupState()); + WellGroupHelpers::updateNetworkLeafNodeProductionRates(schedule(), + reportStepIdx, + well_state_nupcol, + this->groupState()); WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index ab25a065e..0a3f85b16 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -41,7 +41,7 @@ GroupState GroupState::serializationTestObject() { GroupState result(3); result.m_production_rates = {{"test1", {1.0, 2.0}}}; - result.m_network_production_rates={{"test1", {1.0, 20}}}; + result.m_network_leaf_node_production_rates={{"test1", {1.0, 20}}}; result.production_controls = {{"test2", Group::ProductionCMode::LRAT}}; result.prod_red_rates = {{"test3", {3.0, 4.0, 5.0}}}; result.inj_red_rates = {{"test4", {6.0, 7.0}}}; @@ -62,7 +62,7 @@ template bool GroupState::operator==(const GroupState& other) const { return this->m_production_rates == other.m_production_rates && - this->m_network_production_rates == other.m_network_production_rates && + this->m_network_leaf_node_production_rates == other.m_network_leaf_node_production_rates && this->production_controls == other.production_controls && this->prod_red_rates == other.prod_red_rates && this->inj_red_rates == other.inj_red_rates && @@ -86,10 +86,10 @@ bool GroupState::has_production_rates(const std::string& gname) const } template -bool GroupState::has_network_production_rates(const std::string& gname) const +bool GroupState::has_network_leaf_node_production_rates(const std::string& gname) const { - auto group_iter = this->m_network_production_rates.find(gname); - return (group_iter != this->m_network_production_rates.end()); + auto group_iter = this->m_network_leaf_node_production_rates.find(gname); + return (group_iter != this->m_network_leaf_node_production_rates.end()); } template @@ -103,13 +103,13 @@ void GroupState::update_production_rates(const std::string& gname, } template -void GroupState::update_network_production_rates(const std::string& gname, +void GroupState::update_network_leaf_node_production_rates(const std::string& gname, const std::vector& rates) { if (rates.size() != this->num_phases) throw std::logic_error("Wrong number of phases"); - this->m_network_production_rates[gname] = rates; + this->m_network_leaf_node_production_rates[gname] = rates; } template @@ -125,10 +125,10 @@ GroupState::production_rates(const std::string& gname) const template const std::vector& -GroupState::network_production_rates(const std::string& gname) const +GroupState::network_leaf_node_production_rates(const std::string& gname) const { - auto group_iter = this->m_network_production_rates.find(gname); - if (group_iter == this->m_network_production_rates.end()) + auto group_iter = this->m_network_leaf_node_production_rates.find(gname); + if (group_iter == this->m_network_leaf_node_production_rates.end()) throw std::logic_error("No such group"); return group_iter->second; diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index 4a801eedd..b05f21e7a 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -50,13 +50,13 @@ public: bool operator==(const GroupState& other) const; bool has_production_rates(const std::string& gname) const; - bool has_network_production_rates(const std::string& gname) const; + bool has_network_leaf_node_production_rates(const std::string& gname) const; void update_production_rates(const std::string& gname, const std::vector& rates); - void update_network_production_rates(const std::string& gname, + void update_network_leaf_node_production_rates(const std::string& gname, const std::vector& rates); const std::vector& production_rates(const std::string& gname) const; - const std::vector& network_production_rates(const std::string& gname) const; + const std::vector& network_leaf_node_production_rates(const std::string& gname) const; void update_well_group_thp(const std::string& gname, const double& thp); Scalar well_group_thp(const std::string& gname) const; @@ -134,7 +134,7 @@ public: auto forAllGroupData = [&](auto& func) { iterateContainer(m_production_rates, func); - iterateContainer(m_network_production_rates, func); + iterateContainer(m_network_leaf_node_production_rates, func); iterateContainer(prod_red_rates, func); iterateContainer(inj_red_rates, func); iterateContainer(inj_resv_rates, func); @@ -192,7 +192,7 @@ public: { serializer(num_phases); serializer(m_production_rates); - serializer(m_network_production_rates); + serializer(m_network_leaf_node_production_rates); serializer(production_controls); serializer(group_thp); serializer(prod_red_rates); @@ -211,7 +211,7 @@ public: private: std::size_t num_phases{}; std::map> m_production_rates; - std::map> m_network_production_rates; + std::map> m_network_leaf_node_production_rates; std::map production_controls; std::map> prod_red_rates; std::map> inj_red_rates; diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index c0998a047..8fc85cc5b 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -742,16 +742,29 @@ updateGroupProductionRates(const Group& group, rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false); } group_state.update_production_rates(group.name(), rates); +} +template +void WellGroupHelpers:: +updateNetworkLeafNodeProductionRates(const Schedule& schedule, + const int reportStepIdx, + const WellState& wellState, + GroupState& group_state) +{ const auto& network = schedule[reportStepIdx].network(); if (network.active()) { - std::vector network_rates = rates; - if (network.needs_instantaneous_rates(schedule, reportStepIdx)) { - for (int phase = 0; phase < np; ++phase) { - network_rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false, /*network*/ true); + const int np = wellState.numPhases(); + for (const auto& group_name : network.leaf_nodes()) { + assert(schedule[reportStepIdx].groups.has(group_name)); + const auto& group = schedule[reportStepIdx].groups.get(group_name); + std::vector network_rates(np, 0.0); + if (group.numWells() > 0) { + for (int phase = 0; phase < np; ++phase) { + network_rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false, /*network*/ true); + } } + group_state.update_network_leaf_node_production_rates(group_name, network_rates); } - group_state.update_network_production_rates(group.name(), network_rates); } } @@ -938,8 +951,8 @@ computeNetworkPressures(const Network::ExtNetwork& network, node_inflows[node] = zero_rates; continue; } - node_inflows[node] = group_state.network_production_rates(node); + node_inflows[node] = group_state.network_leaf_node_production_rates(node); // Add the ALQ amounts to the gas rates if requested. if (network.node(node).add_gas_lift_gas()) { const auto& group = schedule.getGroup(node, report_time_step); @@ -966,10 +979,7 @@ computeNetworkPressures(const Network::ExtNetwork& network, std::vector& up = node_inflows[(*upbranch).uptree_node()]; const std::vector& down = node_inflows[node]; // @TODO@ Also support NEFAC (for nodes that do not correspond to groups) - double efficiency = 1.0; - if (schedule.hasGroup(node, report_time_step)) { - efficiency = schedule.getGroup(node, report_time_step).getGroupEfficiencyFactor(/*network*/ true); - } + const double efficiency = network.node(node).efficiency(); if (up.empty()) { up = std::vector(down.size(), 0.0); } diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 876bf2420..c5d0e23ca 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -200,6 +200,12 @@ public: const WellState& wellState, GroupState& group_state); + static void updateNetworkLeafNodeProductionRates(const Schedule& schedule, + const int reportStepIdx, + const WellState& wellState, + GroupState& group_state); + + static void updateWellRatesFromGroupTargetScale(const Scalar scale, const Group& group, const Schedule& schedule, From 1913186383a28768064ffc858160e8d491b37e3b Mon Sep 17 00:00:00 2001 From: Vegard Kippe Date: Fri, 14 Feb 2025 10:51:23 +0100 Subject: [PATCH 36/36] Minor updates post review --- opm/simulators/wells/GroupState.cpp | 7 ------- opm/simulators/wells/GroupState.hpp | 1 - opm/simulators/wells/WellGroupHelpers.cpp | 8 ++++---- 3 files changed, 4 insertions(+), 12 deletions(-) diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index 0a3f85b16..73643c69a 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -85,13 +85,6 @@ bool GroupState::has_production_rates(const std::string& gname) const return (group_iter != this->m_production_rates.end()); } -template -bool GroupState::has_network_leaf_node_production_rates(const std::string& gname) const -{ - auto group_iter = this->m_network_leaf_node_production_rates.find(gname); - return (group_iter != this->m_network_leaf_node_production_rates.end()); -} - template void GroupState::update_production_rates(const std::string& gname, const std::vector& rates) diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index b05f21e7a..41cc270bd 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -50,7 +50,6 @@ public: bool operator==(const GroupState& other) const; bool has_production_rates(const std::string& gname) const; - bool has_network_leaf_node_production_rates(const std::string& gname) const; void update_production_rates(const std::string& gname, const std::vector& rates); void update_network_leaf_node_production_rates(const std::string& gname, diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 8fc85cc5b..74425033d 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -127,7 +127,7 @@ namespace Opm { if (wellEcl.getStatus() == Opm::Well::Status::SHUT) continue; - Scalar factor = wellEcl.getEfficiencyFactor(network) * wellState[wellEcl.name()].efficiency_scaling_factor; + const Scalar factor = wellEcl.getEfficiencyFactor(network) * wellState[wellEcl.name()].efficiency_scaling_factor; const auto& ws = wellState.well(well_index.value()); if (res_rates) { const auto& well_rates = ws.reservoir_rates; @@ -978,14 +978,14 @@ computeNetworkPressures(const Network::ExtNetwork& network, // Add downbranch rates to upbranch. std::vector& up = node_inflows[(*upbranch).uptree_node()]; const std::vector& down = node_inflows[node]; - // @TODO@ Also support NEFAC (for nodes that do not correspond to groups) - const double efficiency = network.node(node).efficiency(); + // We now also support NEFAC + const Scalar efficiency = network.node(node).efficiency(); if (up.empty()) { up = std::vector(down.size(), 0.0); } assert (up.size() == down.size()); for (std::size_t ii = 0; ii < up.size(); ++ii) { - up[ii] += (efficiency*down[ii]); + up[ii] += efficiency*down[ii]; } } }