OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-07-04 11:33:06 -05:00
Code Issues Packages Projects Releases Wiki Activity

FlowExp: spring cleaning

Browse Source
This commit is contained in:
Arne Morten Kvarving 2024-04-25 08:50:16 +02:00
parent 1dfdae3892
commit 310625a3fe
5 changed files with 221 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
flowexperimental/BlackOilEnergyIntensiveQuantitiesGlobalIndex.hpp
View File

@ -27,16 +27,19 @@
*/ */
#ifndef OPM_BLACK_OIL_ENERGY_MODULE_GLOBAL_INDEX_HH #ifndef OPM_BLACK_OIL_ENERGY_MODULE_GLOBAL_INDEX_HH
#define OPM_BLACK_OIL_ENERGY_MODULE_GLOBAL_INDEX_HH #define OPM_BLACK_OIL_ENERGY_MODULE_GLOBAL_INDEX_HH
#include <opm/models/blackoil/blackoilenergymodules.hh> #include <opm/models/blackoil/blackoilenergymodules.hh>
namespace Opm { namespace Opm {
/*! /*!
* \ingroup BlackOil * \ingroup BlackOil
* \brief Contains the high level supplements required to extend the black oil * \brief Contains the high level supplements required to extend the black oil
* model by energy using global indices. * model by energy using global indices.
*/ */
template <class TypeTag, bool enableEnergyV = getPropValue<TypeTag, Properties::EnableEnergy>()> template <class TypeTag, bool enableEnergyV = getPropValue<TypeTag, Properties::EnableEnergy>()>
class BlackOilEnergyIntensiveQuantitiesGlobalIndex: public BlackOilEnergyIntensiveQuantities<TypeTag,enableEnergyV> class BlackOilEnergyIntensiveQuantitiesGlobalIndex
: public BlackOilEnergyIntensiveQuantities<TypeTag,enableEnergyV>
{ {
using Parent = BlackOilEnergyIntensiveQuantities<TypeTag, enableEnergyV>; using Parent = BlackOilEnergyIntensiveQuantities<TypeTag, enableEnergyV>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
@ -46,15 +49,17 @@ class BlackOilEnergyIntensiveQuantitiesGlobalIndex: public BlackOilEnergyIntensi
using Scalar = GetPropType<TypeTag, Properties::Scalar>; using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>; using SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>;
using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>; using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
using ParamCache = typename FluidSystem::template ParameterCache<Evaluation>;
static constexpr bool enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>(); static constexpr bool enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>();
using Indices = GetPropType<TypeTag, Properties::Indices>; using Indices = GetPropType<TypeTag, Properties::Indices>;
static constexpr unsigned temperatureIdx = Indices::temperatureIdx; static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
static constexpr unsigned numPhases = FluidSystem::numPhases; static constexpr unsigned numPhases = FluidSystem::numPhases;
public: public:
void updateTemperature_(const Problem& problem, void updateTemperature_([[maybe_unused]] const Problem& problem,
const PrimaryVariables& priVars, const PrimaryVariables& priVars,
unsigned globalSpaceIndex, [[maybe_unused]] unsigned globalSpaceIndex,
unsigned timeIdx) unsigned timeIdx)
{ {
auto& fs = Parent::asImp_().fluidState_; auto& fs = Parent::asImp_().fluidState_;
@ -62,10 +67,10 @@ public:
fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx)); fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx));
} }
void updateEnergyQuantities_(const Problem& problem, void updateEnergyQuantities_(const Problem& problem,
const PrimaryVariables& priVars, [[maybe_unused]] const PrimaryVariables& priVars,
unsigned globalSpaceIndex, unsigned globalSpaceIndex,
unsigned timeIdx, unsigned timeIdx,
const typename FluidSystem::template ParameterCache<Evaluation>& paramCache) const ParamCache& paramCache)
{ {
auto& fs = Parent::asImp_().fluidState_; auto& fs = Parent::asImp_().fluidState_;
@ -93,42 +98,43 @@ public:
// multiplier. This is to avoid negative rock volume for pvmult*porosity > 1 // multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
this->rockFraction_ = problem.rockFraction(globalSpaceIndex, timeIdx); this->rockFraction_ = problem.rockFraction(globalSpaceIndex, timeIdx);
} }
}; };
template <class TypeTag> template <class TypeTag>
class BlackOilEnergyIntensiveQuantitiesGlobalIndex<TypeTag, false>: public BlackOilEnergyIntensiveQuantities<TypeTag, false> class BlackOilEnergyIntensiveQuantitiesGlobalIndex<TypeTag, false>
: public BlackOilEnergyIntensiveQuantities<TypeTag, false>
{
using Parent = BlackOilEnergyIntensiveQuantities<TypeTag, false>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
static constexpr bool enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>();
public:
void updateTemperature_([[maybe_unused]] const Problem& problem,
[[maybe_unused]] const PrimaryVariables& priVars,
[[maybe_unused]] unsigned globalSpaceIdx,
[[maybe_unused]] unsigned timeIdx)
{ {
using Parent = BlackOilEnergyIntensiveQuantities<TypeTag, false>; if constexpr (enableTemperature) {
using Problem = GetPropType<TypeTag, Properties::Problem>; // even if energy is conserved, the temperature can vary over the spatial
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; // domain if the EnableTemperature property is set to true
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; auto& fs = this->asImp_().fluidState_;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>; Scalar T = problem.temperature(globalSpaceIdx, timeIdx);
using Scalar = GetPropType<TypeTag, Properties::Scalar>; fs.setTemperature(T);
static constexpr bool enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>();
public:
void updateTemperature_([[maybe_unused]] const Problem& problem,
[[maybe_unused]] const PrimaryVariables& priVars,
[[maybe_unused]] unsigned globalSpaceIdx,
[[maybe_unused]] unsigned timeIdx)
{
if constexpr (enableTemperature) {
// even if energy is conserved, the temperature can vary over the spatial
// domain if the EnableTemperature property is set to true
auto& fs = this->asImp_().fluidState_;
Scalar T = problem.temperature(globalSpaceIdx, timeIdx);
fs.setTemperature(T);
}
} }
void updateEnergyQuantities_([[maybe_unused]] const Problem& problem, }
[[maybe_unused]] const PrimaryVariables& priVars,
[[maybe_unused]] unsigned globalSpaceIdx,
[[maybe_unused]] unsigned timeIdx,
const typename FluidSystem::template ParameterCache<Evaluation>&)
{ }
void updateEnergyQuantities_([[maybe_unused]] const Problem& problem,
[[maybe_unused]] const PrimaryVariables& priVars,
[[maybe_unused]] unsigned globalSpaceIdx,
[[maybe_unused]] unsigned timeIdx,
const typename FluidSystem::template ParameterCache<Evaluation>&)
{ }
}; };
} } // namespace Opm
#endif #endif

