28 #ifndef EWOMS_DARCY_FLUX_MODULE_HH
29 #define EWOMS_DARCY_FLUX_MODULE_HH
34 #include <dune/common/fvector.hh>
35 #include <dune/common/fmatrix.hh>
40 namespace Properties {
44 template <
class TypeTag>
47 template <
class TypeTag>
50 template <
class TypeTag>
57 template <
class TypeTag>
76 template <
class TypeTag>
77 class DarcyBaseProblem
84 template <
class TypeTag>
85 class DarcyIntensiveQuantities
87 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
89 void update_(const ElementContext &elemCtx,
int dofIdx,
int timeIdx)
124 template <
class TypeTag>
125 class DarcyExtensiveQuantities
127 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
128 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
129 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
130 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
131 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
132 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
133 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
135 enum { dimWorld = GridView::dimensionworld };
138 typedef typename Opm::MathToolbox<Evaluation> Toolbox;
139 typedef Dune::FieldVector<Evaluation, dimWorld> EvalDimVector;
140 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
141 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
197 const auto& gradCalc = elemCtx.gradientCalculator();
200 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
201 const auto &faceNormal = scvf.normal();
207 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
208 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
225 const auto& gIn = elemCtx.problem().gravity(elemCtx,
interiorDofIdx_, timeIdx);
226 const auto& gEx = elemCtx.problem().gravity(elemCtx,
exteriorDofIdx_, timeIdx);
228 const auto &intQuantsIn = elemCtx.intensiveQuantities(
interiorDofIdx_, timeIdx);
229 const auto &intQuantsEx = elemCtx.intensiveQuantities(
exteriorDofIdx_, timeIdx);
233 const auto &posFace = scvf.integrationPos();
236 DimVector distVecIn(posIn);
237 DimVector distVecEx(posEx);
238 DimVector distVecTotal(posEx);
240 distVecIn -= posFace;
241 distVecEx -= posFace;
242 distVecTotal -= posIn;
243 Scalar absDistTotalSquared = distVecTotal.two_norm2();
244 for (
int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
245 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
249 auto rhoIn = intQuantsIn.fluidState().density(phaseIdx);
250 auto pStatIn = - rhoIn*(gIn*distVecIn);
254 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
255 Scalar pStatEx = - rhoEx*(gEx*distVecEx);
261 EvalDimVector f(distVecTotal);
262 f *= (pStatEx - pStatIn)/absDistTotalSquared;
268 if (!std::isfinite(Toolbox::value(
potentialGrad_[phaseIdx][i]))) {
269 OPM_THROW(Opm::NumericalProblem,
270 "Non-finite potential gradient for phase '"
271 << FluidSystem::phaseName(phaseIdx) <<
"'");
277 Valgrind::SetUndefined(
K_);
278 elemCtx.problem().intersectionIntrinsicPermeability(
K_, elemCtx, faceIdx, timeIdx);
279 Valgrind::CheckDefined(
K_);
281 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
282 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
288 Evaluation tmp = Toolbox::createConstant(0.0);
289 for (
unsigned i = 0; i < faceNormal.size(); ++i)
301 const auto &up = elemCtx.intensiveQuantities(
upstreamDofIdx_[phaseIdx], timeIdx);
307 mobility_[phaseIdx] = up.mobility(phaseIdx);
309 mobility_[phaseIdx] = Toolbox::value(up.mobility(phaseIdx));
319 template <
class Flu
idState>
323 const FluidState& fluidState,
324 const typename FluidSystem::ParameterCache ¶mCache)
326 const auto& gradCalc = elemCtx.gradientCalculator();
330 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
331 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
344 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
349 const auto &intQuantsIn = elemCtx.intensiveQuantities(
interiorDofIdx_, timeIdx);
350 K_ = intQuantsIn.intrinsicPermeability();
356 const auto& gIn = elemCtx.problem().gravity(elemCtx,
interiorDofIdx_, timeIdx);
358 const auto& posFace = scvf.integrationPos();
361 DimVector distVecIn(posIn);
362 distVecIn -= posFace;
363 Scalar absDist = distVecIn.two_norm();
364 Scalar gTimesDist = gIn*distVecIn;
366 for (
int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
367 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
371 Evaluation rhoIn = intQuantsIn.fluidState().density(phaseIdx);
372 Evaluation pStatIn = - gTimesDist*rhoIn;
374 Valgrind::CheckDefined(pStatIn);
382 EvalDimVector f(distVecIn);
383 f *= - pStatIn/absDist;
390 if (!std::isfinite(Toolbox::value(
potentialGrad_[phaseIdx][i]))) {
391 OPM_THROW(Opm::NumericalProblem,
392 "Non finite potential gradient for phase '"
393 << FluidSystem::phaseName(phaseIdx) <<
"'");
400 const auto &faceNormal = scvf.normal();
402 const auto &matParams = elemCtx.problem().materialLawParams(elemCtx,
interiorDofIdx_, timeIdx);
404 Scalar kr[numPhases];
405 MaterialLaw::relativePermeabilities(kr, matParams, fluidState);
406 Valgrind::CheckDefined(kr);
408 for (
int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
409 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
412 Evaluation tmp = Toolbox::createConstant(0.0);
413 for (
unsigned i = 0; i < faceNormal.size(); ++i)
428 kr[phaseIdx] / FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
430 mobility_[phaseIdx] = intQuantsIn.mobility(phaseIdx);
431 Valgrind::CheckDefined(
mobility_[phaseIdx]);
443 const auto &scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
444 const DimVector &normal = scvf.normal();
445 Valgrind::CheckDefined(normal);
447 for (
int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
449 volumeFlux_[phaseIdx] = Toolbox::createConstant(0.0);
450 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
454 asImp_().calculateFilterVelocity_(phaseIdx);
456 volumeFlux_[phaseIdx] = Toolbox::createConstant(0.0);
457 for (
unsigned i = 0; i < normal.size(); ++i)
472 const auto &scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
473 const DimVector &normal = scvf.normal();
474 Valgrind::CheckDefined(normal);
476 for (
int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
477 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
479 volumeFlux_[phaseIdx] = Toolbox::createConstant(0);
483 asImp_().calculateFilterVelocity_(phaseIdx);
485 volumeFlux_[phaseIdx] = Toolbox::createConstant(0.0);
486 for (
unsigned i = 0; i < normal.size(); ++i)
494 assert(std::isfinite(Toolbox::value(
mobility_[phaseIdx])));
495 for (
unsigned i = 0; i <
K_.M(); ++ i)
496 for (
unsigned j = 0; j <
K_.N(); ++ j)
497 assert(std::isfinite(
K_[i][j]));
510 Implementation &asImp_()
511 {
return *
static_cast<Implementation*
>(
this); }
513 const Implementation &asImp_()
const
514 {
return *
static_cast<const Implementation*
>(
this); }
const EvalDimVector & filterVelocity(int phaseIdx) const
Return the filter velocity of a fluid phase at the face's integration point [m/s].
Definition: darcyfluxmodule.hh:166
short interiorDofIdx_
Definition: darcyfluxmodule.hh:521
void setPhaseIndex(short phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:155
const DimMatrix & intrinsicPermability() const
Returns the intrinsic permeability tensor for a given sub-control volume face.
Definition: darcyfluxmodule.hh:148
Specifies a flux module which uses the Darcy relation.
Definition: darcyfluxmodule.hh:58
void calculateFluxes_(const ElementContext &elemCtx, int scvfIdx, int timeIdx)
Calculate the volumetric fluxes of all phases.
Definition: darcyfluxmodule.hh:441
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
short upstreamDofIdx_[numPhases]
Definition: darcyfluxmodule.hh:523
const EvalDimVector & potentialGrad(int phaseIdx) const
Return the pressure potential gradient of a fluid phase at the face's integration point [Pa/m]...
Definition: darcyfluxmodule.hh:157
void update_(const ElementContext &elemCtx, int dofIdx, int timeIdx)
Definition: darcyfluxmodule.hh:89
const Evaluation & volumeFlux(int phaseIdx) const
Return the volume flux of a fluid phase at the face's integration point .
Definition: darcyfluxmodule.hh:178
DarcyIntensiveQuantities< TypeTag > FluxIntensiveQuantities
Definition: darcyfluxmodule.hh:60
NEW_PROP_TAG(Grid)
The type of the DUNE grid.
Provides the defaults for the parameters required by the Darcy velocity approach. ...
Definition: darcyfluxmodule.hh:51
void setPhaseIndex(short phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:101
short upstreamIndex_(int phaseIdx) const
Definition: darcyfluxmodule.hh:182
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:125
short downstreamIndex_(int phaseIdx) const
Definition: darcyfluxmodule.hh:185
Evaluation mobility_[numPhases]
Definition: darcyfluxmodule.hh:527
void calculateBoundaryGradients_(const ElementContext &elemCtx, int boundaryFaceIdx, int timeIdx, const FluidState &fluidState, const typename FluidSystem::ParameterCache ¶mCache)
Calculate the gradients at the grid boundary which are required to determine the volumetric fluxes...
Definition: darcyfluxmodule.hh:320
Definition: baseauxiliarymodule.hh:35
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:78
void calculateBoundaryFluxes_(const ElementContext &elemCtx, int boundaryFaceIdx, int timeIdx)
Calculate the volumetric fluxes at a boundary face of all fluid phases.
Definition: darcyfluxmodule.hh:468
Evaluation volumeFlux_[numPhases]
Definition: darcyfluxmodule.hh:534
DimMatrix K_
Definition: darcyfluxmodule.hh:518
Provides the Darcy flux module.
Definition: darcyfluxmodule.hh:48
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: darcyfluxmodule.hh:67
void calculateFilterVelocity_(int phaseIdx)
Definition: darcyfluxmodule.hh:491
This method contains all callback classes for quantities that are required by some extensive quantiti...
Provides the intensive quantities for the Darcy flux module.
Definition: darcyfluxmodule.hh:45
EvalDimVector filterVelocity_[numPhases]
Definition: darcyfluxmodule.hh:530
DarcyExtensiveQuantities< TypeTag > FluxExtensiveQuantities
Definition: darcyfluxmodule.hh:61
void calculateGradients_(const ElementContext &elemCtx, int faceIdx, int timeIdx)
Calculate the gradients which are required to determine the volumetric fluxes.
Definition: darcyfluxmodule.hh:193
short downstreamDofIdx_[numPhases]
Definition: darcyfluxmodule.hh:524
Defines the common properties required by the porous medium multi-phase models.
EvalDimVector potentialGrad_[numPhases]
Definition: darcyfluxmodule.hh:537
DarcyBaseProblem< TypeTag > FluxBaseProblem
Definition: darcyfluxmodule.hh:62
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:95
short exteriorDofIdx_
Definition: darcyfluxmodule.hh:522