OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #628 from osae/extboModels

Browse Source 
Alternative solvent extension for the black oil model.
This commit is contained in:
Tor Harald Sandve 2020-11-18 14:01:11 +01:00 committed by GitHub
commit 5d9783f9b1
10 changed files with 1156 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

992
opm/models/blackoil/blackoilextbomodules.hh Normal file
View File

@ -0,0 +1,992 @@
// -*- 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
*
* \brief Contains the classes required to extend the black-oil model by solvent component.
* For details, refer:
* [*] T.H. Sandve, O. Sævareid and I. Aavatsmark: “Improved Extended Blackoil Formulation
* for CO2 EOR Simulations.” in ECMOR XVII – The 17th European Conference on the
* Mathematics of Oil Recovery, September 2020.
*/
#ifndef EWOMS_BLACK_OIL_EXTBO_MODULE_HH
#define EWOMS_BLACK_OIL_EXTBO_MODULE_HH
#include "blackoilproperties.hh"
//#include <opm/models/io/vtkBlackOilExtboModule.hh> //TODO: Missing ...
#include <opm/models/common/quantitycallbacks.hh>
#include <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/common/UniformXTabulated2DFunction.hpp>
#if HAVE_ECL_INPUT
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SsfnTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/MsfnTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PmiscTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/MiscTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SorwmisTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SgcwmisTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
#endif
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/Unused.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <dune/common/fvector.hh>
#include <string>
namespace Opm {
/*!
* \ingroup BlackOil
* \brief Contains the high level supplements required to extend the black oil
* model.
*/
template <class TypeTag, bool enableExtboV = getPropValue<TypeTag, Properties::EnableExtbo>()>
class BlackOilExtboModule
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Model = GetPropType<TypeTag, Properties::Model>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using EqVector = GetPropType<TypeTag, Properties::EqVector>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef typename Opm::Tabulated1DFunction<Scalar> TabulatedFunction;
typedef typename Opm::UniformXTabulated2DFunction<Scalar> Tabulated2DFunction;
static constexpr unsigned zFractionIdx = Indices::zFractionIdx;
static constexpr unsigned contiZfracEqIdx = Indices::contiZfracEqIdx;
static constexpr unsigned enableExtbo = enableExtboV;
static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
static constexpr unsigned numPhases = FluidSystem::numPhases;
static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
public:
#if HAVE_ECL_INPUT
/*!
* \brief Initialize all internal data structures needed by the solvent module
*/
static void initFromState(const Opm::EclipseState& eclState)
{
// some sanity checks: if extended BO is enabled, the PVTSOL keyword must be
// present, if extended BO is disabled the keyword must not be present.
if (enableExtbo && !eclState.runspec().phases().active(Phase::ZFRACTION))
throw std::runtime_error("Extended black oil treatment requested at compile "
"time, but the deck does not contain the PVTSOL keyword");
else if (!enableExtbo && eclState.runspec().phases().active(Phase::ZFRACTION))
throw std::runtime_error("Extended black oil treatment disabled at compile time, but the deck "
"contains the PVTSOL keyword");
if (!eclState.runspec().phases().active(Phase::ZFRACTION))
return; // solvent treatment is supposed to be disabled
// pvt properties from kw PVTSOL:
const auto& tableManager = eclState.getTableManager();
const auto& pvtsolTables = tableManager.getPvtsolTables();
size_t numPvtRegions = pvtsolTables.size();
BO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
BG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
RS_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
RV_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
X_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
Y_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
VISCO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
VISCG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
PBUB_RS_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
PBUB_RV_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
zLim_.resize(numPvtRegions);
const bool extractCmpFromPvt = true; //<false>: Default values used in [*]
oilCmp_.resize(numPvtRegions);
gasCmp_.resize(numPvtRegions);
for (unsigned regionIdx = 0; regionIdx < numPvtRegions; ++ regionIdx) {
const auto& pvtsolTable = pvtsolTables[regionIdx];
const auto& saturatedTable = pvtsolTable.getSaturatedTable();
assert(saturatedTable.numRows() > 1);
std::vector<Scalar> oilCmp(saturatedTable.numRows(), -4.0e-9); //Default values used in [*]
std::vector<Scalar> gasCmp(saturatedTable.numRows(), -0.08); //-------------"-------------
zLim_[regionIdx] = 0.7; //-------------"-------------
std::vector<Scalar> zArg(saturatedTable.numRows(), 0.0);
for (unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
Scalar ZCO2 = saturatedTable.get("ZCO2", outerIdx);
zArg[outerIdx] = ZCO2;
BO_[regionIdx].appendXPos(ZCO2);
BG_[regionIdx].appendXPos(ZCO2);
RS_[regionIdx].appendXPos(ZCO2);
RV_[regionIdx].appendXPos(ZCO2);
X_[regionIdx].appendXPos(ZCO2);
Y_[regionIdx].appendXPos(ZCO2);
VISCO_[regionIdx].appendXPos(ZCO2);
VISCG_[regionIdx].appendXPos(ZCO2);
PBUB_RS_[regionIdx].appendXPos(ZCO2);
PBUB_RV_[regionIdx].appendXPos(ZCO2);
const auto& underSaturatedTable = pvtsolTable.getUnderSaturatedTable(outerIdx);
size_t numRows = underSaturatedTable.numRows();
Scalar bo0=0.0;
Scalar po0=0.0;
for (unsigned innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
Scalar po = underSaturatedTable.get("P", innerIdx);
Scalar bo = underSaturatedTable.get("B_O", innerIdx);
Scalar bg = underSaturatedTable.get("B_G", innerIdx);
Scalar rs = underSaturatedTable.get("RS", innerIdx)+innerIdx*1.0e-10;
Scalar rv = underSaturatedTable.get("RV", innerIdx)+innerIdx*1.0e-10;
Scalar xv = underSaturatedTable.get("XVOL", innerIdx);
Scalar yv = underSaturatedTable.get("YVOL", innerIdx);
Scalar mo = underSaturatedTable.get("MU_O", innerIdx);
Scalar mg = underSaturatedTable.get("MU_G", innerIdx);
if (bo0 > bo) { // This is undersaturated oil-phase for ZCO2 <= zLim ...
// Here we assume tabulated bo to decay beyond boiling point
if (extractCmpFromPvt) {
Scalar cmpFactor = (bo-bo0)/(po-po0);
oilCmp[outerIdx] = cmpFactor;
zLim_[regionIdx] = ZCO2;
//std::cout << "### cmpFactorOil: " << cmpFactor << " zLim: " << zLim_[regionIdx] << std::endl;
}
break;
} else if (bo0 == bo) { // This is undersaturated gas-phase for ZCO2 > zLim ...
// Here we assume tabulated bo to be constant extrapolated beyond dew point
if (innerIdx+1 < numRows && ZCO2<1.0 && extractCmpFromPvt) {
Scalar rvNxt = underSaturatedTable.get("RV", innerIdx+1)+innerIdx*1.0e-10;
Scalar bgNxt = underSaturatedTable.get("B_G", innerIdx+1);
Scalar cmpFactor = (bgNxt-bg)/(rvNxt-rv);
gasCmp[outerIdx] = cmpFactor;
//std::cout << "### cmpFactorGas: " << cmpFactor << " zLim: " << zLim_[regionIdx] << std::endl;
}
BO_[regionIdx].appendSamplePoint(outerIdx,po,bo);
BG_[regionIdx].appendSamplePoint(outerIdx,po,bg);
RS_[regionIdx].appendSamplePoint(outerIdx,po,rs);
RV_[regionIdx].appendSamplePoint(outerIdx,po,rv);
X_[regionIdx].appendSamplePoint(outerIdx,po,xv);
Y_[regionIdx].appendSamplePoint(outerIdx,po,yv);
VISCO_[regionIdx].appendSamplePoint(outerIdx,po,mo);
VISCG_[regionIdx].appendSamplePoint(outerIdx,po,mg);
break;
}
bo0=bo;
po0=po;
BO_[regionIdx].appendSamplePoint(outerIdx,po,bo);
BG_[regionIdx].appendSamplePoint(outerIdx,po,bg);
RS_[regionIdx].appendSamplePoint(outerIdx,po,rs);
RV_[regionIdx].appendSamplePoint(outerIdx,po,rv);
X_[regionIdx].appendSamplePoint(outerIdx,po,xv);
Y_[regionIdx].appendSamplePoint(outerIdx,po,yv);
VISCO_[regionIdx].appendSamplePoint(outerIdx,po,mo);
VISCG_[regionIdx].appendSamplePoint(outerIdx,po,mg);
// rs,rv -> pressure
PBUB_RS_[regionIdx].appendSamplePoint(outerIdx, rs, po);
PBUB_RV_[regionIdx].appendSamplePoint(outerIdx, rv, po);
}
}
oilCmp_[regionIdx].setXYContainers(zArg, oilCmp, /*sortInput=*/false);
gasCmp_[regionIdx].setXYContainers(zArg, gasCmp, /*sortInput=*/false);
}
// Reference density for pure z-component taken from kw SDENSITY
const auto& sdensityTables = eclState.getTableManager().getSolventDensityTables();
if (sdensityTables.size() == numPvtRegions) {
zReferenceDensity_.resize(numPvtRegions);
for (unsigned regionIdx = 0; regionIdx < numPvtRegions; ++ regionIdx) {
Scalar rhoRefS = sdensityTables[regionIdx].getSolventDensityColumn().front();
zReferenceDensity_[regionIdx]=rhoRefS;
}
}
else
throw std::runtime_error("Extbo: kw SDENSITY is missing or not aligned with NTPVT\n");
}
#endif
/*!
* \brief Register all run-time parameters for the black-oil solvent module.
*/
static void registerParameters()
{
if (!enableExtbo)
// extBO have disabled at compile time
return;
//Opm::VtkBlackOilExtboModule<TypeTag>::registerParameters();
}
/*!
* \brief Register all solvent specific VTK and ECL output modules.
*/
static void registerOutputModules(Model& model,
Simulator& simulator)
{
if (!enableExtbo)
// extBO have disabled at compile time
return;
//model.addOutputModule(new Opm::VtkBlackOilExtboModule<TypeTag>(simulator));
}
static bool primaryVarApplies(unsigned pvIdx)
{
if (!enableExtbo)
// extBO have disabled at compile time
return false;
return pvIdx == zFractionIdx;
}
static std::string primaryVarName(unsigned pvIdx OPM_OPTIM_UNUSED)
{
assert(primaryVarApplies(pvIdx));
return "z_fraction";
}
static Scalar primaryVarWeight(unsigned pvIdx OPM_OPTIM_UNUSED)
{
assert(primaryVarApplies(pvIdx));
// TODO: it may be beneficial to chose this differently.
return static_cast<Scalar>(1.0);
}
static bool eqApplies(unsigned eqIdx)
{
if (!enableExtbo)
return false;
return eqIdx == contiZfracEqIdx;
}
static std::string eqName(unsigned eqIdx OPM_OPTIM_UNUSED)
{
assert(eqApplies(eqIdx));
return "conti^solvent";
}
static Scalar eqWeight(unsigned eqIdx OPM_OPTIM_UNUSED)
{
assert(eqApplies(eqIdx));
// TODO: it may be beneficial to chose this differently.
return static_cast<Scalar>(1.0);
}
template <class LhsEval>
static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
const IntensiveQuantities& intQuants)
{
if (!enableExtbo)
return;
if (blackoilConserveSurfaceVolume) {
storage[contiZfracEqIdx] =
Toolbox::template decay<LhsEval>(intQuants.porosity())
* Toolbox::template decay<LhsEval>(intQuants.yVolume())
* Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(gasPhaseIdx))
* Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(gasPhaseIdx));
if (FluidSystem::enableDissolvedGas()) { // account for dissolved z in oil phase
storage[contiZfracEqIdx] +=
Toolbox::template decay<LhsEval>(intQuants.porosity())
* Toolbox::template decay<LhsEval>(intQuants.xVolume())
* Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs())
* Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(oilPhaseIdx))
* Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(oilPhaseIdx));
}
// Reg. terms: Preliminary attempt to avoid singular behaviour when solvent is invading a pure water
// region. Results seems insensitive to the weighting factor.
// TODO: Further investigations ...
const Scalar regWghtFactor = 1.0e-6;
storage[contiZfracEqIdx] += regWghtFactor*(1.0-Toolbox::template decay<LhsEval>(intQuants.zFraction()))
+ regWghtFactor*Toolbox::template decay<LhsEval>(intQuants.porosity())
* Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(gasPhaseIdx))
* Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(gasPhaseIdx));
storage[contiZfracEqIdx-1] += regWghtFactor*Toolbox::template decay<LhsEval>(intQuants.zFraction());
}
else {
throw std::runtime_error("Only component conservation in terms of surface volumes is implemented. ");
}
}
static void computeFlux(RateVector& flux,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
{
if (!enableExtbo)
return;
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
if (blackoilConserveSurfaceVolume) {
unsigned inIdx = extQuants.interiorIndex();
unsigned upIdxGas = static_cast<unsigned>(extQuants.upstreamIndex(gasPhaseIdx));
const auto& upGas = elemCtx.intensiveQuantities(upIdxGas, timeIdx);
const auto& fsGas = upGas.fluidState();
if (upIdxGas == inIdx) {
flux[contiZfracEqIdx] =
extQuants.volumeFlux(gasPhaseIdx)
* (upGas.yVolume())
* fsGas.invB(gasPhaseIdx);
}
else {
flux[contiZfracEqIdx] =
extQuants.volumeFlux(gasPhaseIdx)
* (Opm::decay<Scalar>(upGas.yVolume()))
* Opm::decay<Scalar>(fsGas.invB(gasPhaseIdx));
}
if (FluidSystem::enableDissolvedGas()) { // account for dissolved z in oil phase
unsigned upIdxOil = static_cast<unsigned>(extQuants.upstreamIndex(oilPhaseIdx));
const auto& upOil = elemCtx.intensiveQuantities(upIdxOil, timeIdx);
const auto& fsOil = upOil.fluidState();
if (upIdxOil == inIdx) {
flux[contiZfracEqIdx] +=
extQuants.volumeFlux(oilPhaseIdx)
* upOil.xVolume()
* fsOil.Rs()
* fsOil.invB(oilPhaseIdx);
}
else {
flux[contiZfracEqIdx] +=
extQuants.volumeFlux(oilPhaseIdx)
* Opm::decay<Scalar>(upOil.xVolume())
* Opm::decay<Scalar>(fsOil.Rs())
* Opm::decay<Scalar>(fsOil.invB(oilPhaseIdx));
}
}
}
else {
throw std::runtime_error("Only component conservation in terms of surface volumes is implemented. ");
}
}
/*!
* \brief Assign the solvent specific primary variables to a PrimaryVariables object
*/
static void assignPrimaryVars(PrimaryVariables& priVars,
Scalar zFraction)
{
if (!enableExtbo)
return;
priVars[zFractionIdx] = zFraction;
}
/*!
* \brief Do a Newton-Raphson update the primary variables of the solvents.
*/
static void updatePrimaryVars(PrimaryVariables& newPv,
const PrimaryVariables& oldPv,
const EqVector& delta)
{
if (!enableExtbo)
return;
// do a plain unchopped Newton update
newPv[zFractionIdx] = oldPv[zFractionIdx] - delta[zFractionIdx];
}
/*!
* \brief Return how much a Newton-Raphson update is considered an error
*/
static Scalar computeUpdateError(const PrimaryVariables& oldPv OPM_UNUSED,
const EqVector& delta OPM_UNUSED)
{
// do not consider consider the cange of solvent primary variables for
// convergence
// TODO: maybe this should be changed
return static_cast<Scalar>(0.0);
}
/*!
* \brief Return how much a residual is considered an error
*/
static Scalar computeResidualError(const EqVector& resid)
{
// do not weight the residual of solvents when it comes to convergence
return std::abs(Toolbox::scalarValue(resid[contiZfracEqIdx]));
}
template <class DofEntity>
static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
{
if (!enableExtbo)
return;
unsigned dofIdx = model.dofMapper().index(dof);
const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
outstream << priVars[zFractionIdx];
}
template <class DofEntity>
static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
{
if (!enableExtbo)
return;
unsigned dofIdx = model.dofMapper().index(dof);
PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
instream >> priVars0[zFractionIdx];
// set the primary variables for the beginning of the current time step.
priVars1 = priVars0[zFractionIdx];
}
static const Scalar xVolume(unsigned pvtRegionIdx, const Scalar& pressure, const Scalar& z) {
const auto& xVolumeTable = X_[pvtRegionIdx];
return xVolumeTable.eval(z, pressure);
}
static const Evaluation xVolume(unsigned pvtRegionIdx, const Evaluation& pressure, const Evaluation& z) {
const auto& xVolumeTable = X_[pvtRegionIdx];
return xVolumeTable.eval(z, pressure);
}
static const Scalar yVolume(unsigned pvtRegionIdx, const Scalar& pressure, const Scalar& z) {
const auto& yVolumeTable = Y_[pvtRegionIdx];
return yVolumeTable.eval(z, pressure);
}
static const Evaluation yVolume(unsigned pvtRegionIdx, const Evaluation& pressure, const Evaluation& z) {
const auto& yVolumeTable = Y_[pvtRegionIdx];
return yVolumeTable.eval(z, pressure);
}
static const Scalar pbubRs(unsigned pvtRegionIdx, const Scalar& z, const Scalar& rs) {
const auto& pbubRsTable = PBUB_RS_[pvtRegionIdx];
return pbubRsTable.eval(z, rs);
}
static const Evaluation pbubRs(unsigned pvtRegionIdx, const Evaluation& z, const Evaluation& rs) {
const auto& pbubRsTable = PBUB_RS_[pvtRegionIdx];
return pbubRsTable.eval(z, rs);
}
static const Scalar pbubRv(unsigned pvtRegionIdx, const Scalar& z, const Scalar& rv) {
const auto& pbubRvTable = PBUB_RV_[pvtRegionIdx];
return pbubRvTable.eval(z, rv);
}
static const Evaluation pbubRv(unsigned pvtRegionIdx, const Evaluation& z, const Evaluation& rv) {
const auto& pbubRvTable = PBUB_RV_[pvtRegionIdx];
return pbubRvTable.eval(z, rv);
}
static const Scalar oilViscosity(unsigned pvtRegionIdx, const Scalar& pressure, const Scalar& z) {
const auto& oilViscosityTable = VISCO_[pvtRegionIdx];
return oilViscosityTable.eval(z, pressure);
}
static const Evaluation oilViscosity(unsigned pvtRegionIdx, const Evaluation& pressure, const Evaluation& z) {
const auto& oilViscosityTable = VISCO_[pvtRegionIdx];
return oilViscosityTable.eval(z, pressure);
}
static const Scalar gasViscosity(unsigned pvtRegionIdx, const Scalar& pressure, const Scalar& z) {
const auto& gasViscosityTable = VISCG_[pvtRegionIdx];
return gasViscosityTable.eval(z, pressure);
}
static const Evaluation gasViscosity(unsigned pvtRegionIdx, const Evaluation& pressure, const Evaluation& z) {
const auto& gasViscosityTable = VISCG_[pvtRegionIdx];
return gasViscosityTable.eval(z, pressure);
}
static const Scalar bo(unsigned pvtRegionIdx, const Scalar& pressure, const Scalar& z) {
const auto& boTable = BO_[pvtRegionIdx];
return boTable.eval(z, pressure);
}
static const Evaluation bo(unsigned pvtRegionIdx, const Evaluation& pressure, const Evaluation& z) {
const auto& boTable = BO_[pvtRegionIdx];
return boTable.eval(z, pressure);
}
static const Scalar bg(unsigned pvtRegionIdx, const Scalar& pressure, const Scalar& z) {
const auto& bgTable = BG_[pvtRegionIdx];
return bgTable.eval(z, pressure);
}
static const Evaluation bg(unsigned pvtRegionIdx, const Evaluation& pressure, const Evaluation& z) {
const auto& bgTable = BG_[pvtRegionIdx];
return bgTable.eval(z, pressure);
}
static const Scalar rs(unsigned pvtRegionIdx, const Scalar& pressure, const Scalar& z) {
const auto& rsTable = RS_[pvtRegionIdx];
return rsTable.eval(z, pressure);
}
static const Evaluation rs(unsigned pvtRegionIdx, const Evaluation& pressure, const Evaluation& z) {
const auto& rsTable = RS_[pvtRegionIdx];
return rsTable.eval(z, pressure);
}
static const Scalar rv(unsigned pvtRegionIdx, const Scalar& pressure, const Scalar& z) {
const auto& rvTable = RV_[pvtRegionIdx];
return rvTable.eval(z, pressure);
}
static const Evaluation rv(unsigned pvtRegionIdx, const Evaluation& pressure, const Evaluation& z) {
const auto& rvTable = RV_[pvtRegionIdx];
return rvTable.eval(z, pressure);
}
static const Scalar referenceDensity(unsigned regionIdx) {
return zReferenceDensity_[regionIdx];
}
static const Scalar zLim(unsigned regionIdx) {
return zLim_[regionIdx];
}
static const Scalar oilCmp(unsigned pvtRegionIdx, const Scalar& z) {
const auto& oilCmpTable = oilCmp_[pvtRegionIdx];
return oilCmpTable.eval(z);
}
static const Evaluation oilCmp(unsigned pvtRegionIdx, const Evaluation& z) {
const auto& oilCmpTable = oilCmp_[pvtRegionIdx];
return oilCmpTable.eval(z);
}
static const Scalar gasCmp(unsigned pvtRegionIdx, const Scalar& z) {
const auto& gasCmpTable = gasCmp_[pvtRegionIdx];
return gasCmpTable.eval(z);
}
static const Evaluation gasCmp(unsigned pvtRegionIdx, const Evaluation& z) {
const auto& gasCmpTable = gasCmp_[pvtRegionIdx];
return gasCmpTable.eval(z);
}
private:
static std::vector<Tabulated2DFunction> X_;
static std::vector<Tabulated2DFunction> Y_;
static std::vector<Tabulated2DFunction> PBUB_RS_;
static std::vector<Tabulated2DFunction> PBUB_RV_;
static std::vector<Tabulated2DFunction> VISCO_;
static std::vector<Tabulated2DFunction> VISCG_;
static std::vector<Tabulated2DFunction> BO_;
static std::vector<Tabulated2DFunction> BG_;
static std::vector<Tabulated2DFunction> RS_;
static std::vector<Tabulated2DFunction> RV_;
static std::vector<Scalar> zReferenceDensity_;
static std::vector<Scalar> zLim_;
static std::vector<TabulatedFunction> oilCmp_;
static std::vector<TabulatedFunction> gasCmp_;
};
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::X_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::Y_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::PBUB_RS_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::PBUB_RV_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::VISCO_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::VISCG_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::BO_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::BG_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::RS_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Tabulated2DFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::RV_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Scalar>
BlackOilExtboModule<TypeTag, enableExtboV>::zReferenceDensity_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::Scalar>
BlackOilExtboModule<TypeTag, enableExtboV>::zLim_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::TabulatedFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::oilCmp_;
template <class TypeTag, bool enableExtboV>
std::vector<typename BlackOilExtboModule<TypeTag, enableExtboV>::TabulatedFunction>
BlackOilExtboModule<TypeTag, enableExtboV>::gasCmp_;
/*!
* \ingroup BlackOil
* \class Opm::BlackOilExtboIntensiveQuantities
*
* \brief Provides the volumetric quantities required for the equations needed by the
* solvents extension of the black-oil model.
*/
template <class TypeTag, bool enableExtboV = getPropValue<TypeTag, Properties::EnableExtbo>()>
class BlackOilExtboIntensiveQuantities
{
using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using ExtboModule = BlackOilExtboModule<TypeTag>;
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
static constexpr int zFractionIdx = Indices::zFractionIdx;
static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
static constexpr double cutOff = 1e-12;
public:
/*!
* \brief Compute extended pvt properties from table lookups.
*
* At this point the pressures of the fluid state are correct.
*/
void zFractionUpdate_(const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
{
const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
unsigned pvtRegionIdx = priVars.pvtRegionIndex();
auto& fs = asImp_().fluidState_;
zFraction_ = priVars.makeEvaluation(zFractionIdx, timeIdx);
oilViscosity_ = ExtboModule::oilViscosity(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
gasViscosity_ = ExtboModule::gasViscosity(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
bo_ = ExtboModule::bo(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
bg_ = ExtboModule::bg(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
bz_ = ExtboModule::bg(pvtRegionIdx, fs.pressure(oilPhaseIdx), 0.99);
if (FluidSystem::enableDissolvedGas())
rs_ = ExtboModule::rs(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
else
rs_ = 0.0;
if (FluidSystem::enableVaporizedOil())
rv_ = ExtboModule::rv(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
else
rv_ = 0.0;
xVolume_ = ExtboModule::xVolume(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
yVolume_ = ExtboModule::yVolume(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
Evaluation pbub = fs.pressure(oilPhaseIdx);
if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
static const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
Scalar Sw = 0.0;
if (Indices::waterEnabled)
Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx).value();
if (Sw >= thresholdWaterFilledCell)
rs_ = 0.0; // water only, zero rs_ ...
}
if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
rs_ = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
const Evaluation zLim = ExtboModule::zLim(pvtRegionIdx);
if (zFraction_ > zLim) {
pbub = ExtboModule::pbubRs(pvtRegionIdx, zLim, rs_);
} else {
pbub = ExtboModule::pbubRs(pvtRegionIdx, zFraction_, rs_);
}
bo_ = ExtboModule::bo(pvtRegionIdx, pbub, zFraction_) + ExtboModule::oilCmp(pvtRegionIdx, zFraction_)*(fs.pressure(oilPhaseIdx)-pbub);
xVolume_ = ExtboModule::xVolume(pvtRegionIdx, pbub, zFraction_);
}
if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
rv_ = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
Evaluation rvsat = ExtboModule::rv(pvtRegionIdx, pbub, zFraction_);
bg_ = ExtboModule::bg(pvtRegionIdx, pbub, zFraction_) + ExtboModule::gasCmp(pvtRegionIdx, zFraction_)*(rv_-rvsat);
yVolume_ = ExtboModule::yVolume(pvtRegionIdx, pbub, zFraction_);
}
}
/*!
* \brief Re-compute face densities to account for zFraction dependency.
*
* At this point the pressures and saturations of the fluid state are correct.
*/
void zPvtUpdate_()
{
const auto& iq = asImp_();
auto& fs = asImp_().fluidState_;
unsigned pvtRegionIdx = iq.pvtRegionIndex();
zRefDensity_ = ExtboModule::referenceDensity(pvtRegionIdx);
fs.setInvB(oilPhaseIdx, 1.0/bo_);
fs.setInvB(gasPhaseIdx, 1.0/bg_);
fs.setDensity(oilPhaseIdx,
fs.invB(oilPhaseIdx)
*(FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx)
+ (1.0-xVolume_)*fs.Rs()*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx)
+ xVolume_*fs.Rs()*zRefDensity_ ));
fs.setDensity(gasPhaseIdx,
fs.invB(gasPhaseIdx)
*(FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx)*(1.0-yVolume_)+yVolume_*zRefDensity_
+ FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx)*fs.Rv()));
}
const Evaluation& zFraction() const
{ return zFraction_; }
const Evaluation& xVolume() const
{ return xVolume_; }
const Evaluation& yVolume() const
{ return yVolume_; }
const Evaluation& oilViscosity() const
{ return oilViscosity_; }
const Evaluation& gasViscosity() const
{ return gasViscosity_; }
const Evaluation& bo() const
{ return bo_; }
const Evaluation& bg() const
{ return bg_; }
const Evaluation& rs() const
{ return rs_; }
const Evaluation& rv() const
{ return rv_; }
const Evaluation zPureInvFormationVolumeFactor() const
{ return 1.0/bz_; }
const Scalar& zRefDensity() const
{ return zRefDensity_; }
private:
protected:
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
// Abstract "mass fraction" accounting for the solvent component. The relation between this
// quantity and the actual mass fraction of solvent, is implicitly defined from the specific
// pvt measurements as provided by kw PVTSOL.
Evaluation zFraction_;
// The solvent component is assumed gas at surface conditions
Evaluation xVolume_; // Solvent volume fraction of Rs
Evaluation yVolume_; // Solvent volume fraction of Sg/Bg
// Standard black oil parameters modified for presence of solvent
Evaluation oilViscosity_;
Evaluation gasViscosity_;
Evaluation bo_;
Evaluation bg_;
Evaluation rs_;
Evaluation rv_;
// Properties of pure solvent
Evaluation bz_;
Scalar zRefDensity_;
};
template <class TypeTag>
class BlackOilExtboIntensiveQuantities<TypeTag, false>
{
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
void zPvtUpdate_()
{ }
void zFractionUpdate_(const ElementContext& elemCtx OPM_UNUSED,
unsigned scvIdx OPM_UNUSED,
unsigned timeIdx OPM_UNUSED)
{ }
const Evaluation& xVolume() const
{ throw std::runtime_error("xVolume() called but extbo is disabled"); }
const Evaluation& yVolume() const
{ throw std::runtime_error("yVolume() called but extbo is disabled"); }
const Evaluation& oilViscosity() const
{ throw std::runtime_error("oilViscosity() called but extbo is disabled"); }
const Evaluation& gasViscosity() const
{ throw std::runtime_error("gasViscosity() called but extbo is disabled"); }
const Evaluation& rs() const
{ throw std::runtime_error("rs() called but extbo is disabled"); }
const Evaluation& rv() const
{ throw std::runtime_error("rv() called but extbo is disabled"); }
const Evaluation& zPureInvFormationVolumeFactor() const
{ throw std::runtime_error("zPureInvFormationVolumeFactor() called but extbo is disabled"); }
const Evaluation& zFraction() const
{ throw std::runtime_error("zFraction() called but extbo is disabled"); }
const Evaluation& zInverseFormationVolumeFactor() const
{ throw std::runtime_error("zInverseFormationVolumeFactor() called but extbo is disabled"); }
const Scalar& zRefDensity() const
{ throw std::runtime_error("zRefDensity() called but extbo is disabled"); }
};
/*!
* \ingroup BlackOil
* \class Opm::BlackOilExtboExtensiveQuantities
*
* \brief Provides the solvent specific extensive quantities to the generic black-oil
* module's extensive quantities.
*/
template <class TypeTag, bool enableExtboV = getPropValue<TypeTag, Properties::EnableExtbo>()>
class BlackOilExtboExtensiveQuantities
{
using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Toolbox = Opm::MathToolbox<Evaluation>;
static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
static constexpr int dimWorld = GridView::dimensionworld;
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
public:
private:
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
};
template <class TypeTag>
class BlackOilExtboExtensiveQuantities<TypeTag, false>
{
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
public:
};
} // namespace Opm
#endif