142
flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp
View File

@ -27,6 +27,19 @@
*/ */
#ifndef OPM_BLACK_OIL_INTENSIVE_QUANTITIES_GLOBAL_INDEX_HPP #ifndef OPM_BLACK_OIL_INTENSIVE_QUANTITIES_GLOBAL_INDEX_HPP
#define OPM_BLACK_OIL_INTENSIVE_QUANTITIES_GLOBAL_INDEX_HPP #define OPM_BLACK_OIL_INTENSIVE_QUANTITIES_GLOBAL_INDEX_HPP
#include "BlackOilEnergyIntensiveQuantitiesGlobalIndex.hpp"
#include <dune/common/fmatrix.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/models/blackoil/blackoilproperties.hh> #include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/blackoil/blackoilsolventmodules.hh> #include <opm/models/blackoil/blackoilsolventmodules.hh>
#include <opm/models/blackoil/blackoilextbomodules.hh> #include <opm/models/blackoil/blackoilextbomodules.hh>
@ -34,24 +47,18 @@
#include <opm/models/blackoil/blackoilfoammodules.hh> #include <opm/models/blackoil/blackoilfoammodules.hh>
#include <opm/models/blackoil/blackoilbrinemodules.hh> #include <opm/models/blackoil/blackoilbrinemodules.hh>
#include <opm/models/blackoil/blackoilenergymodules.hh> #include <opm/models/blackoil/blackoilenergymodules.hh>
#include "BlackOilEnergyIntensiveQuantitiesGlobalIndex.hpp"
#include <opm/models/blackoil/blackoildiffusionmodule.hh> #include <opm/models/blackoil/blackoildiffusionmodule.hh>
#include <opm/models/blackoil/blackoilmicpmodules.hh> #include <opm/models/blackoil/blackoilmicpmodules.hh>
#include <opm/material/fluidstates/BlackOilFluidState.hpp> #include <opm/models/common/directionalmobility.hh>
#include <opm/material/common/Valgrind.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/utility/CopyablePtr.hpp> #include <opm/utility/CopyablePtr.hpp>
#include <dune/common/fmatrix.hh>
#include <cstring>
#include <utility> #include <utility>
#include <fmt/format.h> #include <fmt/format.h>
namespace Opm { namespace Opm {
/*! /*!
* \ingroup BlackOilModel * \ingroup BlackOilModel
* \ingroup IntensiveQuantities * \ingroup IntensiveQuantities
@ -120,7 +127,6 @@ class BlackOilIntensiveQuantitiesGlobalIndex
using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>; using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
public: public:
using FluidState = BlackOilFluidState<Evaluation, using FluidState = BlackOilFluidState<Evaluation,
FluidSystem, FluidSystem,
@ -145,10 +151,6 @@ public:
} }
} }
BlackOilIntensiveQuantitiesGlobalIndex(const BlackOilIntensiveQuantitiesGlobalIndex& other) = default;
BlackOilIntensiveQuantitiesGlobalIndex& operator=(const BlackOilIntensiveQuantitiesGlobalIndex& other) = default;
/*! /*!
* \copydoc IntensiveQuantities::update * \copydoc IntensiveQuantities::update
*/ */
@ -189,7 +191,7 @@ public:
if constexpr (waterEnabled) { if constexpr (waterEnabled) {
if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) { if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
Sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx); Sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
} else if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Disabled){ } else if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Disabled) {
// water is enabled but is not a primary variable i.e. one phase case // water is enabled but is not a primary variable i.e. one phase case
Sw = 1.0; Sw = 1.0;
} }
@ -215,31 +217,39 @@ public:
Evaluation So = 1.0 - Sw - Sg; Evaluation So = 1.0 - Sw - Sg;
// deal with solvent // deal with solvent
if constexpr (enableSolvent) if constexpr (enableSolvent) {
So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx); So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
}
if (FluidSystem::phaseIsActive(waterPhaseIdx)) if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
fluidState_.setSaturation(waterPhaseIdx, Sw); fluidState_.setSaturation(waterPhaseIdx, Sw);
}
if (FluidSystem::phaseIsActive(gasPhaseIdx)) if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
fluidState_.setSaturation(gasPhaseIdx, Sg); fluidState_.setSaturation(gasPhaseIdx, Sg);
}
if (FluidSystem::phaseIsActive(oilPhaseIdx)) if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
fluidState_.setSaturation(oilPhaseIdx, So); fluidState_.setSaturation(oilPhaseIdx, So);
}
std::array<Evaluation, numPhases> pC;// = {0, 0, 0}; std::array<Evaluation, numPhases> pC{};
computeRelpermAndPC(mobility_, pC, problem, fluidState_, globalSpaceIdx); computeRelpermAndPC(mobility_, pC, problem, fluidState_, globalSpaceIdx);
// oil is the reference phase for pressure // oil is the reference phase for pressure
if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) { if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx); const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (FluidSystem::phaseIsActive(phaseIdx)) if (FluidSystem::phaseIsActive(phaseIdx)) {
fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx])); fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
}
}
} else { } else {
const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx); const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (FluidSystem::phaseIsActive(phaseIdx)) if (FluidSystem::phaseIsActive(phaseIdx)) {
fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx])); fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
}
}
} }
Evaluation SoMax = 0.0; Evaluation SoMax = 0.0;
@ -263,8 +273,9 @@ public:
SoMax); SoMax);
fluidState_.setRs(min(RsMax, RsSat)); fluidState_.setRs(min(RsMax, RsSat));
} }
else if constexpr (compositionSwitchEnabled) else if constexpr (compositionSwitchEnabled) {
fluidState_.setRs(0.0); fluidState_.setRs(0.0);
}
} }
if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) { if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
@ -281,8 +292,9 @@ public:
SoMax); SoMax);
fluidState_.setRv(min(RvMax, RvSat)); fluidState_.setRv(min(RvMax, RvSat));
} }
else if constexpr (compositionSwitchEnabled) else if constexpr (compositionSwitchEnabled) {
fluidState_.setRv(0.0); fluidState_.setRv(0.0);
}
} }
if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rvw) { if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rvw) {
@ -300,8 +312,9 @@ public:
} }
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue; continue;
}
computeInverseFormationVolumeFactorAndViscosity(fluidState_, phaseIdx, pvtRegionIdx, SoMax); computeInverseFormationVolumeFactorAndViscosity(fluidState_, phaseIdx, pvtRegionIdx, SoMax);
} }
Valgrind::CheckDefined(mobility_); Valgrind::CheckDefined(mobility_);
@ -320,16 +333,14 @@ public:
rho = fluidState_.invB(gasPhaseIdx); rho = fluidState_.invB(gasPhaseIdx);
rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx); rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableVaporizedOil()) { if (FluidSystem::enableVaporizedOil()) {
rho += rho += fluidState_.invB(gasPhaseIdx) *
fluidState_.invB(gasPhaseIdx) * fluidState_.Rv() *
fluidState_.Rv() * FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
} }
if (FluidSystem::enableVaporizedWater()) { if (FluidSystem::enableVaporizedWater()) {
rho += rho += fluidState_.invB(gasPhaseIdx) *
fluidState_.invB(gasPhaseIdx) * fluidState_.Rvw() *
fluidState_.Rvw() * FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
} }
fluidState_.setDensity(gasPhaseIdx, rho); fluidState_.setDensity(gasPhaseIdx, rho);
} }
@ -339,10 +350,9 @@ public:
rho = fluidState_.invB(oilPhaseIdx); rho = fluidState_.invB(oilPhaseIdx);
rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx); rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableDissolvedGas()) { if (FluidSystem::enableDissolvedGas()) {
rho += rho += fluidState_.invB(oilPhaseIdx) *
fluidState_.invB(oilPhaseIdx) * fluidState_.Rs() *
fluidState_.Rs() * FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
} }
fluidState_.setDensity(oilPhaseIdx, rho); fluidState_.setDensity(oilPhaseIdx, rho);
} }
@ -359,12 +369,12 @@ public:
Evaluation x; Evaluation x;
if (FluidSystem::phaseIsActive(oilPhaseIdx)) { if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure); x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
} else if (FluidSystem::phaseIsActive(waterPhaseIdx)){ } else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure); x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
} else { } else {
x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure); x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
} }
porosity_ *= 1.0 + x + 0.5*x*x; porosity_ *= 1.0 + x + 0.5 * x * x;
} }
// deal with water induced rock compaction // deal with water induced rock compaction
@ -375,8 +385,9 @@ public:
#ifndef NDEBUG #ifndef NDEBUG
// some safety checks in debug mode // some safety checks in debug mode
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue; continue;
}
assert(isfinite(fluidState_.density(phaseIdx))); assert(isfinite(fluidState_.density(phaseIdx)));
assert(isfinite(fluidState_.saturation(phaseIdx))); assert(isfinite(fluidState_.saturation(phaseIdx)));
@ -389,7 +400,6 @@ public:
#endif #endif
} }
/*! /*!
* \copydoc ImmiscibleIntensiveQuantities::fluidState * \copydoc ImmiscibleIntensiveQuantities::fluidState
*/ */
@ -407,20 +417,18 @@ public:
using Dir = FaceDir::DirEnum; using Dir = FaceDir::DirEnum;
if (dirMob_) { if (dirMob_) {
switch(facedir) { switch(facedir) {
case Dir::XPlus: case Dir::XPlus:
return dirMob_->mobilityX_[phaseIdx]; return dirMob_->mobilityX_[phaseIdx];
case Dir::YPlus: case Dir::YPlus:
return dirMob_->mobilityY_[phaseIdx]; return dirMob_->mobilityY_[phaseIdx];
case Dir::ZPlus: case Dir::ZPlus:
return dirMob_->mobilityZ_[phaseIdx]; return dirMob_->mobilityZ_[phaseIdx];
default: default:
throw std::runtime_error("Unexpected face direction"); throw std::runtime_error("Unexpected face direction");
} }
} } else {
else {
return mobility_[phaseIdx]; return mobility_[phaseIdx];
} }
} }
void computeInverseFormationVolumeFactorAndViscosity(FluidState& fluidState, void computeInverseFormationVolumeFactorAndViscosity(FluidState& fluidState,
@ -429,21 +437,21 @@ public:
const Evaluation& SoMax){ const Evaluation& SoMax){
OPM_TIMEBLOCK_LOCAL(UpdateInverseFormationFactorAndViscosity); OPM_TIMEBLOCK_LOCAL(UpdateInverseFormationFactorAndViscosity);
{ {
OPM_TIMEBLOCK_LOCAL(UpdateFormationFactor); OPM_TIMEBLOCK_LOCAL(UpdateFormationFactor);
const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx); const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
fluidState_.setInvB(phaseIdx, b); fluidState_.setInvB(phaseIdx, b);
} }
{ {
OPM_TIMEBLOCK_LOCAL(UpdateViscosity); OPM_TIMEBLOCK_LOCAL(UpdateViscosity);
typename FluidSystem::template ParameterCache<Evaluation> paramCache; typename FluidSystem::template ParameterCache<Evaluation> paramCache;
paramCache.setRegionIndex(pvtRegionIdx); paramCache.setRegionIndex(pvtRegionIdx);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
paramCache.setMaxOilSat(SoMax); paramCache.setMaxOilSat(SoMax);
} }
paramCache.updateAll(fluidState_); paramCache.updateAll(fluidState_);
const auto& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx); const auto& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
mobility_[phaseIdx] /= mu; mobility_[phaseIdx] /= mu;
} }
} }
@ -457,6 +465,7 @@ public:
MaterialLaw::capillaryPressures(pC, materialParams, fluidState_); MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
problem.updateRelperms(mobility, dirMob_, fluidState, globalSpaceIdx); problem.updateRelperms(mobility, dirMob_, fluidState, globalSpaceIdx);
} }
/*! /*!
* \copydoc ImmiscibleIntensiveQuantities::porosity * \copydoc ImmiscibleIntensiveQuantities::porosity
*/ */
@ -482,8 +491,7 @@ public:
* support it as well in our black-oil model. (Note that, if it is not explicitly * support it as well in our black-oil model. (Note that, if it is not explicitly
* specified, the PVT region index is 0.) * specified, the PVT region index is 0.)
*/ */
auto pvtRegionIndex() const auto pvtRegionIndex() const -> decltype(std::declval<FluidState>().pvtRegionIndex())
-> decltype(std::declval<FluidState>().pvtRegionIndex())
{ return fluidState_.pvtRegionIndex(); } { return fluidState_.pvtRegionIndex(); }
/*! /*!

54
flowexperimental/FIBlackOilModelNoCache.hpp
View File

@ -1,32 +1,38 @@
#ifndef OPM_FI_BLACK_OIL_MODEL_NOCACHE_HPP #ifndef OPM_FI_BLACK_OIL_MODEL_NOCACHE_HPP
#define OPM_FI_BLACK_OIL_MODEL_NOCACHE_HPP #define OPM_FI_BLACK_OIL_MODEL_NOCACHE_HPP
#include <opm/simulators/flow/FIBlackoilModel.hpp> #include <opm/simulators/flow/FIBlackoilModel.hpp>
namespace Opm{
template<typename TypeTag>
class FIBlackOilModelNoCache: public FIBlackOilModel<TypeTag>{
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
public:
explicit FIBlackOilModelNoCache(Simulator& simulator)
:FIBlackOilModel<TypeTag>(simulator)
{ }
IntensiveQuantities intensiveQuantities(unsigned globalIdx, unsigned timeIdx) const{ namespace Opm {
OPM_TIMEBLOCK_LOCAL(intensiveQuantitiesNoCache);
const auto& primaryVar = this->solution(timeIdx)[globalIdx];
const auto& problem = this->simulator_.problem();
if (!(this->enableIntensiveQuantityCache_) ||
!(this->intensiveQuantityCacheUpToDate_[timeIdx][globalIdx])){
IntensiveQuantities intQuants;
intQuants.update(problem,primaryVar, globalIdx, timeIdx);
return intQuants;// reqiored for updating extrution factor
}else{
IntensiveQuantities intQuants = (this->intensiveQuantityCache_[timeIdx][globalIdx]);
return intQuants;
}
template<typename TypeTag>
class FIBlackOilModelNoCache: public FIBlackOilModel<TypeTag>
{
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
public:
explicit FIBlackOilModelNoCache(Simulator& simulator)
: FIBlackOilModel<TypeTag>(simulator)
{}
IntensiveQuantities intensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
{
OPM_TIMEBLOCK_LOCAL(intensiveQuantitiesNoCache);
const auto& primaryVar = this->solution(timeIdx)[globalIdx];
const auto& problem = this->simulator_.problem();
if (!(this->enableIntensiveQuantityCache_) ||
!(this->intensiveQuantityCacheUpToDate_[timeIdx][globalIdx])) {
IntensiveQuantities intQuants;
intQuants.update(problem,primaryVar, globalIdx, timeIdx);
return intQuants;// reqiored for updating extrution factor
} else {
IntensiveQuantities intQuants = (this->intensiveQuantityCache_[timeIdx][globalIdx]);
return intQuants;
} }
}
};
} // namespace Opm
};
}
#endif #endif

59
flowexperimental/FlowExpNewtonMethod.hpp
View File

@ -37,23 +37,32 @@
namespace Opm::Properties { namespace Opm::Properties {
template<class TypeTag, class MyTypeTag> template<class TypeTag, class MyTypeTag>
struct EclNewtonSumTolerance { struct EclNewtonSumTolerance
{
using type = UndefinedProperty; using type = UndefinedProperty;
}; };
template<class TypeTag, class MyTypeTag> template<class TypeTag, class MyTypeTag>
struct EclNewtonStrictIterations { struct EclNewtonStrictIterations
{
using type = UndefinedProperty; using type = UndefinedProperty;
}; };
template<class TypeTag, class MyTypeTag> template<class TypeTag, class MyTypeTag>
struct EclNewtonRelaxedVolumeFraction { struct EclNewtonRelaxedVolumeFraction
{
using type = UndefinedProperty; using type = UndefinedProperty;
}; };
template<class TypeTag, class MyTypeTag> template<class TypeTag, class MyTypeTag>
struct EclNewtonSumToleranceExponent { struct EclNewtonSumToleranceExponent
{
using type = UndefinedProperty; using type = UndefinedProperty;
}; };
template<class TypeTag, class MyTypeTag> template<class TypeTag, class MyTypeTag>
struct EclNewtonRelaxedTolerance { struct EclNewtonRelaxedTolerance
{
using type = UndefinedProperty; using type = UndefinedProperty;
}; };
@ -133,10 +142,11 @@ public:
*/ */
bool converged() const bool converged() const
{ {
if (errorPvFraction_ < relaxedMaxPvFraction_) if (errorPvFraction_ < relaxedMaxPvFraction_) {
return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ;
else if (this->numIterations() > numStrictIterations_) } else if (this->numIterations() > numStrictIterations_) {
return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ;
}
return this->error_ <= this->tolerance() && errorSum_ <= sumTolerance_; return this->error_ <= this->tolerance() && errorSum_ <= sumTolerance_;
} }
@ -159,22 +169,24 @@ public:
for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) { for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
// do not consider auxiliary DOFs for the error // do not consider auxiliary DOFs for the error
if (dofIdx >= this->model().numGridDof() if (dofIdx >= this->model().numGridDof()
|| this->model().dofTotalVolume(dofIdx) <= 0.0) || this->model().dofTotalVolume(dofIdx) <= 0.0) {
continue; continue;
}
if (!this->model().isLocalDof(dofIdx)) if (!this->model().isLocalDof(dofIdx)) {
continue; continue;
}
// also do not consider DOFs which are constraint // also do not consider DOFs which are constraint
if (this->enableConstraints_()) { if (this->enableConstraints_()) {
if (constraintsMap.count(dofIdx) > 0) if (constraintsMap.count(dofIdx) > 0) {
continue; continue;
}
} }
const auto& r = currentResidual[dofIdx]; const auto& r = currentResidual[dofIdx];
Scalar pvValue = Scalar pvValue = this->simulator_.problem().referencePorosity(dofIdx, /*timeIdx=*/0) *
this->simulator_.problem().referencePorosity(dofIdx, /*timeIdx=*/0) this->model().dofTotalVolume(dofIdx);
* this->model().dofTotalVolume(dofIdx);
sumPv += pvValue; sumPv += pvValue;
bool cnvViolated = false; bool cnvViolated = false;
@ -193,13 +205,15 @@ public:
this->error_ = max(std::abs(tmpError), this->error_); this->error_ = max(std::abs(tmpError), this->error_);
if (std::abs(tmpError) > this->tolerance_) if (std::abs(tmpError) > this->tolerance_) {
cnvViolated = true; cnvViolated = true;
}
componentSumError[eqIdx] += std::abs(tmpError2); componentSumError[eqIdx] += std::abs(tmpError2);
} }
if (cnvViolated) if (cnvViolated) {
errorPvFraction_ += pvValue; errorPvFraction_ += pvValue;
}
} }
// take the other processes into account // take the other processes into account
@ -214,8 +228,9 @@ public:
errorPvFraction_ /= sumPv; errorPvFraction_ /= sumPv;
errorSum_ = 0; errorSum_ = 0;
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
errorSum_ = std::max(std::abs(componentSumError[eqIdx]), errorSum_); errorSum_ = std::max(std::abs(componentSumError[eqIdx]), errorSum_);
}
// scale the tolerance for the total error with the pore volume. by default, the // scale the tolerance for the total error with the pore volume. by default, the
// exponent is 1/3, i.e., cubic root. // exponent is 1/3, i.e., cubic root.
@ -223,21 +238,26 @@ public:
Scalar y = Parameters::get<TypeTag, Properties::EclNewtonSumToleranceExponent>(); Scalar y = Parameters::get<TypeTag, Properties::EclNewtonSumToleranceExponent>();
sumTolerance_ = x*std::pow(sumPv, y); sumTolerance_ = x*std::pow(sumPv, y);
this->endIterMsg() << " (max: " << this->tolerance_ << ", violated for " << errorPvFraction_*100 << "% of the pore volume), aggegate error: " << errorSum_ << " (max: " << sumTolerance_ << ")"; this->endIterMsg() << " (max: " << this->tolerance_
<< ", violated for " << errorPvFraction_ * 100
<< "% of the pore volume), aggegate error: "
<< errorSum_ << " (max: " << sumTolerance_ << ")";
// make sure that the error never grows beyond the maximum // make sure that the error never grows beyond the maximum
// allowed one // allowed one
if (this->error_ > newtonMaxError) if (this->error_ > newtonMaxError) {
throw NumericalProblem("Newton: Error "+std::to_string(double(this->error_)) throw NumericalProblem("Newton: Error "+std::to_string(double(this->error_))
+ " is larger than maximum allowed error of " + " is larger than maximum allowed error of "
+ std::to_string(double(newtonMaxError))); + std::to_string(double(newtonMaxError)));
}
// make sure that the error never grows beyond the maximum // make sure that the error never grows beyond the maximum
// allowed one // allowed one
if (errorSum_ > newtonMaxError) if (errorSum_ > newtonMaxError) {
throw NumericalProblem("Newton: Sum of the error "+std::to_string(double(errorSum_)) throw NumericalProblem("Newton: Sum of the error "+std::to_string(double(errorSum_))
+ " is larger than maximum allowed error of " + " is larger than maximum allowed error of "
+ std::to_string(double(newtonMaxError))); + std::to_string(double(newtonMaxError)));
}
} }
void endIteration_(SolutionVector& nextSolution, void endIteration_(SolutionVector& nextSolution,
@ -261,6 +281,7 @@ private:
int numStrictIterations_; int numStrictIterations_;
}; };
} // namespace Opm } // namespace Opm
#endif #endif

