mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
338 lines
12 KiB
C++
338 lines
12 KiB
C++
// -*- 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 <http://www.gnu.org/licenses/>.
|
|
|
|
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::FlowProblem
|
|
*/
|
|
|
|
#include <config.h>
|
|
#include <opm/simulators/flow/MixingRateControls.hpp>
|
|
|
|
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
|
|
|
|
#include <algorithm>
|
|
#include <limits>
|
|
|
|
namespace Opm {
|
|
|
|
template<class FluidSystem>
|
|
MixingRateControls<FluidSystem>::
|
|
MixingRateControls(const Schedule& schedule)
|
|
: schedule_(schedule)
|
|
{}
|
|
|
|
template<class FluidSystem>
|
|
MixingRateControls<FluidSystem>::
|
|
MixingRateControls(const MixingRateControls& rhs)
|
|
: schedule_(rhs.schedule_)
|
|
{
|
|
*this = rhs;
|
|
}
|
|
|
|
template<class FluidSystem>
|
|
MixingRateControls<FluidSystem>
|
|
MixingRateControls<FluidSystem>::
|
|
serializationTestObject(const Schedule& schedule)
|
|
{
|
|
MixingRateControls<FluidSystem> 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<class FluidSystem>
|
|
bool MixingRateControls<FluidSystem>::
|
|
operator==(const MixingRateControls& 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<class FluidSystem>
|
|
MixingRateControls<FluidSystem>&
|
|
MixingRateControls<FluidSystem>::
|
|
operator=(const MixingRateControls& 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<class FluidSystem>
|
|
void MixingRateControls<FluidSystem>::
|
|
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
|
|
|
|
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);
|
|
if (this->drsdtConvective(episodeIdx)) {
|
|
convectiveDrs_.resize(numDof, 1.0);
|
|
}
|
|
}
|
|
|
|
template<class FluidSystem>
|
|
bool MixingRateControls<FluidSystem>::
|
|
drsdtActive(int episodeIdx) const
|
|
{
|
|
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
|
|
return (oilVaporizationControl.drsdtActive());
|
|
}
|
|
|
|
template<class FluidSystem>
|
|
bool MixingRateControls<FluidSystem>::
|
|
drvdtActive(int episodeIdx) const
|
|
{
|
|
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
|
|
return (oilVaporizationControl.drvdtActive());
|
|
}
|
|
|
|
template<class FluidSystem>
|
|
bool MixingRateControls<FluidSystem>::
|
|
drsdtConvective(int episodeIdx) const
|
|
{
|
|
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
|
|
return (oilVaporizationControl.drsdtConvective());
|
|
}
|
|
|
|
template<class FluidSystem>
|
|
void MixingRateControls<FluidSystem>::
|
|
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<class FluidSystem>
|
|
void MixingRateControls<FluidSystem>::
|
|
updateLastValues(const unsigned elemIdx,
|
|
const Scalar Rs,
|
|
const Scalar Rv)
|
|
{
|
|
if (!lastRs_.empty()) {
|
|
lastRs_[elemIdx] = Rs;
|
|
}
|
|
|
|
if (!lastRv_.empty()) {
|
|
lastRv_[elemIdx] = Rv;
|
|
}
|
|
}
|
|
|
|
template<class FluidSystem>
|
|
void MixingRateControls<FluidSystem>::
|
|
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<class FluidSystem>
|
|
typename MixingRateControls<FluidSystem>::Scalar
|
|
MixingRateControls<FluidSystem>::
|
|
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<class FluidSystem>
|
|
typename MixingRateControls<FluidSystem>::Scalar
|
|
MixingRateControls<FluidSystem>::
|
|
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<Scalar>::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<class FluidSystem>
|
|
typename MixingRateControls<FluidSystem>::Scalar
|
|
MixingRateControls<FluidSystem>::
|
|
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<Scalar>::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<class FluidSystem>
|
|
void MixingRateControls<FluidSystem>::
|
|
updateConvectiveDRsDt_(const unsigned compressedDofIdx,
|
|
const Scalar t,
|
|
const Scalar p,
|
|
const Scalar pg,
|
|
const Scalar rs,
|
|
const Scalar sg,
|
|
const Scalar poro,
|
|
const Scalar permz,
|
|
const Scalar distZ,
|
|
const Scalar gravity,
|
|
const Scalar salt,
|
|
const Scalar Xhi,
|
|
const Scalar Psi,
|
|
const Scalar omegainn,
|
|
const int pvtRegionIndex)
|
|
{
|
|
const Scalar rssat = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
|
|
FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p, salt) :
|
|
FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p);
|
|
const Scalar saturatedInvB = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
|
|
FluidSystem::waterPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p, salt) :
|
|
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p);
|
|
const Scalar rsZero = 0.0;
|
|
const Scalar sg_max = 1.0;
|
|
const Scalar pureDensity = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
|
|
(FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero, salt)
|
|
* FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIndex)) :
|
|
(FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero)
|
|
* FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex));
|
|
const Scalar saturatedDensity = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
|
|
(saturatedInvB * (FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIndex)
|
|
+ rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex))) :
|
|
(saturatedInvB * (FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex)
|
|
+ rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex)));
|
|
Scalar deltaDensity = saturatedDensity - pureDensity;
|
|
const Scalar visc = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
|
|
FluidSystem::waterPvt().viscosity(pvtRegionIndex, t, p, rs, salt) :
|
|
FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs);
|
|
|
|
// Note that for sLiquid = 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
|
|
|
|
Scalar co2Density =
|
|
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIndex,t,p,Scalar{0.0} /*=Rv*/, Scalar{0.0} /*=Rvw*/) *
|
|
FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex);
|
|
Scalar factor = 1.0;
|
|
Scalar X = (rs - rssat * sg) / (rssat * ( 1.0 - sg));
|
|
Scalar omega = 0.0;
|
|
const Scalar pCap = Opm::abs(pg - p);
|
|
if ((rs >= (rssat * sg)) || (pCap < 1e-12)) {
|
|
if (X > Psi) {
|
|
factor = 0.0;
|
|
omega = omegainn;
|
|
}
|
|
} else {
|
|
factor /= Xhi;
|
|
deltaDensity = (saturatedDensity - co2Density);
|
|
}
|
|
|
|
convectiveDrs_[compressedDofIdx] =
|
|
factor * permz * rssat * max(Scalar{0.0}, deltaDensity) *
|
|
gravity / ( std::max(sg_max - sg, Scalar{0.0}) *
|
|
visc * distZ * poro) + (omega/Xhi);
|
|
}
|
|
|
|
#define INSTANTIATE_TYPE(T) \
|
|
template class MixingRateControls<BlackOilFluidSystem<T,BlackOilDefaultIndexTraits>>;
|
|
|
|
INSTANTIATE_TYPE(double)
|
|
|
|
#if FLOW_INSTANTIATE_FLOAT
|
|
INSTANTIATE_TYPE(float)
|
|
#endif
|
|
|
|
} // namespace Opm
|