32
opm/models/blackoil/blackoilindices.hh
View File

@ -35,7 +35,7 @@ namespace Opm {
*
* \brief The primary variable and equation indices for the black-oil model.
*/
template <unsigned numSolventsV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset>
template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset>
struct BlackOilIndices
{
//! Number of phases active at all times
@ -49,6 +49,9 @@ struct BlackOilIndices
//! Are solvents involved?
static const bool enableSolvent = numSolventsV > 0;
//! Is extbo invoked?
static const bool enableExtbo = numExtbosV > 0;
//! Are polymers involved?
static const bool enablePolymer = numPolymersV > 0;
@ -58,6 +61,9 @@ struct BlackOilIndices
//! Number of solvent components to be considered
static const int numSolvents = enableSolvent ? numSolventsV : 0;
//! Number of components to be considered for extbo
static const int numExtbos = enableExtbo ? numExtbosV : 0;
//! Number of polymer components to be considered
static const int numPolymers = enablePolymer ? numPolymersV : 0;
@ -71,7 +77,7 @@ struct BlackOilIndices
static const int numBrine = enableBrine? 1 : 0;
//! The number of equations
static const int numEq = numPhases + numSolvents + numPolymers + numEnergy + numFoam + numBrine;
static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine;
//! \brief returns the index of "active" component
static constexpr unsigned canonicalToActiveComponentIndex(unsigned compIdx)
@ -104,6 +110,10 @@ struct BlackOilIndices
static const int solventSaturationIdx =
enableSolvent ? PVOffset + numPhases : -1000;
//! Index of the primary variable for the first extbo component
static const int zFractionIdx =
enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
//! Index of the primary variable for the first polymer
static const int polymerConcentrationIdx =
enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
@ -114,15 +124,15 @@ struct BlackOilIndices
//! Index of the primary variable for the foam
static const int foamConcentrationIdx =
enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers : -1000;
//! Index of the primary variable for the brine
static const int saltConcentrationIdx =
enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numExtbos + numPolymers + numFoam : -1000;
//! Index of the primary variable for temperature
static const int temperatureIdx =
enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine : - 1000;
enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
////////
@ -137,9 +147,13 @@ struct BlackOilIndices
static const int contiSolventEqIdx =
enableSolvent ? PVOffset + numPhases : -1000;
//! Index of the continuity equation for the first extbo component
static const int contiZfracEqIdx =
enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
//! Index of the continuity equation for the first polymer component
static const int contiPolymerEqIdx =
enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
enablePolymer ? PVOffset + numPhases + numSolvents + numExtbos : -1000;
//! Index of the continuity equation for the second polymer component (molecular weight)
static const int contiPolymerMWEqIdx =
@ -147,16 +161,16 @@ struct BlackOilIndices
//! Index of the continuity equation for the foam component
static const int contiFoamEqIdx =
enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers : -1000;
//! Index of the continuity equation for the salt water component
static const int contiBrineEqIdx =
enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000;
//! Index of the continuity equation for energy
static const int contiEnergyEqIdx =
enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine: -1000;
enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: -1000;
};
} // namespace Opm

