31#ifndef OPM_NEWTRAN_FLUX_MODULE_HPP
32#define OPM_NEWTRAN_FLUX_MODULE_HPP
34#include <dune/common/fvector.hh>
35#include <dune/common/fmatrix.hh>
37#include <opm/common/OpmLog/OpmLog.hpp>
39#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
41#include <opm/material/common/MathToolbox.hpp>
42#include <opm/material/common/Valgrind.hpp>
53template <
class TypeTag>
54class NewTranIntensiveQuantities;
56template <
class TypeTag>
57class NewTranExtensiveQuantities;
59template <
class TypeTag>
60class NewTranBaseProblem;
66template <
class TypeTag>
85template <
class TypeTag>
93template <
class TypeTag>
98 void update_(
const ElementContext&,
unsigned,
unsigned)
106template <
class TypeTag>
119 enum { dimWorld = GridView::dimensionworld };
120 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
121 enum { numPhases = FluidSystem::numPhases };
122 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
123 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
124 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
126 static constexpr bool enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>();
129 using Toolbox = MathToolbox<Evaluation>;
130 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
131 using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
132 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
142 throw std::invalid_argument(
"The ECL transmissibility module does not provide an explicit intrinsic permeability");
153 throw std::invalid_argument(
"The ECL transmissibility module does not provide explicit potential gradients");
163 {
return pressureDifference_[phaseIdx]; }
173 throw std::invalid_argument(
"The ECL transmissibility module does not provide explicit filter velocities");
186 {
return volumeFlux_[phaseIdx]; }
198 assert(phaseIdx < numPhases);
200 return upIdx_[phaseIdx];
212 assert(phaseIdx < numPhases);
214 return dnIdx_[phaseIdx];
217 void updateSolvent(
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
218 { asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx); }
220 void updatePolymer(
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
221 { asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx); }
226 std::array<short, numPhases>& dnIdx,
228 Evaluation (&pressureDifferences)[numPhases],
229 const ElementContext& elemCtx,
233 const auto& problem = elemCtx.problem();
234 const auto& stencil = elemCtx.stencil(timeIdx);
235 const auto& scvf = stencil.interiorFace(scvfIdx);
236 unsigned interiorDofIdx = scvf.interiorIndex();
237 unsigned exteriorDofIdx = scvf.exteriorIndex();
239 assert(interiorDofIdx != exteriorDofIdx);
241 unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
242 unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
243 Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
244 Scalar faceArea = scvf.area();
245 Scalar thpres = problem.thresholdPressure(I, J);
250 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
252 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
253 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
260 Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
261 Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
265 Scalar distZ = zIn - zEx;
267 Scalar Vin = elemCtx.dofVolume(interiorDofIdx, 0);
268 Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, 0);
270 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
271 if (!FluidSystem::phaseIsActive(phaseIdx))
275 pressureDifferences[phaseIdx],
287 problem.moduleParams());
289 const bool upwindIsInterior = (
static_cast<unsigned>(upIdx[phaseIdx]) == interiorDofIdx);
290 const IntensiveQuantities& up = upwindIsInterior ? intQuantsIn : intQuantsEx;
292 const Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2;
294 const auto& materialLawManager = problem.materialLawManager();
295 FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown;
296 if (materialLawManager->hasDirectionalRelperms()) {
297 facedir = scvf.faceDirFromDirId();
299 if (upwindIsInterior)
301 pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea);
304 pressureDifferences[phaseIdx]*
305 (Toolbox::value(up.mobility(phaseIdx, facedir))*transMult*(-trans/faceArea));
309 template<
class EvalType>
313 const IntensiveQuantities& intQuantsIn,
314 const IntensiveQuantities& intQuantsEx,
315 const unsigned phaseIdx,
316 const unsigned interiorDofIdx,
317 const unsigned exteriorDofIdx,
320 const unsigned globalIndexIn,
321 const unsigned globalIndexEx,
324 const ModuleParams& moduleParams)
329 if (intQuantsIn.mobility(phaseIdx) <= 0.0 &&
330 intQuantsEx.mobility(phaseIdx) <= 0.0)
332 upIdx = interiorDofIdx;
333 dnIdx = exteriorDofIdx;
340 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
341 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
342 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
344 if constexpr(enableConvectiveMixing) {
345 ConvectiveMixingModule::modifyAvgDensity(rhoAvg, intQuantsIn, intQuantsEx, phaseIdx, moduleParams.convectiveMixingModuleParam);
348 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
349 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
352 pressureExterior += Toolbox::value(rhoAvg)*(distZg);
354 pressureExterior += rhoAvg*(distZg);
363 upIdx = exteriorDofIdx;
364 dnIdx = interiorDofIdx;
367 upIdx = interiorDofIdx;
368 dnIdx = exteriorDofIdx;
374 upIdx = interiorDofIdx;
375 dnIdx = exteriorDofIdx;
377 else if (Vin < Vex) {
378 upIdx = exteriorDofIdx;
379 dnIdx = interiorDofIdx;
385 if (globalIndexIn < globalIndexEx) {
386 upIdx = interiorDofIdx;
387 dnIdx = exteriorDofIdx;
390 upIdx = exteriorDofIdx;
391 dnIdx = interiorDofIdx;
419 Valgrind::SetUndefined(*
this);
427 template <
class Flu
idState>
431 const FluidState& exFluidState)
433 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx);
434 const Scalar faceArea = scvf.area();
435 const Scalar zEx = scvf.integrationPos()[dimWorld - 1];
436 const auto& problem = elemCtx.problem();
437 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(0, timeIdx);
438 const auto& intQuantsIn = elemCtx.intensiveQuantities(0, timeIdx);
450 pressureDifference_);
455 if constexpr (enableSolvent) {
456 if (upIdx_[gasPhaseIdx] == 0) {
457 const Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, scvfIdx);
458 const Scalar transModified = trans * Toolbox::value(intQuantsIn.rockCompTransMultiplier());
459 const auto solventFlux = pressureDifference_[gasPhaseIdx] * intQuantsIn.mobility(gasPhaseIdx) * (-transModified/faceArea);
460 asImp_().setSolventVolumeFlux(solventFlux);
462 asImp_().setSolventVolumeFlux(0.0);
471 template <
class Problem,
class Flu
idState,
class EvaluationContainer>
473 const unsigned globalSpaceIdx,
474 const IntensiveQuantities& intQuantsIn,
475 const unsigned bfIdx,
476 const double faceArea,
478 const FluidState& exFluidState,
479 std::array<short, numPhases>& upIdx,
480 std::array<short, numPhases>& dnIdx,
485 bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions();
486 if (!enableBoundaryMassFlux)
489 Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx);
494 Scalar g = problem.gravity()[dimWorld - 1];
501 Scalar zIn = problem.dofCenterDepth(globalSpaceIdx);
505 Scalar distZ = zIn - zEx;
507 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
508 if (!FluidSystem::phaseIsActive(phaseIdx))
513 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
514 const auto& rhoEx = exFluidState.density(phaseIdx);
515 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
517 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
518 Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
519 pressureExterior += rhoAvg*(distZ*g);
527 const unsigned interiorDofIdx = 0;
529 upIdx[phaseIdx] = -1;
530 dnIdx[phaseIdx] = interiorDofIdx;
533 upIdx[phaseIdx] = interiorDofIdx;
534 dnIdx[phaseIdx] = -1;
537 Evaluation transModified = trans;
539 if (upIdx[phaseIdx] == interiorDofIdx) {
543 const auto& up = intQuantsIn;
546 const Scalar transMult = Toolbox::value(up.rockCompTransMultiplier());
547 transModified *= transMult;
555 const auto& matParams = problem.materialLawParams(globalSpaceIdx);
556 std::array<typename FluidState::Scalar,numPhases> kr;
557 MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
559 const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
578 Implementation& asImp_()
579 {
return *
static_cast<Implementation*
>(
this); }
581 const Implementation& asImp_()
const
582 {
return *
static_cast<const Implementation*
>(
this); }
585 Evaluation volumeFlux_[numPhases];
589 Evaluation pressureDifference_[numPhases];
592 std::array<short, numPhases> upIdx_;
593 std::array<short, numPhases> dnIdx_;
Declares the properties required by the black oil model.
Definition: blackoilconvectivemixingmodule.hh:58
Provides the defaults for the parameters required by the transmissibility based volume flux calculati...
Definition: NewTranFluxModule.hpp:87
Provides the ECL flux module.
Definition: NewTranFluxModule.hpp:108
const Evaluation & volumeFlux(unsigned phaseIdx) const
Return the volume flux of a fluid phase at the face's integration point .
Definition: NewTranFluxModule.hpp:185
void updateSolvent(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: NewTranFluxModule.hpp:217
void calculateBoundaryFluxes_(const ElementContext &, unsigned, unsigned)
Definition: NewTranFluxModule.hpp:574
const DimMatrix & intrinsicPermeability() const
Return the intrinsic permeability tensor at a face [m^2].
Definition: NewTranFluxModule.hpp:140
static void calculatePhasePressureDiff_(short &upIdx, short &dnIdx, EvalType &pressureDifference, const IntensiveQuantities &intQuantsIn, const IntensiveQuantities &intQuantsEx, const unsigned phaseIdx, const unsigned interiorDofIdx, const unsigned exteriorDofIdx, const Scalar Vin, const Scalar Vex, const unsigned globalIndexIn, const unsigned globalIndexEx, const Scalar distZg, const Scalar thpres, const ModuleParams &moduleParams)
Definition: NewTranFluxModule.hpp:310
void calculateBoundaryGradients_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx, const FluidState &exFluidState)
Update the required gradients for boundary faces.
Definition: NewTranFluxModule.hpp:428
static void calculateBoundaryGradients_(const Problem &problem, const unsigned globalSpaceIdx, const IntensiveQuantities &intQuantsIn, const unsigned bfIdx, const double faceArea, const double zEx, const FluidState &exFluidState, std::array< short, numPhases > &upIdx, std::array< short, numPhases > &dnIdx, EvaluationContainer &volumeFlux, EvaluationContainer &pressureDifference)
Update the required gradients for boundary faces.
Definition: NewTranFluxModule.hpp:472
unsigned upstreamIndex_(unsigned phaseIdx) const
Returns the local index of the degree of freedom in which is in upstream direction.
Definition: NewTranFluxModule.hpp:196
void updatePolymer(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: NewTranFluxModule.hpp:220
void calculateFluxes_(const ElementContext &, unsigned, unsigned)
Update the volumetric fluxes for all fluid phases on the interior faces of the context.
Definition: NewTranFluxModule.hpp:571
const EvalDimVector & potentialGrad(unsigned) const
Return the pressure potential gradient of a fluid phase at the face's integration point [Pa/m].
Definition: NewTranFluxModule.hpp:151
static void volumeAndPhasePressureDifferences(std::array< short, numPhases > &upIdx, std::array< short, numPhases > &dnIdx, Evaluation(&volumeFlux)[numPhases], Evaluation(&pressureDifferences)[numPhases], const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: NewTranFluxModule.hpp:225
unsigned downstreamIndex_(unsigned phaseIdx) const
Returns the local index of the degree of freedom in which is in downstream direction.
Definition: NewTranFluxModule.hpp:210
const EvalDimVector & filterVelocity(unsigned) const
Return the filter velocity of a fluid phase at the face's integration point [m/s].
Definition: NewTranFluxModule.hpp:171
const Evaluation & pressureDifference(unsigned phaseIdx) const
Return the gravity corrected pressure difference between the interior and the exterior of a face.
Definition: NewTranFluxModule.hpp:162
void calculateGradients_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Update the required gradients for interior faces.
Definition: NewTranFluxModule.hpp:417
Provides the intensive quantities for the ECL flux module.
Definition: NewTranFluxModule.hpp:95
void update_(const ElementContext &, unsigned, unsigned)
Definition: NewTranFluxModule.hpp:98
Declare the properties used by the infrastructure code of the finite volume discretizations.
Definition: blackoilboundaryratevector.hh:37
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:235
Definition: blackoillocalresidualtpfa.hh:133
Specifies a flux module which uses ECL transmissibilities.
Definition: NewTranFluxModule.hpp:68
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: NewTranFluxModule.hpp:76