30#ifndef EWOMS_DARCY_FLUX_MODULE_HH
31#define EWOMS_DARCY_FLUX_MODULE_HH
35#include <opm/common/Exceptions.hpp>
39#include <opm/material/common/Valgrind.hpp>
41#include <dune/common/fvector.hh>
42#include <dune/common/fmatrix.hh>
48template <
class TypeTag>
49class DarcyIntensiveQuantities;
51template <
class TypeTag>
52class DarcyExtensiveQuantities;
54template <
class TypeTag>
55class DarcyBaseProblem;
61template <
class TypeTag>
80template <
class TypeTag>
88template <
class TypeTag>
115template <
class TypeTag>
126 enum { dimWorld = GridView::dimensionworld };
127 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
129 using Toolbox = MathToolbox<Evaluation>;
130 using ParameterCache =
typename FluidSystem::template ParameterCache<Evaluation>;
131 using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
132 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
133 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
189 const auto& gradCalc = elemCtx.gradientCalculator();
192 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
193 const auto& faceNormal = scvf.normal();
195 unsigned i = scvf.interiorIndex();
196 unsigned j = scvf.exteriorIndex();
199 unsigned focusDofIdx = elemCtx.focusDofIndex();
202 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
203 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
217 if (Parameters::get<TypeTag, Properties::EnableGravity>()) {
220 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
221 const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
223 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
224 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
226 const auto& posIn = elemCtx.pos(i, timeIdx);
227 const auto& posEx = elemCtx.pos(j, timeIdx);
228 const auto& posFace = scvf.integrationPos();
231 DimVector distVecIn(posIn);
232 DimVector distVecEx(posEx);
233 DimVector distVecTotal(posEx);
235 distVecIn -= posFace;
236 distVecEx -= posFace;
237 distVecTotal -= posIn;
238 Scalar absDistTotalSquared = distVecTotal.two_norm2();
239 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
240 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
246 if (std::is_same<Scalar, Evaluation>::value ||
249 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
250 pStatIn = - rhoIn*(gIn*distVecIn);
253 Scalar rhoIn = Toolbox::value(intQuantsIn.fluidState().density(phaseIdx));
254 pStatIn = - rhoIn*(gIn*distVecIn);
261 if (std::is_same<Scalar, Evaluation>::value ||
264 const Evaluation& rhoEx = intQuantsEx.fluidState().density(phaseIdx);
265 pStatEx = - rhoEx*(gEx*distVecEx);
268 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
269 pStatEx = - rhoEx*(gEx*distVecEx);
276 Dune::FieldVector<Evaluation, dimWorld> f(distVecTotal);
277 f *= (pStatEx - pStatIn)/absDistTotalSquared;
280 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
283 for (
unsigned dimIdx = 0; dimIdx <
potentialGrad_[phaseIdx].size(); ++dimIdx) {
285 throw NumericalProblem(
"Non-finite potential gradient for phase '"
286 + std::string(FluidSystem::phaseName(phaseIdx))+
"'");
292 Valgrind::SetUndefined(
K_);
293 elemCtx.problem().intersectionIntrinsicPermeability(
K_, elemCtx, faceIdx, timeIdx);
294 Valgrind::CheckDefined(
K_);
296 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
297 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
303 Evaluation tmp = 0.0;
304 for (
unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
318 const auto& up = elemCtx.intensiveQuantities(
upstreamDofIdx_[phaseIdx], timeIdx);
320 mobility_[phaseIdx] = up.mobility(phaseIdx);
322 mobility_[phaseIdx] = Toolbox::value(up.mobility(phaseIdx));
332 template <
class Flu
idState>
334 unsigned boundaryFaceIdx,
336 const FluidState& fluidState)
338 const auto& gradCalc = elemCtx.gradientCalculator();
342 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
343 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
356 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
357 auto i = scvf.interiorIndex();
360 int focusDofIdx = elemCtx.focusDofIndex();
363 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
364 K_ = intQuantsIn.intrinsicPermeability();
367 if (Parameters::get<TypeTag, Properties::EnableGravity>()) {
370 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
371 const auto& posIn = elemCtx.pos(i, timeIdx);
372 const auto& posFace = scvf.integrationPos();
375 DimVector distVecIn(posIn);
376 distVecIn -= posFace;
377 Scalar absDistSquared = distVecIn.two_norm2();
378 Scalar gTimesDist = gIn*distVecIn;
380 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
381 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
385 Evaluation rhoIn = intQuantsIn.fluidState().density(phaseIdx);
386 Evaluation pStatIn = - gTimesDist*rhoIn;
388 Valgrind::CheckDefined(pStatIn);
394 EvalDimVector f(distVecIn);
395 f *= pStatIn / absDistSquared;
398 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
402 for (
unsigned dimIdx = 0; dimIdx <
potentialGrad_[phaseIdx].size(); ++dimIdx) {
404 throw NumericalProblem(
"Non-finite potential gradient for phase '"
405 + std::string(FluidSystem::phaseName(phaseIdx))+
"'");
412 const auto& faceNormal = scvf.normal();
414 const auto& matParams = elemCtx.problem().materialLawParams(elemCtx, i, timeIdx);
416 Scalar kr[numPhases];
417 MaterialLaw::relativePermeabilities(kr, matParams, fluidState);
418 Valgrind::CheckDefined(kr);
420 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
421 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
424 Evaluation tmp = 0.0;
425 for (
unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
441 kr[phaseIdx] / fluidState.viscosity(phaseIdx);
444 Toolbox::value(kr[phaseIdx])
445 / Toolbox::value(fluidState.viscosity(phaseIdx));
448 mobility_[phaseIdx] = Toolbox::value(intQuantsIn.mobility(phaseIdx));
450 mobility_[phaseIdx] = intQuantsIn.mobility(phaseIdx);
451 Valgrind::CheckDefined(
mobility_[phaseIdx]);
463 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
464 const DimVector& normal = scvf.normal();
465 Valgrind::CheckDefined(normal);
467 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
470 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
473 asImp_().calculateFilterVelocity_(phaseIdx);
477 for (
unsigned i = 0; i < normal.size(); ++i)
489 unsigned boundaryFaceIdx,
492 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
493 const DimVector& normal = scvf.normal();
494 Valgrind::CheckDefined(normal);
496 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
497 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
503 asImp_().calculateFilterVelocity_(phaseIdx);
506 for (
unsigned i = 0; i < normal.size(); ++i)
515 for (
unsigned i = 0; i <
K_.M(); ++ i)
516 for (
unsigned j = 0; j <
K_.N(); ++ j)
517 assert(std::isfinite(
K_[i][j]));
530 Implementation& asImp_()
531 {
return *
static_cast<Implementation*
>(
this); }
533 const Implementation& asImp_()
const
534 {
return *
static_cast<const Implementation*
>(
this); }
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:133
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:163
Provides the defaults for the parameters required by the Darcy velocity approach.
Definition: darcyfluxmodule.hh:82
Provides the Darcy flux module.
Definition: darcyfluxmodule.hh:117
Evaluation mobility_[numPhases]
Definition: darcyfluxmodule.hh:541
void calculateBoundaryGradients_(const ElementContext &elemCtx, unsigned boundaryFaceIdx, unsigned timeIdx, const FluidState &fluidState)
Calculate the gradients at the grid boundary which are required to determine the volumetric fluxes.
Definition: darcyfluxmodule.hh:333
short interiorDofIdx_
Definition: darcyfluxmodule.hh:556
const EvalDimVector & filterVelocity(unsigned phaseIdx) const
Return the filter velocity of a fluid phase at the face's integration point [m/s].
Definition: darcyfluxmodule.hh:158
void calculateBoundaryFluxes_(const ElementContext &elemCtx, unsigned boundaryFaceIdx, unsigned timeIdx)
Calculate the volumetric fluxes at a boundary face of all fluid phases.
Definition: darcyfluxmodule.hh:488
short upstreamIndex_(unsigned phaseIdx) const
Definition: darcyfluxmodule.hh:174
EvalDimVector potentialGrad_[numPhases]
Definition: darcyfluxmodule.hh:551
EvalDimVector filterVelocity_[numPhases]
Definition: darcyfluxmodule.hh:544
const EvalDimVector & potentialGrad(unsigned phaseIdx) const
Return the pressure potential gradient of a fluid phase at the face's integration point [Pa/m].
Definition: darcyfluxmodule.hh:149
short upstreamDofIdx_[numPhases]
Definition: darcyfluxmodule.hh:554
void calculateFluxes_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Calculate the volumetric fluxes of all phases.
Definition: darcyfluxmodule.hh:461
short downstreamDofIdx_[numPhases]
Definition: darcyfluxmodule.hh:555
const Evaluation & volumeFlux(unsigned phaseIdx) const
Return the volume flux of a fluid phase at the face's integration point .
Definition: darcyfluxmodule.hh:170
short downstreamIndex_(unsigned phaseIdx) const
Definition: darcyfluxmodule.hh:177
Evaluation volumeFlux_[numPhases]
Definition: darcyfluxmodule.hh:548
void calculateFilterVelocity_(unsigned phaseIdx)
Definition: darcyfluxmodule.hh:511
void calculateGradients_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Calculate the gradients which are required to determine the volumetric fluxes.
Definition: darcyfluxmodule.hh:185
short exteriorDofIdx_
Definition: darcyfluxmodule.hh:557
const DimMatrix & intrinsicPermability() const
Returns the intrinsic permeability tensor for a given sub-control volume face.
Definition: darcyfluxmodule.hh:140
DimMatrix K_
Definition: darcyfluxmodule.hh:538
Provides the intensive quantities for the Darcy flux module.
Definition: darcyfluxmodule.hh:90
void update_(const ElementContext &, unsigned, unsigned)
Definition: darcyfluxmodule.hh:93
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:84
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:108
Defines the common properties required by the porous medium multi-phase models.
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:242
This method contains all callback classes for quantities that are required by some extensive quantiti...
Specifies a flux module which uses the Darcy relation.
Definition: darcyfluxmodule.hh:63
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: darcyfluxmodule.hh:71