52
flowexperimental/flowexp_blackoil.cpp
View File

@ -29,45 +29,54 @@
namespace Opm::Properties { namespace Opm::Properties {
namespace TTag { namespace TTag {
struct FlowExpProblemBlackOil{ struct FlowExpProblemBlackOil
{
using InheritsFrom = std::tuple<FlowExpTypeTag>; using InheritsFrom = std::tuple<FlowExpTypeTag>;
}; };
} }
template<class TypeTag> template<class TypeTag>
struct Model<TypeTag, TTag::FlowExpProblemBlackOil> { struct Model<TypeTag, TTag::FlowExpProblemBlackOil>
{
using type = FIBlackOilModelNoCache<TypeTag>; using type = FIBlackOilModelNoCache<TypeTag>;
}; };
template<class TypeTag> template<class TypeTag>
struct IntensiveQuantities<TypeTag, TTag::FlowExpProblemBlackOil> { struct IntensiveQuantities<TypeTag, TTag::FlowExpProblemBlackOil>
using type = BlackOilIntensiveQuantitiesGlobalIndex<TypeTag>; {
using type = BlackOilIntensiveQuantitiesGlobalIndex<TypeTag>;
}; };
// Set the problem class // Set the problem class
template<class TypeTag> template<class TypeTag>
struct Problem<TypeTag, TTag::FlowExpProblemBlackOil> { struct Problem<TypeTag, TTag::FlowExpProblemBlackOil>
{
using type = FlowExpProblem<TypeTag>; using type = FlowExpProblem<TypeTag>;
}; };
template<class TypeTag> template<class TypeTag>
struct ThreadsPerProcess<TypeTag, TTag::FlowExpProblemBlackOil> { struct ThreadsPerProcess<TypeTag, TTag::FlowExpProblemBlackOil>
{
static constexpr int value = 1; static constexpr int value = 1;
}; };
template<class TypeTag> template<class TypeTag>
struct ContinueOnConvergenceError<TypeTag, TTag::FlowExpProblemBlackOil> { struct ContinueOnConvergenceError<TypeTag, TTag::FlowExpProblemBlackOil>
{
static constexpr bool value = false; static constexpr bool value = false;
}; };
template<class TypeTag> template<class TypeTag>
struct EclNewtonSumTolerance<TypeTag, TTag::FlowExpProblemBlackOil> { struct EclNewtonSumTolerance<TypeTag, TTag::FlowExpProblemBlackOil>
{
using type = GetPropType<TypeTag, Scalar>; using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1e-5; static constexpr type value = 1e-5;
}; };
// the default for the allowed volumetric error for oil per second // the default for the allowed volumetric error for oil per second
template<class TypeTag> template<class TypeTag>
struct NewtonTolerance<TypeTag, TTag::FlowExpProblemBlackOil> { struct NewtonTolerance<TypeTag, TTag::FlowExpProblemBlackOil>
{
using type = GetPropType<TypeTag, Scalar>; using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 1e-2; static constexpr type value = 1e-2;
}; };
@ -75,30 +84,39 @@ struct NewtonTolerance<TypeTag, TTag::FlowExpProblemBlackOil> {
// set fraction of the pore volume where the volumetric residual may be violated during // set fraction of the pore volume where the volumetric residual may be violated during
// strict Newton iterations // strict Newton iterations
template<class TypeTag> template<class TypeTag>
struct EclNewtonRelaxedVolumeFraction<TypeTag, TTag::FlowExpProblemBlackOil> { struct EclNewtonRelaxedVolumeFraction<TypeTag, TTag::FlowExpProblemBlackOil>
{
using type = GetPropType<TypeTag, Scalar>; using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 0.0; static constexpr type value = 0.0;
}; };
template<class TypeTag> template<class TypeTag>
struct EclNewtonRelaxedTolerance<TypeTag, TTag::FlowExpProblemBlackOil> { struct EclNewtonRelaxedTolerance<TypeTag, TTag::FlowExpProblemBlackOil>
{
using type = GetPropType<TypeTag, Scalar>; using type = GetPropType<TypeTag, Scalar>;
static constexpr type value = 10*getPropValue<TypeTag, Properties::NewtonTolerance>(); static constexpr type value = 10*getPropValue<TypeTag, Properties::NewtonTolerance>();
}; };
template<class TypeTag> template<class TypeTag>
struct EnableDiffusion<TypeTag, TTag::FlowExpProblemBlackOil> { static constexpr bool value = false; }; struct EnableDiffusion<TypeTag, TTag::FlowExpProblemBlackOil>
{
static constexpr bool value = false;
};
template<class TypeTag> template<class TypeTag>
struct EnableDisgasInWater<TypeTag, TTag::FlowExpProblemBlackOil> { static constexpr bool value = false; }; struct EnableDisgasInWater<TypeTag, TTag::FlowExpProblemBlackOil>
{
static constexpr bool value = false;
};
template<class TypeTag> template<class TypeTag>
struct Simulator<TypeTag, TTag::FlowExpProblemBlackOil> { using type = Opm::Simulator<TypeTag>; }; struct Simulator<TypeTag, TTag::FlowExpProblemBlackOil>
{
using type = Opm::Simulator<TypeTag>;
};
} }
int main(int argc, char** argv) int main(int argc, char** argv)
{ {
using TypeTag = Opm::Properties::TTag::FlowExpProblemBlackOil; using TypeTag = Opm::Properties::TTag::FlowExpProblemBlackOil;