30 #ifndef EWOMS_DARCY_FLUX_MODULE_HH
31 #define EWOMS_DARCY_FLUX_MODULE_HH
36 #include <opm/common/Valgrind.hpp>
37 #include <opm/common/Unused.hpp>
38 #include <opm/common/ErrorMacros.hpp>
39 #include <opm/common/Exceptions.hpp>
41 #include <dune/common/fvector.hh>
42 #include <dune/common/fmatrix.hh>
47 namespace Properties {
51 template <
class TypeTag>
54 template <
class TypeTag>
57 template <
class TypeTag>
64 template <
class TypeTag>
83 template <
class TypeTag>
84 class DarcyBaseProblem
91 template <
class TypeTag>
92 class DarcyIntensiveQuantities
94 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
96 void update_(const ElementContext& elemCtx OPM_UNUSED,
97 unsigned dofIdx OPM_UNUSED,
98 unsigned timeIdx OPM_UNUSED)
118 template <
class TypeTag>
119 class DarcyExtensiveQuantities
121 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
122 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
123 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
124 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
125 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
126 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
127 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
129 enum { dimWorld = GridView::dimensionworld };
132 typedef typename Opm::MathToolbox<Evaluation> Toolbox;
133 typedef typename FluidSystem::template ParameterCache<Evaluation> ParameterCache;
134 typedef Dune::FieldVector<Evaluation, dimWorld> EvalDimVector;
135 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
136 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
153 {
return potentialGrad_[phaseIdx]; }
162 {
return filterVelocity_[phaseIdx]; }
174 {
return volumeFlux_[phaseIdx]; }
177 short upstreamIndex_(
unsigned phaseIdx)
const
178 {
return upstreamDofIdx_[phaseIdx]; }
180 short downstreamIndex_(
unsigned phaseIdx)
const
181 {
return downstreamDofIdx_[phaseIdx]; }
192 const auto& gradCalc = elemCtx.gradientCalculator();
195 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
196 const auto& faceNormal = scvf.normal();
198 unsigned i = scvf.interiorIndex();
199 unsigned j = scvf.exteriorIndex();
200 interiorDofIdx_ =
static_cast<short>(i);
201 exteriorDofIdx_ =
static_cast<short>(j);
202 unsigned focusDofIdx = elemCtx.focusDofIndex();
205 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
206 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
207 Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
212 gradCalc.calculateGradient(potentialGrad_[phaseIdx],
216 Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
223 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
224 const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
226 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
227 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
229 const auto& posIn = elemCtx.pos(i, timeIdx);
230 const auto& posEx = elemCtx.pos(j, timeIdx);
231 const auto& posFace = scvf.integrationPos();
234 DimVector distVecIn(posIn);
235 DimVector distVecEx(posEx);
236 DimVector distVecTotal(posEx);
238 distVecIn -= posFace;
239 distVecEx -= posFace;
240 distVecTotal -= posIn;
241 Scalar absDistTotalSquared = distVecTotal.two_norm2();
242 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
243 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
249 if (std::is_same<Scalar, Evaluation>::value ||
250 interiorDofIdx_ == static_cast<int>(focusDofIdx))
252 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
253 pStatIn = - rhoIn*(gIn*distVecIn);
256 Scalar rhoIn = Toolbox::value(intQuantsIn.fluidState().density(phaseIdx));
257 pStatIn = - rhoIn*(gIn*distVecIn);
264 if (std::is_same<Scalar, Evaluation>::value ||
265 exteriorDofIdx_ == static_cast<int>(focusDofIdx))
267 const Evaluation& rhoEx = intQuantsEx.fluidState().density(phaseIdx);
268 pStatEx = - rhoEx*(gEx*distVecEx);
271 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
272 pStatEx = - rhoEx*(gEx*distVecEx);
279 Dune::FieldVector<Evaluation, dimWorld> f(distVecTotal);
280 f *= (pStatEx - pStatIn)/absDistTotalSquared;
283 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
284 potentialGrad_[phaseIdx][dimIdx] += f[dimIdx];
286 for (
unsigned dimIdx = 0; dimIdx < potentialGrad_[phaseIdx].size(); ++dimIdx) {
287 if (!std::isfinite(Toolbox::value(potentialGrad_[phaseIdx][dimIdx]))) {
288 OPM_THROW(Opm::NumericalProblem,
289 "Non-finite potential gradient for phase '"
290 << FluidSystem::phaseName(phaseIdx) <<
"'");
296 Opm::Valgrind::SetUndefined(K_);
297 elemCtx.problem().intersectionIntrinsicPermeability(K_, elemCtx, faceIdx, timeIdx);
298 Opm::Valgrind::CheckDefined(K_);
300 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
301 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
302 Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
307 Evaluation tmp = 0.0;
308 for (
unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
309 tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx];
312 upstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
313 downstreamDofIdx_[phaseIdx] = interiorDofIdx_;
316 upstreamDofIdx_[phaseIdx] = interiorDofIdx_;
317 downstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
322 const auto& up = elemCtx.intensiveQuantities(upstreamDofIdx_[phaseIdx], timeIdx);
323 if (upstreamDofIdx_[phaseIdx] == static_cast<int>(focusDofIdx))
324 mobility_[phaseIdx] = up.mobility(phaseIdx);
326 mobility_[phaseIdx] = Toolbox::value(up.mobility(phaseIdx));
336 template <
class Flu
idState>
338 unsigned boundaryFaceIdx,
340 const FluidState& fluidState,
341 const typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache)
343 const auto& gradCalc = elemCtx.gradientCalculator();
347 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
348 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
349 Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
354 gradCalc.calculateBoundaryGradient(potentialGrad_[phaseIdx],
358 Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
361 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
362 auto i = scvf.interiorIndex();
363 interiorDofIdx_ =
static_cast<short>(i);
364 exteriorDofIdx_ = -1;
365 int focusDofIdx = elemCtx.focusDofIndex();
368 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
369 K_ = intQuantsIn.intrinsicPermeability();
375 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
376 const auto& posIn = elemCtx.pos(i, timeIdx);
377 const auto& posFace = scvf.integrationPos();
380 DimVector distVecIn(posIn);
381 distVecIn -= posFace;
382 Scalar absDist = distVecIn.two_norm();
383 Scalar gTimesDist = gIn*distVecIn;
385 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
386 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
390 Evaluation rhoIn = intQuantsIn.fluidState().density(phaseIdx);
391 Evaluation pStatIn = - gTimesDist*rhoIn;
393 Opm::Valgrind::CheckDefined(pStatIn);
401 EvalDimVector f(distVecIn);
402 f *= - pStatIn/absDist;
405 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
406 potentialGrad_[phaseIdx][dimIdx] += f[dimIdx];
408 Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
409 for (
unsigned dimIdx = 0; dimIdx < potentialGrad_[phaseIdx].size(); ++dimIdx) {
410 if (!std::isfinite(Toolbox::value(potentialGrad_[phaseIdx][dimIdx]))) {
411 OPM_THROW(Opm::NumericalProblem,
412 "Non finite potential gradient for phase '"
413 << FluidSystem::phaseName(phaseIdx) <<
"'");
420 const auto& faceNormal = scvf.normal();
422 const auto& matParams = elemCtx.problem().materialLawParams(elemCtx, i, timeIdx);
424 Scalar kr[numPhases];
425 MaterialLaw::relativePermeabilities(kr, matParams, fluidState);
426 Opm::Valgrind::CheckDefined(kr);
428 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
429 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
432 Evaluation tmp = 0.0;
433 for (
unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
434 tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx];
437 upstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
438 downstreamDofIdx_[phaseIdx] = interiorDofIdx_;
441 upstreamDofIdx_[phaseIdx] = interiorDofIdx_;
442 downstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
446 if (upstreamDofIdx_[phaseIdx] < 0)
447 mobility_[phaseIdx] =
448 kr[phaseIdx] / FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
449 else if (upstreamDofIdx_[phaseIdx] != focusDofIdx)
450 mobility_[phaseIdx] = Toolbox::value(intQuantsIn.mobility(phaseIdx));
452 mobility_[phaseIdx] = intQuantsIn.mobility(phaseIdx);
453 Opm::Valgrind::CheckDefined(mobility_[phaseIdx]);
465 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
466 const DimVector& normal = scvf.normal();
467 Opm::Valgrind::CheckDefined(normal);
469 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
470 filterVelocity_[phaseIdx] = 0.0;
471 volumeFlux_[phaseIdx] = 0.0;
472 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
475 asImp_().calculateFilterVelocity_(phaseIdx);
476 Opm::Valgrind::CheckDefined(filterVelocity_[phaseIdx]);
478 volumeFlux_[phaseIdx] = 0.0;
479 for (
unsigned i = 0; i < normal.size(); ++i)
480 volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i];
491 unsigned boundaryFaceIdx,
494 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
495 const DimVector& normal = scvf.normal();
496 Opm::Valgrind::CheckDefined(normal);
498 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
499 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
500 filterVelocity_[phaseIdx] = 0.0;
501 volumeFlux_[phaseIdx] = 0.0;
505 asImp_().calculateFilterVelocity_(phaseIdx);
506 Opm::Valgrind::CheckDefined(filterVelocity_[phaseIdx]);
507 volumeFlux_[phaseIdx] = 0.0;
508 for (
unsigned i = 0; i < normal.size(); ++i)
509 volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i];
513 void calculateFilterVelocity_(
unsigned phaseIdx)
516 assert(std::isfinite(Toolbox::value(mobility_[phaseIdx])));
517 for (
unsigned i = 0; i < K_.M(); ++ i)
518 for (
unsigned j = 0; j < K_.N(); ++ j)
519 assert(std::isfinite(K_[i][j]));
522 K_.mv(potentialGrad_[phaseIdx], filterVelocity_[phaseIdx]);
523 filterVelocity_[phaseIdx] *= - mobility_[phaseIdx];
526 for (
unsigned i = 0; i < filterVelocity_[phaseIdx].size(); ++ i)
527 assert(std::isfinite(Toolbox::value(filterVelocity_[phaseIdx][i])));
532 Implementation& asImp_()
533 {
return *
static_cast<Implementation*
>(
this); }
535 const Implementation& asImp_()
const
536 {
return *
static_cast<const Implementation*
>(
this); }
543 Evaluation mobility_[numPhases];
546 EvalDimVector filterVelocity_[numPhases];
550 Evaluation volumeFlux_[numPhases];
553 EvalDimVector potentialGrad_[numPhases];
556 short upstreamDofIdx_[numPhases];
557 short downstreamDofIdx_[numPhases];
558 short interiorDofIdx_;
559 short exteriorDofIdx_;
void calculateFluxes_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Calculate the volumetric fluxes of all phases.
Definition: darcyfluxmodule.hh:463
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: darcyfluxmodule.hh:74
Specifies a flux module which uses the Darcy relation.
Definition: darcyfluxmodule.hh:65
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:163
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:132
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:161
Provides the intensive quantities for the Darcy flux module.
Definition: darcyfluxmodule.hh:52
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
Defines the common properties required by the porous medium multi-phase models.
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:108
Provides the defaults for the parameters required by the Darcy velocity approach. ...
Definition: darcyfluxmodule.hh:58
Provides the Darcy flux module.
Definition: darcyfluxmodule.hh:55
This method contains all callback classes for quantities that are required by some extensive quantiti...
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:83
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
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:152
void calculateBoundaryGradients_(const ElementContext &elemCtx, unsigned boundaryFaceIdx, unsigned timeIdx, const FluidState &fluidState, const typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache)
Calculate the gradients at the grid boundary which are required to determine the volumetric fluxes...
Definition: darcyfluxmodule.hh:337
const DimMatrix & intrinsicPermability() const
Returns the intrinsic permeability tensor for a given sub-control volume face.
Definition: darcyfluxmodule.hh:143
void calculateBoundaryFluxes_(const ElementContext &elemCtx, unsigned boundaryFaceIdx, unsigned timeIdx)
Calculate the volumetric fluxes at a boundary face of all fluid phases.
Definition: darcyfluxmodule.hh:490
const Evaluation & volumeFlux(unsigned phaseIdx) const
Return the volume flux of a fluid phase at the face's integration point .
Definition: darcyfluxmodule.hh:173
void calculateGradients_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Calculate the gradients which are required to determine the volumetric fluxes.
Definition: darcyfluxmodule.hh:188
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247