23
opm/models/blackoil/blackoilintensivequantities.hh
View File

@ -30,6 +30,7 @@
#include "blackoilproperties.hh"
#include "blackoilsolventmodules.hh"
#include "blackoilextbomodules.hh"
#include "blackoilpolymermodules.hh"
#include "blackoilfoammodules.hh"
#include "blackoilbrinemodules.hh"
@ -55,6 +56,7 @@ class BlackOilIntensiveQuantities
: public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
, public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
, public BlackOilSolventIntensiveQuantities<TypeTag>
, public BlackOilExtboIntensiveQuantities<TypeTag>
, public BlackOilPolymerIntensiveQuantities<TypeTag>
, public BlackOilFoamIntensiveQuantities<TypeTag>
, public BlackOilBrineIntensiveQuantities<TypeTag>
@ -75,6 +77,7 @@ class BlackOilIntensiveQuantities
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
@ -210,6 +213,9 @@ public:
// update the Saturation functions for the blackoil solvent module.
asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
// update extBO parameters
asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
Evaluation SoMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
SoMax = Opm::max(fluidState_.saturation(oilPhaseIdx),
@ -223,7 +229,7 @@ public:
// we use the compositions of the gas-saturated oil and oil-saturated gas.
if (FluidSystem::enableDissolvedGas()) {
Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
const Evaluation& RsSat =
const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
oilPhaseIdx,
pvtRegionIdx,
@ -235,7 +241,7 @@ public:
if (FluidSystem::enableVaporizedOil()) {
Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
const Evaluation& RvSat =
const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx,
@ -257,7 +263,7 @@ public:
// the gas phase is not present, but we need to compute its "composition"
// for the gravity correction anyway
Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
const auto& RvSat =
const auto& RvSat = enableExtbo ? asImp_().rv() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
gasPhaseIdx,
pvtRegionIdx,
@ -276,7 +282,7 @@ public:
// the oil phase is not present, but we need to compute its "composition" for
// the gravity correction anyway
Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
const auto& RsSat =
const auto& RsSat = enableExtbo ? asImp_().rs() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
oilPhaseIdx,
pvtRegionIdx,
@ -306,7 +312,12 @@ public:
fluidState_.setInvB(phaseIdx, b);
const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
mobility_[phaseIdx] /= mu;
if (enableExtbo && phaseIdx == oilPhaseIdx)
mobility_[phaseIdx] /= asImp_().oilViscosity();
else if (enableExtbo && phaseIdx == gasPhaseIdx)
mobility_[phaseIdx] /= asImp_().gasViscosity();
else
mobility_[phaseIdx] /= mu;
}
Opm::Valgrind::CheckDefined(mobility_);
@ -366,6 +377,7 @@ public:
porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
asImp_().zPvtUpdate_();
asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
@ -446,6 +458,7 @@ public:
private:
friend BlackOilSolventIntensiveQuantities<TypeTag>;
friend BlackOilExtboIntensiveQuantities<TypeTag>;
friend BlackOilPolymerIntensiveQuantities<TypeTag>;
friend BlackOilEnergyIntensiveQuantities<TypeTag>;
friend BlackOilFoamIntensiveQuantities<TypeTag>;

8
opm/models/blackoil/blackoillocalresidual.hh
View File

@ -30,6 +30,7 @@
#include "blackoilproperties.hh"
#include "blackoilsolventmodules.hh"
#include "blackoilextbomodules.hh"
#include "blackoilpolymermodules.hh"
#include "blackoilenergymodules.hh"
#include "blackoilfoammodules.hh"
@ -80,6 +81,7 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
using Toolbox = Opm::MathToolbox<Evaluation>;
using SolventModule = BlackOilSolventModule<TypeTag>;
using ExtboModule = BlackOilExtboModule<TypeTag>;
using PolymerModule = BlackOilPolymerModule<TypeTag>;
using EnergyModule = BlackOilEnergyModule<TypeTag>;
using FoamModule = BlackOilFoamModule<TypeTag>;
@ -143,6 +145,9 @@ public:
// deal with solvents (if present)
SolventModule::addStorage(storage, intQuants);
// deal with zFracton (if present)
ExtboModule::addStorage(storage, intQuants);
// deal with polymer (if present)
PolymerModule::addStorage(storage, intQuants);
@ -186,6 +191,9 @@ public:
// deal with solvents (if present)
SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
// deal with zFracton (if present)
ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
// deal with polymer (if present)
PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);

19
opm/models/blackoil/blackoilmodel.hh
View File

@ -45,6 +45,7 @@
#include "blackoilpolymermodules.hh"
#include "blackoilfoammodules.hh"
#include "blackoilbrinemodules.hh"
#include "blackoilextbomodules.hh"
#include "blackoildarcyfluxmodule.hh"
#include <opm/models/common/multiphasebasemodel.hh>
@ -122,6 +123,7 @@ struct FluxModule<TypeTag, TTag::BlackOilModel> { using type = Opm::BlackOilDarc
template<class TypeTag>
struct Indices<TypeTag, TTag::BlackOilModel>
{ using type = Opm::BlackOilIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
getPropValue<TypeTag, Properties::EnableEnergy>(),
getPropValue<TypeTag, Properties::EnableFoam>(),
@ -144,6 +146,8 @@ public:
template<class TypeTag>
struct EnableSolvent<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
template<class TypeTag>
struct EnableExtbo<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
template<class TypeTag>
struct EnablePolymer<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
template<class TypeTag>
struct EnablePolymerMW<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
@ -271,6 +275,7 @@ class BlackOilModel
static const bool waterEnabled = Indices::waterEnabled;
using SolventModule = BlackOilSolventModule<TypeTag>;
using ExtboModule = BlackOilExtboModule<TypeTag>;
using PolymerModule = BlackOilPolymerModule<TypeTag>;
using EnergyModule = BlackOilEnergyModule<TypeTag>;
public:
@ -286,6 +291,7 @@ public:
ParentType::registerParameters();
SolventModule::registerParameters();
ExtboModule::registerParameters();
PolymerModule::registerParameters();
EnergyModule::registerParameters();
@ -315,6 +321,8 @@ public:
oss << "composition_switching";
else if (SolventModule::primaryVarApplies(pvIdx))
return SolventModule::primaryVarName(pvIdx);
else if (ExtboModule::primaryVarApplies(pvIdx))
return ExtboModule::primaryVarName(pvIdx);
else if (PolymerModule::primaryVarApplies(pvIdx))
return PolymerModule::primaryVarName(pvIdx);
else if (EnergyModule::primaryVarApplies(pvIdx))
@ -336,6 +344,8 @@ public:
oss << "conti_" << FluidSystem::phaseName(eqIdx - Indices::conti0EqIdx);
else if (SolventModule::eqApplies(eqIdx))
return SolventModule::eqName(eqIdx);
else if (ExtboModule::eqApplies(eqIdx))
return ExtboModule::eqName(eqIdx);
else if (PolymerModule::eqApplies(eqIdx))
return PolymerModule::eqName(eqIdx);
else if (EnergyModule::eqApplies(eqIdx))
@ -369,6 +379,10 @@ public:
else if (SolventModule::primaryVarApplies(pvIdx))
return SolventModule::primaryVarWeight(pvIdx);
// deal with primary variables stemming from the extBO module
else if (ExtboModule::primaryVarApplies(pvIdx))
return ExtboModule::primaryVarWeight(pvIdx);
// deal with primary variables stemming from the polymer module
else if (PolymerModule::primaryVarApplies(pvIdx))
return PolymerModule::primaryVarWeight(pvIdx);
@ -424,6 +438,9 @@ public:
if (SolventModule::eqApplies(eqIdx))
return SolventModule::eqWeight(eqIdx);
if (ExtboModule::eqApplies(eqIdx))
return ExtboModule::eqWeight(eqIdx);
else if (PolymerModule::eqApplies(eqIdx))
return PolymerModule::eqWeight(eqIdx);
@ -463,6 +480,7 @@ public:
outstream << priVars.pvtRegionIndex() << " ";
SolventModule::serializeEntity(*this, outstream, dof);
ExtboModule::serializeEntity(*this, outstream, dof);
PolymerModule::serializeEntity(*this, outstream, dof);
EnergyModule::serializeEntity(*this, outstream, dof);
}
@ -500,6 +518,7 @@ public:
throw std::runtime_error("Could not deserialize degree of freedom "+std::to_string(dofIdx));
SolventModule::deserializeEntity(*this, instream, dof);
ExtboModule::deserializeEntity(*this, instream, dof);
PolymerModule::deserializeEntity(*this, instream, dof);
EnergyModule::deserializeEntity(*this, instream, dof);

12
opm/models/blackoil/blackoilnewtonmethod.hh
View File

@ -212,6 +212,7 @@ protected:
const EqVector& currentResidual)
{
static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
static constexpr bool enableExtbo = Indices::zFractionIdx >= 0;
static constexpr bool enablePolymer = Indices::polymerConcentrationIdx >= 0;
static constexpr bool enablePolymerWeight = Indices::polymerMoleWeightIdx >= 0;
static constexpr bool enableEnergy = Indices::temperatureIdx >= 0;
@ -287,6 +288,13 @@ protected:
else if (enableSolvent && pvIdx == Indices::solventSaturationIdx)
// solvent saturation updates are also subject to the Appleyard chop
delta *= satAlpha;
else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
// z fraction updates are also subject to the Appleyard chop
if (delta > currentValue[Indices::zFractionIdx])
delta = currentValue[Indices::zFractionIdx];
if (delta < currentValue[Indices::zFractionIdx]-1.0)
delta = currentValue[Indices::zFractionIdx]-1.0;
}
else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
const double sign = delta >= 0. ? 1. : -1.;
// maximum change of polymer molecular weight, the unit is MDa.
@ -303,6 +311,10 @@ protected:
if (enableSolvent && pvIdx == Indices::solventSaturationIdx)
nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
// keep the z fraction between 0 and 1
if (enableExtbo && pvIdx == Indices::zFractionIdx)
nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
// keep the polymer concentration above 0
if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx)
nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);

