diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 18bd866a83..7a0f8f3410 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -33,6 +33,7 @@ list (APPEND MAIN_SOURCE_FILES ebos/eclgenericvanguard.cc ebos/eclgenericwriter.cc ebos/eclinterregflows.cc + ebos/eclmixingratecontrols.cc ebos/eclsolutioncontainers.cc ebos/ecltransmissibility.cc ebos/equil/equilibrationhelpers.cc @@ -386,6 +387,7 @@ list (APPEND PUBLIC_HEADER_FILES ebos/eclgenericwriter.hh ebos/eclgenericwriter_impl.hh ebos/eclinterregflows.hh + ebos/eclmixingratecontrols.hh ebos/eclmpiserializer.hh ebos/eclnewtonmethod.hh ebos/ecloutputblackoilmodule.hh diff --git a/ebos/eclgenericproblem.hh b/ebos/eclgenericproblem.hh index 3d607f7957..2f80bd6153 100644 --- a/ebos/eclgenericproblem.hh +++ b/ebos/eclgenericproblem.hh @@ -28,6 +28,7 @@ #ifndef EWOMS_GENERIC_ECL_PROBLEM_HH #define EWOMS_GENERIC_ECL_PROBLEM_HH +#include #include #include @@ -275,8 +276,8 @@ public: { return maxFails_; } bool vapparsActive(int episodeIdx) const; - - int numPressurePointsEquil() const + + int numPressurePointsEquil() const { return numPressurePointsEquil_; } bool operator==(const EclGenericProblem& rhs) const; @@ -291,22 +292,11 @@ public: serializer(overburdenPressure_); serializer(solventSaturation_); serializer(micp_); - serializer(lastRv_); - serializer(maxDRv_); - serializer(convectiveDrs_); - serializer(lastRs_); - serializer(maxDRs_); - serializer(dRsDtOnlyFreeGas_); + serializer(mixControls_); } protected: - bool drsdtActive_(int episodeIdx) const; - bool drvdtActive_(int episodeIdx) const; - bool drsdtConvective_(int episodeIdx) const; - void initFluidSystem_(); - void initDRSDT_(size_t numDof, - int episodeIdx); /*! * \brief Always returns true. The ecl output writer takes care of the rest @@ -383,13 +373,7 @@ protected: std::vector solventSaturation_; MICPSolutionContainer micp_; - std::vector lastRv_; - std::vector maxDRv_; - - std::vector convectiveDrs_; - std::vector lastRs_; - std::vector maxDRs_; - std::vector dRsDtOnlyFreeGas_; // apply the DRSDT rate limit only to cells that exhibit free gas + EclMixingRateControls mixControls_; // time stepping parameters bool enableTuning_; diff --git a/ebos/eclgenericproblem_impl.hh b/ebos/eclgenericproblem_impl.hh index c523cb75af..a897fb3913 100644 --- a/ebos/eclgenericproblem_impl.hh +++ b/ebos/eclgenericproblem_impl.hh @@ -88,6 +88,7 @@ EclGenericProblem(const EclipseState& eclState, : eclState_(eclState) , schedule_(schedule) , gridView_(gridView) + , mixControls_(schedule) { } @@ -106,12 +107,7 @@ serializationTestObject(const EclipseState& eclState, result.solventSaturation_ = {15.0}; result.polymer_ = PolymerSolutionContainer::serializationTestObject(); result.micp_ = MICPSolutionContainer::serializationTestObject(); - result.lastRv_ = {21.0}; - result.maxDRv_ = {22.0, 23.0}; - result.convectiveDrs_ = {24.0, 25.0, 26.0}; - result.lastRs_ = {27.0}; - result.maxDRs_ = {28.0}; - result.dRsDtOnlyFreeGas_ = {false, true}; + result.mixControls_ = EclMixingRateControls::serializationTestObject(schedule); return result; } @@ -444,36 +440,6 @@ vapparsActive(int episodeIdx) const return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS); } -template -bool EclGenericProblem:: -drsdtActive_(int episodeIdx) const -{ - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drsdtActive() && bothOilGasActive); -} - -template -bool EclGenericProblem:: -drvdtActive_(int episodeIdx) const -{ - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drvdtActive() && bothOilGasActive); -} - -template -bool EclGenericProblem:: -drsdtConvective_(int episodeIdx) const -{ - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drsdtConvective() && bothOilGasActive); -} - template bool EclGenericProblem:: beginEpisode_(bool enableExperiments, @@ -539,16 +505,7 @@ beginTimeStep_(bool enableExperiments, } // update explicit quantities between timesteps. - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - if (drsdtActive_(episodeIdx)) - // DRSDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) - maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*timeStepSize; - - if (drvdtActive_(episodeIdx)) - // DRVDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) - maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*timeStepSize; + this->mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize); } template @@ -615,7 +572,7 @@ readBlackoilExtentionsInitialConditions_(size_t numDof, } else { micp_.calciteConcentration.resize(numDof, 0.0); } -} + } } @@ -663,15 +620,8 @@ template Scalar EclGenericProblem:: drsdtcon(unsigned elemIdx, int episodeIdx) const { - if (convectiveDrs_.empty()) - return 0; - - // The episode index is set to -1 in the initialization phase. - // Output drsdt value for index 0 - episodeIdx = std::max(episodeIdx, 0); - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - int pvtRegionIdx = pvtRegionIndex(elemIdx); - return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*convectiveDrs_[elemIdx]; + return this->mixControls_.drsdtcon(elemIdx, episodeIdx, + this->pvtRegionIndex(elemIdx)); } template @@ -802,48 +752,17 @@ maxPolymerAdsorption(unsigned elemIdx) const return polymer_.maxAdsorption[elemIdx]; } -template -void EclGenericProblem:: -initDRSDT_(size_t numDof, - int episodeIdx) -{ - // deal with DRSDT - unsigned ntpvt = eclState_.runspec().tabdims().getNumPVTTables(); - //TODO We may want to only allocate these properties only if active. - //But since they may be activated at later time we need some more - //intrastructure to handle it - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - maxDRv_.resize(ntpvt, 1e30); - lastRv_.resize(numDof, 0.0); - maxDRs_.resize(ntpvt, 1e30); - dRsDtOnlyFreeGas_.resize(ntpvt, false); - lastRs_.resize(numDof, 0.0); - maxDRv_.resize(ntpvt, 1e30); - lastRv_.resize(numDof, 0.0); - maxOilSaturation_.resize(numDof, 0.0); - if (drsdtConvective_(episodeIdx)) { - convectiveDrs_.resize(numDof, 1.0); - } - } -} - template bool EclGenericProblem:: operator==(const EclGenericProblem& rhs) const { - return this->maxOilSaturation_ == rhs.maxOilSaturation_ && - this->maxWaterSaturation_ == rhs.maxWaterSaturation_ && + return this->maxWaterSaturation_ == rhs.maxWaterSaturation_ && this->minOilPressure_ == rhs.minOilPressure_ && this->overburdenPressure_ == rhs.overburdenPressure_ && this->solventSaturation_ == rhs.solventSaturation_ && this->polymer_ == rhs.polymer_ && this->micp_ == rhs.micp_ && - this->lastRv_ == rhs.lastRv_ && - this->maxDRv_ == rhs.maxDRv_ && - this->convectiveDrs_ == rhs.convectiveDrs_ && - this->lastRs_ == rhs.lastRs_ && - this->maxDRs_ == rhs.maxDRs_ && - this->dRsDtOnlyFreeGas_ == rhs.dRsDtOnlyFreeGas_; + this->mixControls_ == rhs.mixControls_; } } // namespace Opm diff --git a/ebos/eclmixingratecontrols.cc b/ebos/eclmixingratecontrols.cc new file mode 100644 index 0000000000..94d3526c6e --- /dev/null +++ b/ebos/eclmixingratecontrols.cc @@ -0,0 +1,299 @@ +// -*- 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::EclProblem + */ + +#include +#include + +#include + +#include +#include + +namespace Opm { + +template +EclMixingRateControls:: +EclMixingRateControls(const Schedule& schedule) + : schedule_(schedule) +{} + +template +EclMixingRateControls:: +EclMixingRateControls(const EclMixingRateControls& rhs) + : schedule_(rhs.schedule_) +{ + *this = rhs; +} + +template +EclMixingRateControls +EclMixingRateControls:: +serializationTestObject(const Schedule& schedule) +{ + EclMixingRateControls result(schedule); + result.lastRv_ = {21.0}; + result.maxDRv_ = {22.0, 23.0}; + result.convectiveDrs_ = {24.0, 25.0, 26.0}; + result.lastRs_ = {27.0}; + result.maxDRs_ = {28.0}; + result.dRsDtOnlyFreeGas_ = {false, true}; + + return result; +} + +template +bool EclMixingRateControls:: +operator==(const EclMixingRateControls& rhs) const +{ + return this->lastRv_ == rhs.lastRv_ && + this->maxDRv_ == rhs.maxDRv_ && + this->convectiveDrs_ == rhs.convectiveDrs_ && + this->lastRs_ == rhs.lastRs_ && + this->maxDRs_ == rhs.maxDRs_ && + this->dRsDtOnlyFreeGas_ == rhs.dRsDtOnlyFreeGas_; +} + +template +EclMixingRateControls& +EclMixingRateControls:: +operator=(const EclMixingRateControls& rhs) +{ + this->lastRv_ = rhs.lastRv_; + this->maxDRv_ = rhs.maxDRv_; + this->convectiveDrs_ = rhs.convectiveDrs_; + this->lastRs_ = rhs.lastRs_; + this->maxDRs_ = rhs.maxDRs_; + this->dRsDtOnlyFreeGas_ = rhs.dRsDtOnlyFreeGas_; + + return *this; +} + +template +void EclMixingRateControls:: +init(std::size_t numDof, int episodeIdx, const unsigned ntpvt) +{ + // deal with DRSDT + //TODO We may want to only allocate these properties only if active. + //But since they may be activated at later time we need some more + //intrastructure to handle it + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + maxDRv_.resize(ntpvt, 1e30); + lastRv_.resize(numDof, 0.0); + maxDRs_.resize(ntpvt, 1e30); + dRsDtOnlyFreeGas_.resize(ntpvt, false); + lastRs_.resize(numDof, 0.0); + maxDRv_.resize(ntpvt, 1e30); + lastRv_.resize(numDof, 0.0); + if (this->drsdtConvective(episodeIdx)) { + convectiveDrs_.resize(numDof, 1.0); + } + } +} + +template +bool EclMixingRateControls:: +drsdtActive(int episodeIdx) const +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + return (oilVaporizationControl.drsdtActive() && bothOilGasActive); +} + +template +bool EclMixingRateControls:: +drvdtActive(int episodeIdx) const +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + return (oilVaporizationControl.drvdtActive() && bothOilGasActive); +} + +template +bool EclMixingRateControls:: +drsdtConvective(int episodeIdx) const +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + return (oilVaporizationControl.drsdtConvective() && bothOilGasActive); +} + +template +void EclMixingRateControls:: +updateExplicitQuantities(const int episodeIdx, + const Scalar timeStepSize) +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + if (this->drsdtActive(episodeIdx)) { + // DRSDT is enabled + for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) + maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * timeStepSize; + } + + if (this->drvdtActive(episodeIdx)) { + // DRVDT is enabled + for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) + maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx) * timeStepSize; + } +} + +template +void EclMixingRateControls:: +updateLastValues(const unsigned elemIdx, + const Scalar Rs, + const Scalar Rv) +{ + if (!lastRs_.empty()) { + lastRs_[elemIdx] = Rs; + } + + if (!lastRv_.empty()) { + lastRv_[elemIdx] = Rv; + } +} + +template +void EclMixingRateControls:: +updateMaxValues(const int episodeIdx, + const Scalar timeStepSize) +{ + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + if (this->drsdtActive(episodeIdx)) { + // DRSDT is enabled + for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx) { + maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * timeStepSize; + } + } + + if (this->drvdtActive(episodeIdx)) { + // DRVDT is enabled + for (std::size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx) { + maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx) * timeStepSize; + } + } +} + +template +Scalar EclMixingRateControls:: +drsdtcon(const unsigned elemIdx, + int episodeIdx, + const int pvtRegionIdx) const +{ + if (convectiveDrs_.empty()) { + return 0; + } + + // The episode index is set to -1 in the initialization phase. + // Output drsdt value for index 0 + episodeIdx = std::max(episodeIdx, 0); + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + return oilVaporizationControl.getMaxDRSDT(pvtRegionIdx) * convectiveDrs_[elemIdx]; +} + +template +Scalar EclMixingRateControls:: +maxGasDissolutionFactor(const unsigned timeIdx, + const unsigned globalDofIdx, + const int episodeIdx, + const int pvtRegionIdx) const +{ + if (!this->drsdtActive(episodeIdx) || maxDRs_[pvtRegionIdx] < 0.0) { + return std::numeric_limits::max() / 2.0; + } + + Scalar scaling = 1.0; + if (this->drsdtConvective(episodeIdx)) { + scaling = convectiveDrs_[globalDofIdx]; + } + + // this is a bit hacky because it assumes that a time discretization with only + // two time indices is used. + if (timeIdx == 0) { + return lastRs_[globalDofIdx] + maxDRs_[pvtRegionIdx] * scaling; + } else { + return lastRs_[globalDofIdx]; + } +} + +template +Scalar EclMixingRateControls:: +maxOilVaporizationFactor(const unsigned timeIdx, + const unsigned globalDofIdx, + const int episodeIdx, + const int pvtRegionIdx) const +{ + if (!this->drvdtActive(episodeIdx) || maxDRv_[pvtRegionIdx] < 0.0) { + return std::numeric_limits::max() / 2.0; + } + + // this is a bit hacky because it assumes that a time discretization with only + // two time indices is used. + if (timeIdx == 0) { + return lastRv_[globalDofIdx] + maxDRv_[pvtRegionIdx]; + } else { + return lastRv_[globalDofIdx]; + } +} +template +void EclMixingRateControls:: +updateConvectiveDRsDt_(const unsigned compressedDofIdx, + const Scalar t, + const Scalar p, + const Scalar rs, + const Scalar so, + const Scalar poro, + const Scalar permz, + const Scalar distZ, + const Scalar gravity, + const int pvtRegionIndex) +{ + const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p); + const Scalar saturatedInvB + = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p); + const Scalar rsZero = 0.0; + const Scalar pureDensity + = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero) + * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex); + const Scalar saturatedDensity = saturatedInvB + * (FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex) + + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex)); + const Scalar deltaDensity = saturatedDensity - pureDensity; + const Scalar visc = FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs); + // Note that for so = 0 this gives no limits (inf) for the dissolution rate + // Also we restrict the effect of convective mixing to positive density differences + // i.e. we only allow for fingers moving downward + convectiveDrs_[compressedDofIdx] + = permz * rssat * max(0.0, deltaDensity) * gravity / (so * visc * distZ * poro); +} + +template class EclMixingRateControls, double>; + +} // namespace Opm diff --git a/ebos/eclmixingratecontrols.hh b/ebos/eclmixingratecontrols.hh new file mode 100644 index 0000000000..89d3d8eb2a --- /dev/null +++ b/ebos/eclmixingratecontrols.hh @@ -0,0 +1,183 @@ +// -*- 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::EclProblem + */ +#ifndef ECL_MIXING_RATE_CONTROLS_HH +#define ECL_MIXING_RATE_CONTROLS_HH + +#include + +#include +#include + +#include +#include + +namespace Opm { + +class EclipseState; + +//! \brief Class handling mixing rate controls for an EclProblem. +template +class EclMixingRateControls { +public: + EclMixingRateControls(const Schedule& schedule); + EclMixingRateControls(const EclMixingRateControls& rhs); + + static EclMixingRateControls serializationTestObject(const Schedule& schedule); + + bool operator==(const EclMixingRateControls& rhs) const; + EclMixingRateControls& operator=(const EclMixingRateControls& rhs); + + void init(std::size_t numDof, int episodeIdx, const unsigned ntpvt); + + bool drsdtActive(int episodeIdx) const; + bool drvdtActive(int episodeIdx) const; + bool drsdtConvective(int episodeIdx) const; + + /*! + * \brief Returns the dynamic drsdt convective mixing value + */ + Scalar drsdtcon(const unsigned elemIdx, + int episodeIdx, + const int pvtRegionIdx) const; + + /*! + * \brief Returns the maximum value of the gas dissolution factor at the current time + * for a given degree of freedom. + */ + Scalar maxGasDissolutionFactor(unsigned timeIdx, + unsigned globalDofIdx, + const int episodeIdx, + const int pvtRegionIdx) const; + + /*! + * \brief Returns the maximum value of the oil vaporization factor at the current + * time for a given degree of freedom. + */ + Scalar maxOilVaporizationFactor(const unsigned timeIdx, + const unsigned globalDofIdx, + const int episodeIdx, + const int pvtRegionIdx) const; + + void updateExplicitQuantities(const int episodeIdx, + const Scalar timeStepSize); + + void updateLastValues(const unsigned elemIdx, + const Scalar Rs, + const Scalar Rv); + + void updateMaxValues(const int episodeIdx, + const Scalar timeStepSize); + + template + void serializeOp(Serializer& serializer) + { + serializer(lastRv_); + serializer(maxDRv_); + serializer(convectiveDrs_); + serializer(lastRs_); + serializer(maxDRs_); + serializer(dRsDtOnlyFreeGas_); + } + + template + void update(unsigned compressedDofIdx, + const IntensiveQuantities& iq, + const int episodeIdx, + const Scalar gravity, + const Scalar permZ, + const Scalar distZ, + const int pvtRegionIdx, + const std::array& active) + { + if (active[0]) { + // This implements the convective DRSDT as described in + // Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow + // simulator" Submitted to TCCS 11, 2021 + const auto& fs = iq.fluidState(); + this->updateConvectiveDRsDt_(compressedDofIdx, + getValue(fs.temperature(FluidSystem::oilPhaseIdx)), + getValue(fs.pressure(FluidSystem::oilPhaseIdx)), + getValue(fs.Rs()), + getValue(fs.saturation(FluidSystem::oilPhaseIdx)), + getValue(iq.porosity()), + permZ, + distZ, + gravity, + fs.pvtRegionIndex()); + } + + if (active[1]) { + const auto& fs = iq.fluidState(); + + using FluidState = typename std::decay::type; + + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); + constexpr Scalar freeGasMinSaturation_ = 1e-7; + if (oilVaporizationControl.getOption(pvtRegionIdx) || + fs.saturation(FluidSystem::gasPhaseIdx) > freeGasMinSaturation_) { + lastRs_[compressedDofIdx] + = BlackOil::template getRs_(fs, iq.pvtRegionIndex()); + } + else + lastRs_[compressedDofIdx] = std::numeric_limits::infinity(); + } + + if (active[2]) { + const auto& fs = iq.fluidState(); + using FluidState = typename std::decay::type; + lastRv_[compressedDofIdx] + = BlackOil::template getRv_(fs, iq.pvtRegionIndex()); + } + } + +private: + void updateConvectiveDRsDt_(const unsigned compressedDofIdx, + const Scalar t, + const Scalar p, + const Scalar rs, + const Scalar so, + const Scalar poro, + const Scalar permz, + const Scalar distZ, + const Scalar gravity, + const int pvtRegionIndex); + + std::vector lastRv_; + std::vector maxDRv_; + + std::vector convectiveDrs_; + std::vector lastRs_; + std::vector maxDRs_; + std::vector dRsDtOnlyFreeGas_; // apply the DRSDT rate limit only to cells that exhibit free gas + + const Schedule& schedule_; +}; + +} // namespace Opm + +#endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 0170666c87..b9c72ba3b6 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -375,7 +375,14 @@ public: this->initFluidSystem_(); // deal with DRSDT - this->initDRSDT_(this->model().numGridDof(), this->episodeIndex()); + this->mixControls_.init(this->model().numGridDof(), + this->episodeIndex(), + eclState.runspec().tabdims().getNumPVTTables()); + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0); + } this->readRockParameters_(simulator.vanguard().cellCenterDepths(), [&simulator](const unsigned idx) @@ -1166,22 +1173,9 @@ public: */ Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const { - int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx); - int episodeIdx = this->episodeIndex(); - if (!this->drsdtActive_(episodeIdx) || this->maxDRs_[pvtRegionIdx] < 0.0) - return std::numeric_limits::max()/2.0; - - Scalar scaling = 1.0; - if (this->drsdtConvective_(episodeIdx)) { - scaling = this->convectiveDrs_[globalDofIdx]; - } - - // this is a bit hacky because it assumes that a time discretization with only - // two time indices is used. - if (timeIdx == 0) - return this->lastRs_[globalDofIdx] + this->maxDRs_[pvtRegionIdx] * scaling; - else - return this->lastRs_[globalDofIdx]; + return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx, + this->episodeIndex(), + this->pvtRegionIndex(globalDofIdx)); } /*! @@ -1190,17 +1184,9 @@ public: */ Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const { - int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx); - int episodeIdx = this->episodeIndex(); - if (!this->drvdtActive_(episodeIdx) || this->maxDRv_[pvtRegionIdx] < 0.0) - return std::numeric_limits::max()/2.0; - - // this is a bit hacky because it assumes that a time discretization with only - // two time indices is used. - if (timeIdx == 0) - return this->lastRv_[globalDofIdx] + this->maxDRv_[pvtRegionIdx]; - else - return this->lastRv_[globalDofIdx]; + return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx, + this->episodeIndex(), + this->pvtRegionIndex(globalDofIdx)); } /*! @@ -1214,8 +1200,8 @@ public: bool recycleFirstIterationStorage() const { int episodeIdx = this->episodeIndex(); - return !this->drsdtActive_(episodeIdx) && - !this->drvdtActive_(episodeIdx) && + return !this->mixControls_.drsdtActive(episodeIdx) && + !this->mixControls_.drvdtActive(episodeIdx) && this->rockCompPoroMultWc_.empty() && this->rockCompPoroMult_.empty(); } @@ -1712,82 +1698,32 @@ protected: // update the "last Rs" values for all elements, including the ones in the ghost // and overlap regions int episodeIdx = this->episodeIndex(); - std::array active{this->drsdtConvective_(episodeIdx), - this->drsdtActive_(episodeIdx), - this->drvdtActive_(episodeIdx)}; - if (!active[0] && !active[1] && !active[2]) - return; + std::array active{this->mixControls_.drsdtConvective(episodeIdx), + this->mixControls_.drsdtActive(episodeIdx), + this->mixControls_.drvdtActive(episodeIdx)}; + if (!active[0] && !active[1] && !active[2]) { + return; + } this->updateProperty_("EclProblem::updateCompositionChangeLimits_()) failed:", - [this,episodeIdx,active](unsigned compressedDofIdx, const IntensiveQuantities& iq) + [this,episodeIdx,active](unsigned compressedDofIdx, + const IntensiveQuantities& iq) { - this->updateCompositionChangeLimits_(compressedDofIdx, - iq, - episodeIdx, - active); + const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx); + const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0; + const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx); + this->mixControls_.update(compressedDofIdx, + iq, + episodeIdx, + this->gravity_[dim - 1], + perm[dim - 1][dim - 1], + distZ, + pvtRegionIdx, + active); } ); } - void updateCompositionChangeLimits_(unsigned compressedDofIdx, const IntensiveQuantities& iq,int episodeIdx, const std::array& active) - { - auto& simulator = this->simulator(); - auto& vanguard = simulator.vanguard(); - if (active[0]) { - // This implements the convective DRSDT as described in - // Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow - // simulator" Submitted to TCCS 11, 2021 - const Scalar g = this->gravity_[dim - 1]; - const DimMatrix& perm = intrinsicPermeability(compressedDofIdx); - const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability - const Scalar distZ = vanguard.cellThickness(compressedDofIdx); - const auto& fs = iq.fluidState(); - const Scalar t = getValue(fs.temperature(FluidSystem::oilPhaseIdx)); - const Scalar p = getValue(fs.pressure(FluidSystem::oilPhaseIdx)); - const Scalar so = getValue(fs.saturation(FluidSystem::oilPhaseIdx)); - const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), t, p); - const Scalar saturatedInvB - = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), t, p); - const Scalar rsZero = 0.0; - const Scalar pureDensity - = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), t, p, rsZero) - * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex()); - const Scalar saturatedDensity = saturatedInvB - * (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex()) - + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex())); - const Scalar deltaDensity = saturatedDensity - pureDensity; - const Scalar rs = getValue(fs.Rs()); - const Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(), t, p, rs); - const Scalar poro = getValue(iq.porosity()); - // Note that for so = 0 this gives no limits (inf) for the dissolution rate - // Also we restrict the effect of convective mixing to positive density differences - // i.e. we only allow for fingers moving downward - this->convectiveDrs_[compressedDofIdx] - = permz * rssat * max(0.0, deltaDensity) * g / (so * visc * distZ * poro); - } - - if (active[1]) { - const auto& fs = iq.fluidState(); - - using FluidState = typename std::decay::type; - - int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx); - const auto& oilVaporizationControl = vanguard.schedule()[episodeIdx].oilvap(); - if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_) - this->lastRs_[compressedDofIdx] - = BlackOil::template getRs_(fs, iq.pvtRegionIndex()); - else - this->lastRs_[compressedDofIdx] = std::numeric_limits::infinity(); - } - - if (active[2]) { - const auto& fs = iq.fluidState(); - using FluidState = typename std::decay::type; - this->lastRv_[compressedDofIdx] - = BlackOil::template getRv_(fs, iq.pvtRegionIndex()); - } - } - bool updateMaxOilSaturation_() { OPM_TIMEBLOCK(updateMaxOilSaturation); @@ -2065,13 +2001,7 @@ protected: this->solventSaturation_[elemIdx] = ssol; } - if (! this->lastRs_.empty()) { - this->lastRs_[elemIdx] = elemFluidState.Rs(); - } - - if (! this->lastRv_.empty()) { - this->lastRv_[elemIdx] = elemFluidState.Rv(); - } + this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv()); if constexpr (enablePolymer) this->polymer_.concentration[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx); @@ -2086,16 +2016,7 @@ protected: } const int episodeIdx = this->episodeIndex(); - const auto& oilVaporizationControl = simulator.vanguard().schedule()[episodeIdx].oilvap(); - if (this->drsdtActive_(episodeIdx)) - // DRSDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < this->maxDRs_.size(); ++pvtRegionIdx) - this->maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*simulator.timeStepSize(); - - if (this->drvdtActive_(episodeIdx)) - // DRVDT is enabled - for (size_t pvtRegionIdx = 0; pvtRegionIdx < this->maxDRv_.size(); ++pvtRegionIdx) - this->maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*simulator.timeStepSize(); + this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize()); // assign the restart solution to the current solution. note that we still need // to compute real initial solution after this because the initial fluid states @@ -2554,8 +2475,6 @@ private: std::vector initialFluidStates_; - constexpr static Scalar freeGasMinSaturation_ = 1e-7; - bool enableDriftCompensation_; GlobalEqVector drift_;