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

fix for flag when DRSDTCON is active, adressed warnings and fixed initialization

Browse Source
This commit is contained in:
Trine Mykkeltvedt 2024-06-24 12:41:53 +02:00 committed by Tor Harald Sandve
parent 3d943770a5
commit 9696ff0824
5 changed files with 57 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

86
opm/models/blackoil/blackoilconvectivemixingmodule.hh
View File

@ -43,13 +43,13 @@ namespace Opm {
/*!
* \copydoc Opm::BlackOilConvectiveMixingModule
* \brief Provides the convective term in the transport flux for the brine
* when convective mixing (enhanced dissolution of CO2 in brine) occurs.
* when convective mixing (enhanced dissolution of CO2 in brine) occurs.
* Controlled by the regimes for a controlvolume:
* i) initial phase (CO2 dissolves in brine due to diffusion)
* ii) linear phase (Convective fingers of CO2-rich brine propagate downwards)
* iii) steady-state-phase (fingers have passed through the bottom of a control
* -volume but the larger scale convective process is still active)
* iv) decline phase (Convection ceases at the large-scale when the CO2
* iv) decline phase (Convection ceases at the large-scale when the CO2
* has been completely dissolved)
*/
@ -84,25 +84,32 @@ public:
#endif
template <class Context>
static void addConvectiveMixingFlux(RateVector& flux,
const Context& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
static bool active(const Context&) {
return false;
}
template <class Context>
static void addConvectiveMixingFlux(RateVector&,
const Context&,
unsigned,
unsigned)
{}
/*!
* \brief Adds the convective mixing mass flux flux to the flux vector over a flux
* integration point.
*/
static void addConvectiveMixingFlux(RateVector& flux,
const IntensiveQuantities& intQuantsIn,
const IntensiveQuantities& intQuantsEx,
const unsigned globalIndexIn,
const unsigned globalIndexEx,
const Scalar distZg,
const Scalar trans,
const Scalar faceArea,
const ConvectiveMixingModuleParam& info)
static void addConvectiveMixingFlux(RateVector&,
const IntensiveQuantities&,
const IntensiveQuantities&,
const unsigned,
const unsigned,
const Scalar,
const Scalar,
const Scalar,
const ConvectiveMixingModuleParam&)
{}
};
template <class TypeTag>
class BlackOilConvectiveMixingModule<TypeTag, /*enableConvectiveMixing=*/true>
@ -122,7 +129,7 @@ public:
struct ConvectiveMixingModuleParam
{
bool active_;
bool active_{false};
std::vector<Scalar> Xhi_;
std::vector<Scalar> Psi_;
};
@ -138,8 +145,8 @@ public:
return;
}
if (info.Xhi_.empty()) {
info.Xhi_.resize(numRegions);
info.Psi_.resize(numRegions);
info.Xhi_.resize(numRegions);
info.Psi_.resize(numRegions);
}
for (size_t i = 0; i < numRegions; ++i ) {
info.Xhi_[i] = control.getMaxDRSDT(i);
@ -149,12 +156,17 @@ public:
#endif
template <class Context>
static void addConvectiveMixingFlux(RateVector& flux,
static bool active(const Context& elemCtx) {
const auto& problem = elemCtx.problem();
return problem.moduleParams().convectiveMixingModuleParam.active_;
}
template <class Context>
static void addConvectiveMixingFlux(RateVector& flux,
const Context& elemCtx,
unsigned scvfIdx,
unsigned scvfIdx,
unsigned timeIdx)
{
// need for dary flux calculation
const auto& problem = elemCtx.problem();
const auto& stencil = elemCtx.stencil(timeIdx);
@ -185,10 +197,6 @@ public:
problem.moduleParams().convectiveMixingModuleParam);
}
/*!
* \brief Adds the convective mixing mass flux flux to the flux vector over a flux
* integration point.
@ -198,19 +206,17 @@ public:
const IntensiveQuantities& intQuantsEx,
const unsigned globalIndexIn,
const unsigned globalIndexEx,
const Scalar distZg,
const Scalar distZg,
const Scalar trans,
const Scalar faceArea,
const ConvectiveMixingModuleParam& info)
{
if (!info.active_) {
return;
}
const auto& liquidPhaseIdx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPhaseIdx :
FluidSystem::waterPhaseIdx :
FluidSystem::oilPhaseIdx;
const Evaluation SoMax = 0.0;
@ -225,11 +231,11 @@ public:
const auto& salt_in = Opm::getValue(intQuantsIn.fluidState().saltSaturation());
const auto& bLiquidIn = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
const auto& bLiquidIn = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_in, p_in, rs_zero_in, salt_in):
FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_in, p_in, rs_zero_in);
const auto& bLiquidSatIn = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
const auto& bLiquidSatIn = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPvt().saturatedInverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_in, p_in, salt_in) :
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_in, p_in);
@ -240,7 +246,7 @@ public:
const auto rho_in = bLiquidIn * densityLiquidIn;
const auto rho_sat_in = bLiquidSatIn
* (densityLiquidIn + rssat_in * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, intQuantsIn.pvtRegionIndex()));
//exteriour
const Scalar rs_zero_ex = 0.0;
const auto& t_ex = Opm::getValue(intQuantsEx.fluidState().temperature(liquidPhaseIdx));
@ -250,24 +256,21 @@ public:
intQuantsEx.pvtRegionIndex(),
SoMax);
const auto& salt_ex = Opm::getValue(intQuantsEx.fluidState().saltSaturation());
const auto& bLiquidEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
const auto& bLiquidEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_ex, p_ex, rs_zero_ex, salt_ex) :
FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_ex, p_ex, rs_zero_ex);
const auto& bLiquidSatEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
const auto& bLiquidSatEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPvt().saturatedInverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_ex, p_ex, salt_ex) :
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_ex, p_ex);
const auto& densityLiquidEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
const auto& densityLiquidEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPvt().waterReferenceDensity(intQuantsEx.pvtRegionIndex()) :
FluidSystem::oilPvt().oilReferenceDensity(intQuantsEx.pvtRegionIndex());
const auto rho_ex = bLiquidEx * densityLiquidEx;
const auto rho_sat_ex = bLiquidSatEx
* (densityLiquidEx + rssat_ex * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, intQuantsEx.pvtRegionIndex()));
//rho difference approximation
const auto delta_rho = (rho_sat_ex + rho_sat_in - rho_in -rho_ex)/2;
const auto pressure_difference_convective_mixing = delta_rho * distZg;
@ -290,7 +293,6 @@ public:
unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
const auto& down = (downIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
unsigned globalDownIndex = (downIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
const auto& Rsup = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
up.fluidState().Rsw() :
@ -298,13 +300,11 @@ public:
const auto& Rsdown = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
down.fluidState().Rsw() :
down.fluidState().Rs();
const auto& RsSat = FluidSystem::saturatedDissolutionFactor(up.fluidState(),
liquidPhaseIdx,
up.pvtRegionIndex(),
SoMax);
const Evaluation& transMult = up.rockCompTransMultiplier();
Evaluation sg = up.fluidState().saturation(FluidSystem::gasPhaseIdx);
Evaluation X = (Rsup - RsSat * sg) / (RsSat * ( 1.0 - sg));
if ( (X > info.Psi_[up.pvtRegionIndex()] || Rsdown > 0) ) {
@ -312,13 +312,13 @@ public:
const auto& visc = up.fluidState().viscosity(liquidPhaseIdx);
// what will be the flux when muliplied with trans_mob
const auto convectiveFlux = -trans*transMult*info.Xhi_[up.pvtRegionIndex()]*invB*pressure_difference_convective_mixing*Rsup/(visc*faceArea);
unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
if (globalUpIndex == globalIndexIn)
flux[conti0EqIdx + activeGasCompIdx] += convectiveFlux;
else
flux[conti0EqIdx + activeGasCompIdx] += Opm::getValue(convectiveFlux);
}
}
}
};
};

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

@ -38,6 +38,7 @@
#include "blackoildiffusionmodule.hh"
#include "blackoildispersionmodule.hh"
#include "blackoilmicpmodules.hh"
#include "blackoilconvectivemixingmodule.hh"
#include <opm/common/TimingMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
@ -132,6 +133,7 @@ class BlackOilIntensiveQuantities
using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
using BrineModule = BlackOilBrineModule<TypeTag>;
using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
public:
@ -304,8 +306,7 @@ public:
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
SoMax = max(fluidState_.saturation(oilPhaseIdx),
problem.maxOilSaturation(globalSpaceIdx));
}
// take the meaning of the switching primary variable into account for the gas
// and oil phase compositions
if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
@ -315,15 +316,14 @@ public:
if (FluidSystem::enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
FluidSystem::saturatedDissolutionFactor(fluidState_,
oilPhaseIdx,
pvtRegionIdx,
SoMax);
oilPhaseIdx,
pvtRegionIdx,
SoMax);
fluidState_.setRs(min(RsMax, RsSat));
}
else if constexpr (compositionSwitchEnabled)
fluidState_.setRs(0.0);
}
}
if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
fluidState_.setRv(Rv);
@ -406,7 +406,8 @@ public:
rho = fluidState_.invB(waterPhaseIdx);
rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableDissolvedGasInWater()) {
if(!enableConvectiveMixing) {
bool ConvectiveMixingActive = ConvectiveMixingModule::active(elemCtx);
if(!ConvectiveMixingActive) {
rho +=
fluidState_.invB(waterPhaseIdx) *
fluidState_.Rsw() *
@ -437,8 +438,9 @@ public:
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
rho = fluidState_.invB(oilPhaseIdx);
rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
if (FluidSystem::enableDissolvedGas()) {
if(!enableConvectiveMixing) {
if (FluidSystem::enableDissolvedGas()) {
bool ConvectiveMixingActive = ConvectiveMixingModule::active(elemCtx);
if(!ConvectiveMixingActive) {
rho +=
fluidState_.invB(oilPhaseIdx) *
fluidState_.Rs() *

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

@ -94,7 +94,6 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
using MICPModule = BlackOilMICPModule<TypeTag>;
using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
public:
/*!
* \copydoc FvBaseLocalResidual::computeStorage
@ -214,7 +213,6 @@ public:
else
evalPhaseFluxes_<Scalar>(flux, phaseIdx, pvtRegionIdx, extQuants, up.fluidState());
}
// deal with solvents (if present)
SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
@ -238,7 +236,7 @@ public:
DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
// deal with convective mixing (if present)
// deal with convective mixing (if enabled)
ConvectiveMixingModule::addConvectiveMixingFlux(flux,elemCtx, scvfIdx, timeIdx);
}

5
opm/models/blackoil/blackoillocalresidualtpfa.hh
View File

@ -409,7 +409,7 @@ public:
}
}
}
}
// deal with solvents (if present)
static_assert(!enableSolvent, "Relevant computeFlux() method must be implemented for this module before enabling.");
@ -470,7 +470,6 @@ public:
moduleParams.convectiveMixingModuleParam);
}
// deal with diffusion (if present). opm-models expects per area flux (added in the tmpdiffusivity).
if constexpr(enableDiffusion){
typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
@ -645,7 +644,7 @@ public:
RateVector& bdyFlux,
const BoundaryConditionData& bdyInfo,
const IntensiveQuantities& insideIntQuants,
unsigned globalSpaceIdx)
[[maybe_unused]] unsigned globalSpaceIdx)
{
OPM_TIMEBLOCK_LOCAL(computeBoundaryThermal);
// only heat is allowed to flow through this boundary

1
opm/models/discretization/common/tpfalinearizer.hh
View File

@ -650,7 +650,6 @@ public:
return;
}
const unsigned int numCells = model_().numTotalDof();
#ifdef _OPENMP
#pragma omp parallel for
#endif