26
opm/models/blackoil/blackoilonephaseindices.hh
View File

@ -37,7 +37,7 @@ namespace Opm {
*
* \brief The primary variable and equation indices for the black-oil model.
*/
template <unsigned numSolventsV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned canonicalCompIdx>
template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned canonicalCompIdx>
struct BlackOilOnePhaseIndices
{
//! Is phase enabled or not
@ -48,6 +48,9 @@ struct BlackOilOnePhaseIndices
//! Are solvents involved?
static const bool enableSolvent = numSolventsV > 0;
//! Is extbo invoked?
static const bool enableExtbo = numExtbosV > 0;
//! Are polymers involved?
static const bool enablePolymer = numPolymersV > 0;
@ -57,6 +60,9 @@ struct BlackOilOnePhaseIndices
//! Number of solvent components to be considered
static const int numSolvents = enableSolvent ? numSolventsV : 0;
//! Number of components to be considered for extbo
static const int numExtbos = enableExtbo ? numExtbosV : 0;
//! Number of polymer components to be considered
static const int numPolymers = enablePolymer ? numPolymersV : 0;
@ -73,7 +79,7 @@ struct BlackOilOnePhaseIndices
static const int numPhases = 1;
//! The number of equations
static const int numEq = numPhases + numSolvents + numPolymers + numEnergy + numFoam + numBrine;
static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine;
//////////////////////////////
// Primary variable indices
@ -97,6 +103,10 @@ struct BlackOilOnePhaseIndices
static const int solventSaturationIdx =
enableSolvent ? PVOffset + numPhases : -1000;
//! Index of the primary variable for the first extbo component
static const int zFractionIdx =
enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
//! Index of the primary variable for the first polymer
static const int polymerConcentrationIdx =
enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
@ -111,11 +121,11 @@ struct BlackOilOnePhaseIndices
//! Index of the primary variable for the salt
static const int saltConcentrationIdx =
enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000;
//! Index of the primary variable for temperature
static const int temperatureIdx =
enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine: - 1000;
enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: - 1000;
//////////////////////
// Equation indices
@ -148,6 +158,10 @@ struct BlackOilOnePhaseIndices
static const int contiSolventEqIdx =
enableSolvent ? PVOffset + numPhases : -1000;
//! Index of the continuity equation for the first extbo component
static const int contiZfracEqIdx =
enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
//! Index of the continuity equation for the first polymer component
static const int contiPolymerEqIdx =
enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
@ -162,11 +176,11 @@ struct BlackOilOnePhaseIndices
//! Index of the continuity equation for the salt component
static const int contiBrineEqIdx =
enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000;
//! Index of the continuity equation for energy
static const int contiEnergyEqIdx =
enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine: -1000;
enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: -1000;
};
} // namespace Opm

