Loading...
Searching...
No Matches
blackoilprimaryvariables.hh
Go to the documentation of this file.
// -*- 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.
*/
#ifndef EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH
#define EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH
#include <dune/common/fvector.hh>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
#include <opm/models/blackoil/blackoilbioeffectsmodules.hh>
#include <opm/models/blackoil/blackoilbrinemodules.hh>
#include <opm/models/blackoil/blackoilenergymodules.hh>
#include <opm/models/blackoil/blackoilextbomodules.hh>
#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/blackoil/blackoilsolventmodules.hh>
#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <stdexcept>
#include <type_traits>
namespace Opm::Parameters {
template<class Scalar>
struct PressureScale
{ static constexpr Scalar value = 1.0; };
} // namespace Opm::Parameters
namespace Opm {
template <class TypeTag>
class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
{
using ParentType = FvBasePrimaryVariables<TypeTag>;
using Implementation = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
// number of equations
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
// primary variable indices
enum { waterSwitchIdx = Indices::waterSwitchIdx };
enum { pressureSwitchIdx = Indices::pressureSwitchIdx };
enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
enum { saltConcentrationIdx = Indices::saltConcentrationIdx };
enum { solventSaturationIdx = Indices::solventSaturationIdx };
static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
static constexpr bool waterEnabled = Indices::waterEnabled;
static constexpr bool gasEnabled = Indices::gasEnabled;
static constexpr bool oilEnabled = Indices::oilEnabled;
// phase indices from the fluid system
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
// 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>() };
enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
enum { enableMICP = Indices::enableMICP };
enum { gasCompIdx = FluidSystem::gasCompIdx };
enum { waterCompIdx = FluidSystem::waterCompIdx };
enum { oilCompIdx = FluidSystem::oilCompIdx };
using Toolbox = MathToolbox<Evaluation>;
using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
using EnergyModule = BlackOilEnergyModule<TypeTag, enableEnergy>;
using BrineModule = BlackOilBrineModule<TypeTag, enableBrine>;
using BioeffectsModule = BlackOilBioeffectsModule<TypeTag>;
static_assert(numPhases == 3, "The black-oil model assumes three phases!");
static_assert(numComponents == 3, "The black-oil model assumes three components!");
public:
enum class WaterMeaning : std::uint8_t {
Sw, // water saturation
Rvw, // vaporized water
Rsw, // dissolved gas in water
Disabled, // The primary variable is not used
};
enum class PressureMeaning : std::uint8_t {
Po, // oil pressure
Pg, // gas pressure
Pw, // water pressure
};
enum class GasMeaning : std::uint8_t {
Sg, // gas saturation
Rs, // dissolved gas in oil
Rv, // vapporized oil
Disabled, // The primary variable is not used
};
enum class BrineMeaning : std::uint8_t {
Cs, // salt concentration
Sp, // (precipitated) salt saturation
Disabled, // The primary variable is not used
};
enum class SolventMeaning : std::uint8_t {
Ss, // solvent saturation
Rsolw, // dissolved solvent in water
Disabled, // The primary variable is not used
};
BlackOilPrimaryVariables()
{
Valgrind::SetUndefined(*this);
pvtRegionIdx_ = 0;
}
BlackOilPrimaryVariables(const BlackOilPrimaryVariables& value) = default;
static BlackOilPrimaryVariables serializationTestObject()
{
BlackOilPrimaryVariables result;
result.pvtRegionIdx_ = 1;
result.primaryVarsMeaningBrine_ = BrineMeaning::Sp;
result.primaryVarsMeaningGas_ = GasMeaning::Rv;
result.primaryVarsMeaningPressure_ = PressureMeaning::Pg;
result.primaryVarsMeaningWater_ = WaterMeaning::Rsw;
result.primaryVarsMeaningSolvent_ = SolventMeaning::Ss;
for (std::size_t i = 0; i < result.size(); ++i) {
result[i] = i + 1;
}
return result;
}
static void init()
{
// TODO: these parameters have undocumented non-trivial dependencies
pressureScale_ = Parameters::Get<Parameters::PressureScale<Scalar>>();
}
static void registerParameters()
{
Parameters::Register<Parameters::PressureScale<Scalar>>
("Scaling of pressure primary variable");
}
void setPressureScale(Scalar val)
{ pressureScale_ = val; }
Evaluation
makeEvaluation(unsigned varIdx, unsigned timeIdx,
LinearizationType linearizationType = LinearizationType()) const
{
const Scalar scale = varIdx == pressureSwitchIdx ? this->pressureScale_ : Scalar{1.0};
if (std::is_same_v<Evaluation, Scalar>) {
return (*this)[varIdx] * scale; // finite differences
}
else {
// automatic differentiation
if (timeIdx == linearizationType.time) {
return Toolbox::createVariable((*this)[varIdx], varIdx) * scale;
}
else {
return Toolbox::createConstant((*this)[varIdx]) * scale;
}
}
}
void setPvtRegionIndex(unsigned value)
{ pvtRegionIdx_ = static_cast<unsigned short>(value); }
unsigned pvtRegionIndex() const
{ return pvtRegionIdx_; }
WaterMeaning primaryVarsMeaningWater() const
{ return primaryVarsMeaningWater_; }
void setPrimaryVarsMeaningWater(WaterMeaning newMeaning)
{ primaryVarsMeaningWater_ = newMeaning; }
PressureMeaning primaryVarsMeaningPressure() const
{ return primaryVarsMeaningPressure_; }
void setPrimaryVarsMeaningPressure(PressureMeaning newMeaning)
{ primaryVarsMeaningPressure_ = newMeaning; }
GasMeaning primaryVarsMeaningGas() const
{ return primaryVarsMeaningGas_; }
void setPrimaryVarsMeaningGas(GasMeaning newMeaning)
{ primaryVarsMeaningGas_ = newMeaning; }
BrineMeaning primaryVarsMeaningBrine() const
{ return primaryVarsMeaningBrine_; }
void setPrimaryVarsMeaningBrine(BrineMeaning newMeaning)
{ primaryVarsMeaningBrine_ = newMeaning; }
SolventMeaning primaryVarsMeaningSolvent() const
{ return primaryVarsMeaningSolvent_; }
void setPrimaryVarsMeaningSolvent(SolventMeaning newMeaning)
{ primaryVarsMeaningSolvent_ = newMeaning; }
template <class FluidState>
void assignNaive(const FluidState& fluidState)
{
using ConstEvaluation = std::remove_reference_t<typename FluidState::Scalar>;
using FsEvaluation = std::remove_const_t<ConstEvaluation>;
using FsToolbox = MathToolbox<FsEvaluation>;
const bool gasPresent =
FluidSystem::phaseIsActive(gasPhaseIdx)
? fluidState.saturation(gasPhaseIdx) > 0.0
: false;
const bool oilPresent =
FluidSystem::phaseIsActive(oilPhaseIdx)
? fluidState.saturation(oilPhaseIdx) > 0.0
: false;
const bool waterPresent =
FluidSystem::phaseIsActive(waterPhaseIdx)
? fluidState.saturation(waterPhaseIdx) > 0.0
: false;
const auto& saltSaturation =
BlackOil::getSaltSaturation_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
const bool precipitatedSaltPresent = enableSaltPrecipitation ? saltSaturation > 0.0 : false;
const bool oneActivePhases = FluidSystem::numActivePhases() == 1;
// deal with the primary variables for the energy extension
EnergyModule::assignPrimaryVars(*this, fluidState);
// Determine the meaning of the pressure primary variables
// Depending on the phases present, this variable is either interpreted as the
// pressure of the oil phase, gas phase (if no oil) or water phase (if only water)
if (gasPresent && FluidSystem::enableVaporizedOil() && !oilPresent) {
primaryVarsMeaningPressure_ = PressureMeaning::Pg;
}
else if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
primaryVarsMeaningPressure_ = PressureMeaning::Po;
}
else if (waterPresent && FluidSystem::enableDissolvedGasInWater() && !gasPresent) {
primaryVarsMeaningPressure_ = PressureMeaning::Pw;
}
else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
primaryVarsMeaningPressure_ = PressureMeaning::Pg;
}
else {
assert(FluidSystem::phaseIsActive(waterPhaseIdx));
primaryVarsMeaningPressure_ = PressureMeaning::Pw;
}
// Determine the meaning of the water primary variables
// Depending on the phases present, this variable is either interpreted as
// water saturation or vapporized water in the gas phase
// For two-phase gas-oil models and one-phase case the variable is disabled.
if (waterPresent && gasPresent) {
primaryVarsMeaningWater_ = WaterMeaning::Sw;
}
else if (gasPresent && FluidSystem::enableVaporizedWater()) {
primaryVarsMeaningWater_ = WaterMeaning::Rvw;
}
else if (waterPresent && FluidSystem::enableDissolvedGasInWater()) {
primaryVarsMeaningWater_ = WaterMeaning::Rsw;
}
else if (FluidSystem::phaseIsActive(waterPhaseIdx) && !oneActivePhases) {
primaryVarsMeaningWater_ = WaterMeaning::Sw;
}
else {
primaryVarsMeaningWater_ = WaterMeaning::Disabled;
}
// Determine the meaning of the gas primary variables
// Depending on the phases present, this variable is either interpreted as the
// saturation of the gas phase, as the fraction of the gas component in the oil
// phase (Rs) or as the fraction of the oil component (Rv) in the gas phase.
// For two-phase water-oil and water-gas models and one-phase case the variable is disabled.
if (gasPresent && oilPresent) {
primaryVarsMeaningGas_ = GasMeaning::Sg;
}
else if (oilPresent && FluidSystem::enableDissolvedGas()) {
primaryVarsMeaningGas_ = GasMeaning::Rs;
}
else if (gasPresent && FluidSystem::enableVaporizedOil()) {
primaryVarsMeaningGas_ = GasMeaning::Rv;
}
else if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx)) {
primaryVarsMeaningGas_ = GasMeaning::Sg;
}
else {
primaryVarsMeaningGas_ = GasMeaning::Disabled;
}
// Determine the meaning of the brine primary variables
if constexpr (enableSaltPrecipitation) {
if (precipitatedSaltPresent) {
primaryVarsMeaningBrine_ = BrineMeaning::Sp;
}
else {
primaryVarsMeaningBrine_ = BrineMeaning::Cs;
}
}
else {
primaryVarsMeaningBrine_ = BrineMeaning::Disabled;
}
// assign the actual primary variables
switch (primaryVarsMeaningPressure()) {
case PressureMeaning::Po:
this->setScaledPressure_(FsToolbox::value(fluidState.pressure(oilPhaseIdx)));
break;
case PressureMeaning::Pg:
this->setScaledPressure_(FsToolbox::value(fluidState.pressure(gasPhaseIdx)));
break;
case PressureMeaning::Pw:
this->setScaledPressure_(FsToolbox::value(fluidState.pressure(waterPhaseIdx)));
break;
default:
throw std::logic_error("No valid primary variable selected for pressure");
}
switch (primaryVarsMeaningWater()) {
case WaterMeaning::Sw:
{
(*this)[waterSwitchIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
break;
}
case WaterMeaning::Rvw:
{
const auto& rvw = BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
(*this)[waterSwitchIdx] = rvw;
break;
}
case WaterMeaning::Rsw:
{
const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
(*this)[waterSwitchIdx] = Rsw;
break;
}
case WaterMeaning::Disabled:
break;
default:
throw std::logic_error("No valid primary variable selected for water");
}
switch (primaryVarsMeaningGas()) {
case GasMeaning::Sg:
(*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
break;
case GasMeaning::Rs:
{
const auto& rs = BlackOil::getRs_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
(*this)[compositionSwitchIdx] = rs;
break;
}
case GasMeaning::Rv:
{
const auto& rv = BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
(*this)[compositionSwitchIdx] = rv;
break;
}
case GasMeaning::Disabled:
break;
default:
throw std::logic_error("No valid primary variable selected for composision");
}
}
bool adaptPrimaryVariables(const Problem& problem,
unsigned globalDofIdx,
[[maybe_unused]] Scalar swMaximum,
Scalar thresholdWaterFilledCell, Scalar eps = 0.0)
{
// this function accesses quite a few black-oil specific low-level functions
// directly for better performance (instead of going the canonical way through
// the IntensiveQuantities). The reason is that most intensive quantities are not
// required to be able to decide if the primary variables needs to be switched or
// not, so it would be a waste to compute them.
// Both the primary variable meaning of water and gas are disabled i.e.
// It is a one-phase case and we no variable meaning switch is needed.
if (primaryVarsMeaningWater() == WaterMeaning::Disabled &&
primaryVarsMeaningGas() == GasMeaning::Disabled)
{
return false;
}
// Read the current saturation from the primary variables
Scalar sw = 0.0;
Scalar sg = 0.0;
Scalar saltConcentration = 0.0;
const Scalar& T = asImp_().temperature_(problem, globalDofIdx);
if (primaryVarsMeaningWater() == WaterMeaning::Sw) {
sw = (*this)[waterSwitchIdx];
}
if (primaryVarsMeaningGas() == GasMeaning::Sg) {
sg = (*this)[compositionSwitchIdx];
}
if (primaryVarsMeaningWater() == WaterMeaning::Rsw) {
sw = 1.0;
}
if (primaryVarsMeaningGas() == GasMeaning::Disabled && gasEnabled) {
sg = 1.0 - sw; // water + gas case
}
// if solid phase disappeares: Sp (Solid salt saturation) -> Cs (salt concentration)
// if solid phase appears: Cs (salt concentration) -> Sp (Solid salt saturation)
if constexpr (enableSaltPrecipitation) {
const Scalar saltSolubility = BrineModule::saltSol(pvtRegionIndex());
if (primaryVarsMeaningBrine() == BrineMeaning::Sp) {
saltConcentration = saltSolubility;
const Scalar saltSat = (*this)[saltConcentrationIdx];
if (saltSat < -eps) { // precipitated salt dissappears
setPrimaryVarsMeaningBrine(BrineMeaning::Cs);
(*this)[saltConcentrationIdx] = saltSolubility; // set salt concentration to solubility limit
}
}
else if (primaryVarsMeaningBrine() == BrineMeaning::Cs) {
saltConcentration = (*this)[saltConcentrationIdx];
if (saltConcentration > saltSolubility + eps) { // salt concentration exceeds solubility limit
setPrimaryVarsMeaningBrine(BrineMeaning::Sp);
(*this)[saltConcentrationIdx] = 0.0;
}
}
}
// if solvent saturation disappeares: Ss (Solvent saturation) -> Rsolw (solvent dissolved in water)
// if solvent saturation appears: Rsolw (solvent dissolved in water) -> Ss (Solvent saturation)
// Scalar rsolw = 0.0; // not needed at the moment since we dont allow for vapwat in combination with rsolw
if constexpr (enableSolvent) {
if (SolventModule::isSolubleInWater()) {
const Scalar p = (*this)[pressureSwitchIdx]; // cap-pressure?
const Scalar solLimit =
SolventModule::solubilityLimit(pvtRegionIndex(), T , p, saltConcentration);
if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
const Scalar solSat = (*this)[solventSaturationIdx];
if (solSat < -eps) { // solvent dissappears
setPrimaryVarsMeaningSolvent(SolventMeaning::Rsolw);
(*this)[solventSaturationIdx] = solLimit; // set rsolw to solubility limit
}
}
else if (primaryVarsMeaningSolvent() == SolventMeaning::Rsolw) {
const Scalar rsolw = (*this)[solventSaturationIdx];
if (rsolw > solLimit + eps) { // solvent appears as phase
setPrimaryVarsMeaningSolvent(SolventMeaning::Ss);
(*this)[solventSaturationIdx] = 0.0;
}
}
}
}
// keep track if any meaning has changed
bool changed = false;
// Special case for cells with almost only water
// for these cells both saturations (if the phase is enabled) is used
// to avoid singular systems.
// If dissolved gas in water is enabled we shouldn't enter
// here but instead switch to Rsw as primary variable
// as sw >= 1.0 -> gas <= 0 (i.e. gas phase disappears)
if (sw >= thresholdWaterFilledCell && !FluidSystem::enableDissolvedGasInWater()) {
// make sure water saturations does not exceed sw_maximum. Default to 1.0
if constexpr (waterEnabled) {
(*this)[Indices::waterSwitchIdx] = std::min(swMaximum, sw);
assert(primaryVarsMeaningWater() == WaterMeaning::Sw);
}
// the hydrocarbon gas saturation is set to 0.0
if constexpr (compositionSwitchEnabled) {
(*this)[Indices::compositionSwitchIdx] = 0.0;
}
changed = primaryVarsMeaningGas() != GasMeaning::Sg;
if (changed) {
if constexpr (compositionSwitchEnabled) {
setPrimaryVarsMeaningGas(GasMeaning::Sg);
}
// use water pressure?
}
return changed;
}
pcFactor_ = 1.0;
if constexpr (enableBrine) {
if (BrineModule::hasPcfactTables() && primaryVarsMeaningBrine() == BrineMeaning::Sp) {
unsigned satnumRegionIdx = problem.satnumRegionIndex(globalDofIdx);
Scalar Sp = saltConcentration_();
Scalar porosityFactor = std::min(1.0 - Sp, 1.0); //phi/phi_0
const auto& pcfactTable = BrineModule::pcfactTable(satnumRegionIdx);
pcFactor_ = pcfactTable.eval(porosityFactor, /*extrapolation=*/true);
}
}
else if constexpr (enableBioeffects) {
if (BioeffectsModule::hasPcfactTables() && problem.referencePorosity(globalDofIdx, 0) > 0) {
unsigned satnumRegionIdx = problem.satnumRegionIndex(globalDofIdx);
Scalar Sb = biofilmVolumeFraction_() /
problem.referencePorosity(globalDofIdx, 0);
Scalar porosityFactor = std::min(1.0 - Sb, 1.0); //phi/phi_0
const auto& pcfactTable = BioeffectsModule::pcfactTable(satnumRegionIdx);
pcFactor_ = pcfactTable.eval(porosityFactor, /*extrapolation=*/true);
}
}
switch (primaryVarsMeaningWater()) {
case WaterMeaning::Sw:
{
// if water phase disappeares: Sw (water saturation) -> Rvw (fraction of water in gas phase)
if (sw < -eps && sg > eps && FluidSystem::enableVaporizedWater()) {
Scalar p = this->pressure_();
if (primaryVarsMeaningPressure() == PressureMeaning::Po) {
std::array<Scalar, numPhases> pC{};
const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
const Scalar so = 1.0 - sg - solventSaturation_();
computeCapillaryPressures_(pC, so, sg + solventSaturation_(), /*sw=*/ 0.0, matParams);
p += pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
}
const Scalar rvwSat =
FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
T,
p,
saltConcentration);
setPrimaryVarsMeaningWater(WaterMeaning::Rvw);
(*this)[Indices::waterSwitchIdx] = rvwSat; // primary variable becomes Rvw
changed = true;
break;
}
// if gas phase disappeares: Sw (water saturation) -> Rsw (fraction of gas in water phase)
// and Pg (gas pressure) -> Pw ( water pressure)
if (sg < -eps && sw > eps && FluidSystem::enableDissolvedGasInWater()) {
const Scalar pg = this->pressure_();
assert(primaryVarsMeaningPressure() == PressureMeaning::Pg);
std::array<Scalar, numPhases> pC = { 0.0 };
const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
const Scalar so = 1.0 - sw - solventSaturation_();
computeCapillaryPressures_(pC, so, /*sg=*/ 0.0, sw, matParams);
const Scalar pw = pg + pcFactor_ * (pC[waterPhaseIdx] - pC[gasPhaseIdx]);
const Scalar rswSat =
FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
T,
pw,
saltConcentration);
setPrimaryVarsMeaningWater(WaterMeaning::Rsw);
const Scalar rswMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
(*this)[Indices::waterSwitchIdx] = std::min(rswSat, rswMax); //primary variable becomes Rsw
setPrimaryVarsMeaningPressure(PressureMeaning::Pw);
this->setScaledPressure_(pw);
changed = true;
break;
}
break;
}
case WaterMeaning::Rvw:
{
const Scalar& rvw = (*this)[waterSwitchIdx];
Scalar p = this->pressure_();
if (primaryVarsMeaningPressure() == PressureMeaning::Po) {
std::array<Scalar, numPhases> pC{};
const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
const Scalar so = 1.0 - sg - solventSaturation_();
computeCapillaryPressures_(pC, so, sg + solventSaturation_(), /*sw=*/ 0.0, matParams);
p += pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
}
const Scalar rvwSat =
FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
T,
p,
saltConcentration);
// if water phase appears: Rvw (fraction of water in gas phase) -> Sw (water saturation)
if (rvw > rvwSat * (1.0 + eps)) {
setPrimaryVarsMeaningWater(WaterMeaning::Sw);
(*this)[Indices::waterSwitchIdx] = 0.0; // water saturation
changed = true;
}
break;
}
case WaterMeaning::Rsw:
{
// Gas phase not present. The hydrocarbon gas phase
// appears as soon as more of the gas component is present in the water phase
// than what saturated water can hold.
const Scalar& pw = this->pressure_();
assert(primaryVarsMeaningPressure() == PressureMeaning::Pw);
const Scalar rswSat =
FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
T,
pw,
saltConcentration);
const Scalar rsw = (*this)[Indices::waterSwitchIdx];
const Scalar rswMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
if (rsw > std::min(rswSat, rswMax)) {
// the gas phase appears, i.e., switch the primary variables to WaterMeaning::Sw
setPrimaryVarsMeaningWater(WaterMeaning::Sw);
(*this)[Indices::waterSwitchIdx] = 1.0; // hydrocarbon water saturation
setPrimaryVarsMeaningPressure(PressureMeaning::Pg);
std::array<Scalar, numPhases> pC{};
const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
computeCapillaryPressures_(pC, /*so=*/ 0.0, /*sg=*/ 0.0, /*sw=*/ 1.0, matParams);
const Scalar pg = pw + pcFactor_ * (pC[gasPhaseIdx] - pC[waterPhaseIdx]);
this->setScaledPressure_(pg);
changed = true;
}
break;
}
case WaterMeaning::Disabled:
break;
default:
throw std::logic_error("No valid primary variable selected for water");
}
// if gas phase disappeares: Sg (gas saturation) -> Rs (fraction of gas in oil phase)
// if oil phase disappeares: Sg (gas saturation) -> Rv (fraction of oil in gas phase)
// Po (oil pressure ) -> Pg (gas pressure)
// if gas phase appears: Rs (fraction of gas in oil phase) -> Sg (gas saturation)
// if oil phase appears: Rv (fraction of oil in gas phase) -> Sg (gas saturation)
// Pg (gas pressure ) -> Po (oil pressure)
switch (primaryVarsMeaningGas()) {
case GasMeaning::Sg:
{
const Scalar s = 1.0 - sw - solventSaturation_();
if (sg < -eps && s > 0.0 && FluidSystem::enableDissolvedGas()) {
const Scalar po = this->pressure_();
setPrimaryVarsMeaningGas(GasMeaning::Rs);
const Scalar soMax = std::max(s, problem.maxOilSaturation(globalDofIdx));
const Scalar rsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
const Scalar rsSat =
enableExtbo
? ExtboModule::rs(pvtRegionIndex(), po, zFraction_())
: FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
T,
po,
s,
soMax);
(*this)[Indices::compositionSwitchIdx] = std::min(rsMax, rsSat);
changed = true;
}
const Scalar so = 1.0 - sw - solventSaturation_() - sg;
if (so < -eps && sg > 0.0 && FluidSystem::enableVaporizedOil()) {
// the oil phase disappeared and some hydrocarbon gas phase is still
// present, i.e., switch the primary variables to GasMeaning::Rv.
// we only have the oil pressure readily available, but we need the gas
// pressure, i.e. we must determine capillary pressure
const Scalar po = this->pressure_();
std::array<Scalar, numPhases> pC{};
const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
computeCapillaryPressures_(pC, /*so=*/0.0, sg + solventSaturation_(), sw, matParams);
const Scalar pg = po + pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
// we start at the GasMeaning::Rv value that corresponds to that of oil-saturated
// hydrocarbon gas
setPrimaryVarsMeaningPressure(PressureMeaning::Pg);
this->setScaledPressure_(pg);
const Scalar soMax = problem.maxOilSaturation(globalDofIdx);
const Scalar rvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
const Scalar rvSat =
enableExtbo
? ExtboModule::rv(pvtRegionIndex(), pg, zFraction_())
: FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
T,
pg,
Scalar(0),
soMax);
setPrimaryVarsMeaningGas(GasMeaning::Rv);
(*this)[Indices::compositionSwitchIdx] = std::min(rvMax, rvSat);
changed = true;
}
break;
}
case GasMeaning::Rs:
{
// Gas phase not present. The hydrocarbon gas phase
// appears as soon as more of the gas component is present in the oil phase
// than what saturated oil can hold.
const Scalar po = this->pressure_();
const Scalar so = 1.0 - sw - solventSaturation_();
const Scalar soMax = std::max(so, problem.maxOilSaturation(globalDofIdx));
const Scalar rsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
const Scalar rsSat =
enableExtbo
? ExtboModule::rs(pvtRegionIndex(), po, zFraction_())
: FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
T,
po,
so,
soMax);
const Scalar rs = (*this)[Indices::compositionSwitchIdx];
if (rs > std::min(rsMax, rsSat * (Scalar{1.0} + eps))) {
// the gas phase appears, i.e., switch the primary variables to GasMeaning::Sg
setPrimaryVarsMeaningGas(GasMeaning::Sg);
(*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
changed = true;
}
break;
}
case GasMeaning::Rv:
{
// The oil phase appears as
// soon as more of the oil component is present in the hydrocarbon gas phase
// than what saturated gas contains. Note that we use the blackoil specific
// low-level PVT objects here for performance reasons.
const Scalar pg = this->pressure_();
const Scalar soMax = problem.maxOilSaturation(globalDofIdx);
const Scalar rvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
const Scalar rvSat =
enableExtbo
? ExtboModule::rv(pvtRegionIndex(), pg, zFraction_())
: FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
T,
pg,
/*so=*/Scalar(0.0),
soMax);
const Scalar rv = (*this)[Indices::compositionSwitchIdx];
if (rv > std::min(rvMax, rvSat * (Scalar{1.0} + eps))) {
// switch to phase equilibrium mode because the oil phase appears. here
// we also need the capillary pressures to calculate the oil phase
// pressure using the gas phase pressure
const Scalar sg2 = 1.0 - sw - solventSaturation_();
std::array<Scalar, numPhases> pC{};
const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
computeCapillaryPressures_(pC,
/*so=*/0.0,
/*sg=*/sg2 + solventSaturation_(),
sw,
matParams);
const Scalar po = pg + pcFactor_ * (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
setPrimaryVarsMeaningGas(GasMeaning::Sg);
setPrimaryVarsMeaningPressure(PressureMeaning::Po);
this->setScaledPressure_(po);
(*this)[Indices::compositionSwitchIdx] = sg2; // hydrocarbon gas saturation
changed = true;
}
break;
}
case GasMeaning::Disabled:
break;
default:
throw std::logic_error("No valid primary variable selected for water");
}
return changed;
}
bool chopAndNormalizeSaturations()
{
if (primaryVarsMeaningWater() == WaterMeaning::Disabled &&
primaryVarsMeaningGas() == GasMeaning::Disabled)
{
return false;
}
Scalar sw = 0.0;
if (primaryVarsMeaningWater() == WaterMeaning::Sw) {
sw = (*this)[Indices::waterSwitchIdx];
}
Scalar sg = 0.0;
if (primaryVarsMeaningGas() == GasMeaning::Sg) {
sg = (*this)[Indices::compositionSwitchIdx];
}
Scalar ssol = 0.0;
if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
ssol =(*this) [Indices::solventSaturationIdx];
}
Scalar so = 1.0 - sw - sg - ssol;
sw = std::min(std::max(sw, Scalar{0.0}), Scalar{1.0});
so = std::min(std::max(so, Scalar{0.0}), Scalar{1.0});
sg = std::min(std::max(sg, Scalar{0.0}), Scalar{1.0});
ssol = std::min(std::max(ssol, Scalar{0.0}), Scalar{1.0});
const Scalar st = sw + so + sg + ssol;
sw = sw / st;
sg = sg / st;
ssol = ssol / st;
assert(st > 0.5);
if (primaryVarsMeaningWater() == WaterMeaning::Sw) {
(*this)[Indices::waterSwitchIdx] = sw;
}
if (primaryVarsMeaningGas() == GasMeaning::Sg) {
(*this)[Indices::compositionSwitchIdx] = sg;
}
if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
(*this) [Indices::solventSaturationIdx] = ssol;
}
return (st != 1);
}
BlackOilPrimaryVariables& operator=(const BlackOilPrimaryVariables& other) = default;
using ParentType::operator=;
void checkDefined() const
{
#ifndef NDEBUG
// check the "real" primary variables
for (unsigned i = 0; i < this->size(); ++i) {
Valgrind::CheckDefined((*this)[i]);
}
// check the "pseudo" primary variables
Valgrind::CheckDefined(primaryVarsMeaningWater_);
Valgrind::CheckDefined(primaryVarsMeaningGas_);
Valgrind::CheckDefined(primaryVarsMeaningPressure_);
Valgrind::CheckDefined(primaryVarsMeaningBrine_);
Valgrind::CheckDefined(primaryVarsMeaningSolvent_);
Valgrind::CheckDefined(pvtRegionIdx_);
#endif // NDEBUG
}
template<class Serializer>
void serializeOp(Serializer& serializer)
{
using FV = Dune::FieldVector<Scalar, getPropValue<TypeTag, Properties::NumEq>()>;
serializer(static_cast<FV&>(*this));
serializer(primaryVarsMeaningWater_);
serializer(primaryVarsMeaningPressure_);
serializer(primaryVarsMeaningGas_);
serializer(primaryVarsMeaningBrine_);
serializer(primaryVarsMeaningSolvent_);
serializer(pvtRegionIdx_);
}
bool operator==(const BlackOilPrimaryVariables& rhs) const
{
return
static_cast<const FvBasePrimaryVariables<TypeTag>&>(*this) == rhs
&& this->primaryVarsMeaningWater_ == rhs.primaryVarsMeaningWater_
&& this->primaryVarsMeaningPressure_ == rhs.primaryVarsMeaningPressure_
&& this->primaryVarsMeaningGas_ == rhs.primaryVarsMeaningGas_
&& this->primaryVarsMeaningBrine_ == rhs.primaryVarsMeaningBrine_
&& this->primaryVarsMeaningSolvent_ == rhs.primaryVarsMeaningSolvent_
&& this->pvtRegionIdx_ == rhs.pvtRegionIdx_;
}
private:
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation& asImp_() const
{ return *static_cast<const Implementation*>(this); }
Scalar solventSaturation_() const
{
if constexpr (enableSolvent) {
if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
return (*this)[Indices::solventSaturationIdx];
}
}
return 0.0;
}
Scalar zFraction_() const
{
if constexpr (enableExtbo) {
return (*this)[Indices::zFractionIdx];
}
else {
return 0.0;
}
}
Scalar saltConcentration_() const
{
if constexpr (enableBrine) {
return (*this)[Indices::saltConcentrationIdx];
}
else {
return 0.0;
}
}
Scalar biofilmVolumeFraction_() const
{
if constexpr (enableBioeffects)
return (*this)[Indices::biofilmVolumeFractionIdx];
else
return 0.0;
}
Scalar temperature_(const Problem& problem, [[maybe_unused]] unsigned globalDofIdx) const
{
if constexpr (enableEnergy) {
return (*this)[Indices::temperatureIdx];
}
else if constexpr (enableTemperature) {
return problem.temperature(globalDofIdx, /*timeIdx*/ 0);
}
else {
return FluidSystem::reservoirTemperature();
}
}
template <class Container>
void computeCapillaryPressures_(Container& result,
Scalar so,
Scalar sg,
Scalar sw,
const MaterialLawParams& matParams) const
{
using SatOnlyFluidState = SimpleModularFluidState<Scalar,
numPhases,
numComponents,
FluidSystem,
/*storePressure=*/false,
/*storeTemperature=*/false,
/*storeComposition=*/false,
/*storeFugacity=*/false,
/*storeSaturation=*/true,
/*storeDensity=*/false,
/*storeViscosity=*/false,
/*storeEnthalpy=*/false>;
SatOnlyFluidState fluidState;
fluidState.setSaturation(waterPhaseIdx, sw);
fluidState.setSaturation(oilPhaseIdx, so);
fluidState.setSaturation(gasPhaseIdx, sg);
MaterialLaw::capillaryPressures(result, matParams, fluidState);
}
Scalar pressure_() const
{ return (*this)[Indices::pressureSwitchIdx] * this->pressureScale_; }
void setScaledPressure_(Scalar pressure)
{ (*this)[Indices::pressureSwitchIdx] = pressure / (this->pressureScale_); }
WaterMeaning primaryVarsMeaningWater_{WaterMeaning::Disabled};
PressureMeaning primaryVarsMeaningPressure_{PressureMeaning::Po};
GasMeaning primaryVarsMeaningGas_{GasMeaning::Disabled};
BrineMeaning primaryVarsMeaningBrine_{BrineMeaning::Disabled};
SolventMeaning primaryVarsMeaningSolvent_{SolventMeaning::Disabled};
unsigned short pvtRegionIdx_;
Scalar pcFactor_;
static inline Scalar pressureScale_ = 1.0;
};
} // namespace Opm
#endif
blackoilbioeffectsmodules.hh
Contains the classes required to extend the black-oil model by bioeffects.
blackoilbrinemodules.hh
Contains the classes required to extend the black-oil model by brine.
blackoilenergymodules.hh
Contains the classes required to extend the black-oil model by energy.
blackoilextbomodules.hh
Contains the classes required to extend the black-oil model by solvent component.
blackoilproperties.hh
Declares the properties required by the black oil model.
blackoilsolventmodules.hh
Contains the classes required to extend the black-oil model by solvents.
Opm::BlackOilBioeffectsModule
Contains the high level supplements required to extend the black oil model by bioeffects.
Definition blackoilbioeffectsmodules.hh:93
Opm::BlackOilBrineModule
Contains the high level supplements required to extend the black oil model by brine.
Definition blackoilbrinemodules.hh:56
Opm::BlackOilEnergyModule
Contains the high level supplements required to extend the black oil model by energy.
Definition blackoilenergymodules.hh:58
Opm::BlackOilEnergyModule< TypeTag, enableEnergy >::assignPrimaryVars
static void assignPrimaryVars(PrimaryVariables &priVars, const FluidState &fluidState)
Definition blackoilenergymodules.hh:263
Opm::BlackOilExtboModule
Contains the high level supplements required to extend the black oil model.
Definition blackoilextbomodules.hh:62
Opm::BlackOilPrimaryVariables::pvtRegionIndex
unsigned pvtRegionIndex() const
Return the index of the region which should be used for PVT properties.
Definition blackoilprimaryvariables.hh:236
Opm::BlackOilPrimaryVariables::primaryVarsMeaningGas
GasMeaning primaryVarsMeaningGas() const
Return the interpretation which should be applied to the switching primary variables.
Definition blackoilprimaryvariables.hh:271
Opm::BlackOilPrimaryVariables::primaryVarsMeaningPressure
PressureMeaning primaryVarsMeaningPressure() const
Return the interpretation which should be applied to the switching primary variables.
Definition blackoilprimaryvariables.hh:257
Opm::BlackOilPrimaryVariables::setPrimaryVarsMeaningSolvent
void setPrimaryVarsMeaningSolvent(SolventMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition blackoilprimaryvariables.hh:298
Opm::BlackOilPrimaryVariables::adaptPrimaryVariables
bool adaptPrimaryVariables(const Problem &problem, unsigned globalDofIdx, Scalar swMaximum, Scalar thresholdWaterFilledCell, Scalar eps=0.0)
Adapt the interpretation of the switching variables to be physically meaningful.
Definition blackoilprimaryvariables.hh:476
Opm::BlackOilPrimaryVariables::setPrimaryVarsMeaningPressure
void setPrimaryVarsMeaningPressure(PressureMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition blackoilprimaryvariables.hh:264
Opm::BlackOilPrimaryVariables::primaryVarsMeaningWater
WaterMeaning primaryVarsMeaningWater() const
Return the interpretation which should be applied to the switching primary variables.
Definition blackoilprimaryvariables.hh:243
Opm::BlackOilPrimaryVariables::setPrimaryVarsMeaningWater
void setPrimaryVarsMeaningWater(WaterMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition blackoilprimaryvariables.hh:250
Opm::BlackOilPrimaryVariables::assignNaive
void assignNaive(const FluidState &fluidState)
Directly retrieve the primary variables from an arbitrary fluid state.
Definition blackoilprimaryvariables.hh:305
Opm::BlackOilPrimaryVariables::setPrimaryVarsMeaningBrine
void setPrimaryVarsMeaningBrine(BrineMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition blackoilprimaryvariables.hh:288
Opm::BlackOilPrimaryVariables::BlackOilPrimaryVariables
BlackOilPrimaryVariables(const BlackOilPrimaryVariables &value)=default
Copy constructor.
Opm::BlackOilPrimaryVariables::setPrimaryVarsMeaningGas
void setPrimaryVarsMeaningGas(GasMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition blackoilprimaryvariables.hh:278
Opm::BlackOilPrimaryVariables::setPvtRegionIndex
void setPvtRegionIndex(unsigned value)
Set the index of the region which should be used for PVT properties.
Definition blackoilprimaryvariables.hh:230
Opm::BlackOilPrimaryVariables::checkDefined
void checkDefined() const
< Import base class assignment operators.
Definition blackoilprimaryvariables.hh:909
Opm::BlackOilSolventModule
Contains the high level supplements required to extend the black oil model by solvents.
Definition blackoilsolventmodules.hh:68
Opm::FvBasePrimaryVariables::init
static void init()
< Import base class assignment operators.
Definition fvbaseprimaryvariables.hh:79
fvbaseprimaryvariables.hh
Represents the primary variables used by the a model.
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
Opm::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:233
Opm::getPropValue
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:240
Opm::Parameters::Get
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition parametersystem.hpp:187
Opm::Parameters::PressureScale
Definition blackoilprimaryvariables.hh:55