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>
44#include <opm/models/discretization/common/fvbaseproperties.hh>
45#include <opm/models/blackoil/blackoilproperties.hh>
46#include <opm/models/utils/signum.hh>
52template <
class TypeTag>
53class NewTranIntensiveQuantities;
55template <
class TypeTag>
56class NewTranExtensiveQuantities;
58template <
class TypeTag>
59class NewTranBaseProblem;
65template <
class TypeTag>
84template <
class TypeTag>
92template <
class TypeTag>
95 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
97 void update_(
const ElementContext&,
unsigned,
unsigned)
105template <
class TypeTag>
108 using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
110 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
111 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
112 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
113 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
114 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
115 using GridView = GetPropType<TypeTag, Properties::GridView>;
116 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
118 enum { dimWorld = GridView::dimensionworld };
119 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
120 enum { numPhases = FluidSystem::numPhases };
121 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
122 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
123 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
125 using Toolbox = MathToolbox<Evaluation>;
126 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
127 using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
128 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
136 throw std::invalid_argument(
"The ECL transmissibility module does not provide an explicit intrinsic permeability");
147 throw std::invalid_argument(
"The ECL transmissibility module does not provide explicit potential gradients");
157 {
return pressureDifference_[phaseIdx]; }
167 throw std::invalid_argument(
"The ECL transmissibility module does not provide explicit filter velocities");
180 {
return volumeFlux_[phaseIdx]; }
192 assert(phaseIdx < numPhases);
194 return upIdx_[phaseIdx];
206 assert(phaseIdx < numPhases);
208 return dnIdx_[phaseIdx];
211 void updateSolvent(
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
212 { asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx); }
214 void updatePolymer(
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
215 { asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx); }
220 std::array<short, numPhases>& dnIdx,
222 Evaluation (&pressureDifferences)[numPhases],
223 const ElementContext& elemCtx,
227 const auto& problem = elemCtx.problem();
228 const auto& stencil = elemCtx.stencil(timeIdx);
229 const auto& scvf = stencil.interiorFace(scvfIdx);
230 unsigned interiorDofIdx = scvf.interiorIndex();
231 unsigned exteriorDofIdx = scvf.exteriorIndex();
233 assert(interiorDofIdx != exteriorDofIdx);
235 unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
236 unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
237 Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
238 Scalar faceArea = scvf.area();
239 Scalar thpres = problem.thresholdPressure(I, J);
244 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
246 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
247 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
254 Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
255 Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
259 Scalar distZ = zIn - zEx;
261 Scalar Vin = elemCtx.dofVolume(interiorDofIdx, 0);
262 Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, 0);
264 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
265 if (!FluidSystem::phaseIsActive(phaseIdx))
269 pressureDifferences[phaseIdx],
281 if (pressureDifferences[phaseIdx] == 0) {
286 const bool upwindIsInterior = (
static_cast<unsigned>(upIdx[phaseIdx]) == interiorDofIdx);
287 const IntensiveQuantities& up = upwindIsInterior ? intQuantsIn : intQuantsEx;
289 const Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2;
291 const auto& materialLawManager = problem.materialLawManager();
292 FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown;
293 if (materialLawManager->hasDirectionalRelperms()) {
294 facedir = scvf.faceDirFromDirId();
296 if (upwindIsInterior)
298 pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea);
301 pressureDifferences[phaseIdx]*
302 (Toolbox::value(up.mobility(phaseIdx, facedir))*transMult*(-trans/faceArea));
306 template<
class EvalType>
310 const IntensiveQuantities& intQuantsIn,
311 const IntensiveQuantities& intQuantsEx,
312 const unsigned phaseIdx,
313 const unsigned interiorDofIdx,
314 const unsigned exteriorDofIdx,
317 const unsigned globalIndexIn,
318 const unsigned globalIndexEx,
326 if (intQuantsIn.mobility(phaseIdx) <= 0.0 &&
327 intQuantsEx.mobility(phaseIdx) <= 0.0)
329 upIdx = interiorDofIdx;
330 dnIdx = exteriorDofIdx;
337 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
338 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
339 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
341 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
342 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
345 pressureExterior += Toolbox::value(rhoAvg)*(distZg);
347 pressureExterior += rhoAvg*(distZg);
356 upIdx = exteriorDofIdx;
357 dnIdx = interiorDofIdx;
360 upIdx = interiorDofIdx;
361 dnIdx = exteriorDofIdx;
367 upIdx = interiorDofIdx;
368 dnIdx = exteriorDofIdx;
370 else if (Vin < Vex) {
371 upIdx = exteriorDofIdx;
372 dnIdx = interiorDofIdx;
378 if (globalIndexIn < globalIndexEx) {
379 upIdx = interiorDofIdx;
380 dnIdx = exteriorDofIdx;
383 upIdx = exteriorDofIdx;
384 dnIdx = interiorDofIdx;
412 Valgrind::SetUndefined(*
this);
420 template <
class Flu
idState>
424 const FluidState& exFluidState)
426 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx);
427 const Scalar faceArea = scvf.area();
428 const Scalar zEx = scvf.integrationPos()[dimWorld - 1];
429 const auto& problem = elemCtx.problem();
430 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(0, timeIdx);
431 const auto& intQuantsIn = elemCtx.intensiveQuantities(0, timeIdx);
443 pressureDifference_);
448 if constexpr (enableSolvent) {
449 if (upIdx_[gasPhaseIdx] == 0) {
450 const Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, scvfIdx);
451 const Scalar transModified = trans * Toolbox::value(intQuantsIn.rockCompTransMultiplier());
452 const auto solventFlux = pressureDifference_[gasPhaseIdx] * intQuantsIn.mobility(gasPhaseIdx) * (-transModified/faceArea);
453 asImp_().setSolventVolumeFlux(solventFlux);
455 asImp_().setSolventVolumeFlux(0.0);
464 template <
class Problem,
class Flu
idState,
class EvaluationContainer>
466 const unsigned globalSpaceIdx,
467 const IntensiveQuantities& intQuantsIn,
468 const unsigned bfIdx,
469 const double faceArea,
471 const FluidState& exFluidState,
472 std::array<short, numPhases>& upIdx,
473 std::array<short, numPhases>& dnIdx,
478 bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions();
479 if (!enableBoundaryMassFlux)
482 Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx);
487 Scalar g = problem.gravity()[dimWorld - 1];
494 Scalar zIn = problem.dofCenterDepth(globalSpaceIdx);
498 Scalar distZ = zIn - zEx;
500 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
501 if (!FluidSystem::phaseIsActive(phaseIdx))
506 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
507 const auto& rhoEx = exFluidState.density(phaseIdx);
508 Evaluation rhoAvg = (rhoIn + rhoEx)/2;
510 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
511 Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
512 pressureExterior += rhoAvg*(distZ*g);
520 const unsigned interiorDofIdx = 0;
522 upIdx[phaseIdx] = -1;
523 dnIdx[phaseIdx] = interiorDofIdx;
526 upIdx[phaseIdx] = interiorDofIdx;
527 dnIdx[phaseIdx] = -1;
530 Evaluation transModified = trans;
532 if (upIdx[phaseIdx] == interiorDofIdx) {
536 const auto& up = intQuantsIn;
539 const Scalar transMult = Toolbox::value(up.rockCompTransMultiplier());
540 transModified *= transMult;
548 const auto& matParams = problem.materialLawParams(globalSpaceIdx);
549 std::array<typename FluidState::Scalar,numPhases> kr;
550 MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
552 const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
571 Implementation& asImp_()
572 {
return *
static_cast<Implementation*
>(
this); }
574 const Implementation& asImp_()
const
575 {
return *
static_cast<const Implementation*
>(
this); }
578 Evaluation volumeFlux_[numPhases];
582 Evaluation pressureDifference_[numPhases];
585 std::array<short, numPhases> upIdx_;
586 std::array<short, numPhases> dnIdx_;
Provides the defaults for the parameters required by the transmissibility based volume flux calculati...
Definition: NewTranFluxModule.hpp:86
Provides the ECL flux module.
Definition: NewTranFluxModule.hpp:107
const Evaluation & volumeFlux(unsigned phaseIdx) const
Return the volume flux of a fluid phase at the face's integration point .
Definition: NewTranFluxModule.hpp:179
void updateSolvent(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: NewTranFluxModule.hpp:211
void calculateBoundaryFluxes_(const ElementContext &, unsigned, unsigned)
Definition: NewTranFluxModule.hpp:567
const DimMatrix & intrinsicPermeability() const
Return the intrinsic permeability tensor at a face [m^2].
Definition: NewTranFluxModule.hpp:134
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)
Definition: NewTranFluxModule.hpp:307
void calculateBoundaryGradients_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx, const FluidState &exFluidState)
Update the required gradients for boundary faces.
Definition: NewTranFluxModule.hpp:421
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:465
unsigned upstreamIndex_(unsigned phaseIdx) const
Returns the local index of the degree of freedom in which is in upstream direction.
Definition: NewTranFluxModule.hpp:190
void updatePolymer(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: NewTranFluxModule.hpp:214
void calculateFluxes_(const ElementContext &, unsigned, unsigned)
Update the volumetric fluxes for all fluid phases on the interior faces of the context.
Definition: NewTranFluxModule.hpp:564
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:145
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:219
unsigned downstreamIndex_(unsigned phaseIdx) const
Returns the local index of the degree of freedom in which is in downstream direction.
Definition: NewTranFluxModule.hpp:204
const EvalDimVector & filterVelocity(unsigned) const
Return the filter velocity of a fluid phase at the face's integration point [m/s].
Definition: NewTranFluxModule.hpp:165
const Evaluation & pressureDifference(unsigned phaseIdx) const
Return the gravity corrected pressure difference between the interior and the exterior of a face.
Definition: NewTranFluxModule.hpp:156
void calculateGradients_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Update the required gradients for interior faces.
Definition: NewTranFluxModule.hpp:410
Provides the intensive quantities for the ECL flux module.
Definition: NewTranFluxModule.hpp:94
void update_(const ElementContext &, unsigned, unsigned)
Definition: NewTranFluxModule.hpp:97
Definition: BlackoilPhases.hpp:27
Specifies a flux module which uses ECL transmissibilities.
Definition: NewTranFluxModule.hpp:67
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: NewTranFluxModule.hpp:75