62
opm/models/blackoil/blackoilprimaryvariables.hh
View File

@ -30,6 +30,7 @@
#include "blackoilproperties.hh"
#include "blackoilsolventmodules.hh"
#include "blackoilextbomodules.hh"
#include "blackoilpolymermodules.hh"
#include "blackoilenergymodules.hh"
#include "blackoilfoammodules.hh"
@ -49,6 +50,9 @@ namespace Opm {
template <class TypeTag, bool enableSolvent>
class BlackOilSolventModule;
template <class TypeTag, bool enableExtbo>
class BlackOilExtboModule;
template <class TypeTag, bool enablePolymer>
class BlackOilPolymerModule;
@ -94,6 +98,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
// component indices from the fluid system
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
@ -105,6 +110,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
using Toolbox = typename Opm::MathToolbox<Evaluation>;
using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
using PolymerModule = BlackOilPolymerModule<TypeTag, enablePolymer>;
using EnergyModule = BlackOilEnergyModule<TypeTag, enableEnergy>;
using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
@ -402,7 +408,10 @@ public:
Scalar T = asImp_().temperature_();
Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
Scalar RsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
Scalar RsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
Scalar RsSat = enableExtbo ? ExtboModule::rs(pvtRegionIndex(),
po,
zFraction_())
: FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
T,
po,
So2,
@ -412,6 +421,7 @@ public:
if (compositionSwitchEnabled)
(*this)[Indices::compositionSwitchIdx] =
std::min(RsMax, RsSat);
return true;
}
@ -433,12 +443,14 @@ public:
Scalar T = asImp_().temperature_();
Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
Scalar RvSat =
FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
T,
pg,
Scalar(0),
SoMax);
Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
pg,
zFraction_())
: FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
T,
pg,
Scalar(0),
SoMax);
setPrimaryVarsMeaning(Sw_pg_Rv);
(*this)[Indices::pressureSwitchIdx] = pg;
if (compositionSwitchEnabled)
@ -473,12 +485,14 @@ public:
Scalar So = 1.0 - Sw - solventSaturation_();
Scalar SoMax = std::max(So, problem.maxOilSaturation(globalDofIdx));
Scalar RsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
Scalar RsSat =
FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
T,
po,
So,
SoMax);
Scalar RsSat = enableExtbo ? ExtboModule::rs(pvtRegionIndex(),
po,
zFraction_())
: FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
T,
po,
So,
SoMax);
Scalar Rs = (*this)[Indices::compositionSwitchIdx];
if (Rs > std::min(RsMax, RsSat*(1.0 + eps))) {
@ -530,12 +544,14 @@ public:
Scalar T = asImp_().temperature_();
Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
Scalar RvSat =
FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
T,
pg,
/*So=*/Scalar(0.0),
SoMax);
Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
pg,
zFraction_())
: FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
T,
pg,
/*So=*/Scalar(0.0),
SoMax);
Scalar Rv = (*this)[Indices::compositionSwitchIdx];
if (Rv > std::min(RvMax, RvSat*(1.0 + eps))) {
@ -686,6 +702,14 @@ private:
return (*this)[Indices::solventSaturationIdx];
}
Scalar zFraction_() const
{
if (!enableExtbo)
return 0.0;
return (*this)[Indices::zFractionIdx];
}
Scalar polymerConcentration_() const
{
if (!enablePolymer)

3
opm/models/blackoil/blackoilproperties.hh
View File

@ -40,6 +40,9 @@ struct EnableEclipseOutput { using type = UndefinedProperty; };
//! Enable the ECL-blackoil extension for solvents. ("Second gas")
template<class TypeTag, class MyTypeTag>
struct EnableSolvent { using type = UndefinedProperty; };
//! Enable the ECL-blackoil extension for extended BO. ("Second gas" - alternative approach)
template<class TypeTag, class MyTypeTag>
struct EnableExtbo { using type = UndefinedProperty; };
//! Enable the ECL-blackoil extension for polymer.
template<class TypeTag, class MyTypeTag>
struct EnablePolymer { using type = UndefinedProperty; };

22
opm/models/blackoil/blackoiltwophaseindices.hh
View File

@ -37,7 +37,7 @@ namespace Opm {
*
* \brief The primary variable and equation indices for the black-oil model.
*/
template <unsigned numSolventsV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned disabledCanonicalCompIdx>
template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned disabledCanonicalCompIdx>
struct BlackOilTwoPhaseIndices
{
//! Is phase enabled or not
@ -48,6 +48,9 @@ struct BlackOilTwoPhaseIndices
//! Are solvents involved?
static const bool enableSolvent = numSolventsV > 0;
//! Is extbo invoked?
static const bool enableExtbo = numExtbosV > 0;
//! Are polymers involved?
static const bool enablePolymer = numPolymersV > 0;
@ -57,6 +60,9 @@ struct BlackOilTwoPhaseIndices
//! Number of solvent components to be considered
static const int numSolvents = enableSolvent ? numSolventsV : 0;
//! Number of components to be considered for extbo
static const int numExtbos = enableExtbo ? numExtbosV : 0;
//! Number of polymer components to be considered
static const int numPolymers = enablePolymer ? numPolymersV : 0;
@ -73,7 +79,7 @@ struct BlackOilTwoPhaseIndices
static const int numPhases = 2;
//! The number of equations
static const int numEq = numPhases + numSolvents + numPolymers + numEnergy + numFoam + numBrine;
static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine;
//////////////////////////////
// Primary variable indices
@ -97,6 +103,10 @@ struct BlackOilTwoPhaseIndices
static const int solventSaturationIdx =
enableSolvent ? PVOffset + numPhases : -1000;
//! Index of the primary variable for the first extbo component
static const int zFractionIdx =
enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
//! Index of the primary variable for the first polymer
static const int polymerConcentrationIdx =
enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
@ -115,7 +125,7 @@ struct BlackOilTwoPhaseIndices
//! Index of the primary variable for temperature
static const int temperatureIdx =
enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine : - 1000;
enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
//////////////////////
// Equation indices
@ -166,6 +176,10 @@ struct BlackOilTwoPhaseIndices
static const int contiSolventEqIdx =
enableSolvent ? PVOffset + numPhases : -1000;
//! Index of the continuity equation for the first extbo component
static const int contiZfracEqIdx =
enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
//! Index of the continuity equation for the first polymer component
static const int contiPolymerEqIdx =
enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
@ -184,7 +198,7 @@ struct BlackOilTwoPhaseIndices
//! Index of the continuity equation for energy
static const int contiEnergyEqIdx =
enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine : -1000;
enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : -1000;
};
} // namespace Opm