28 #ifndef EWOMS_FORCHHEIMER_FLUX_MODULE_HH
29 #define EWOMS_FORCHHEIMER_FLUX_MODULE_HH
35 #include <dune/common/fvector.hh>
36 #include <dune/common/fmatrix.hh>
41 namespace Properties {
47 template <
class TypeTag>
50 template <
class TypeTag>
53 template <
class TypeTag>
60 template <
class TypeTag>
79 template <
class TypeTag>
80 class ForchheimerBaseProblem
82 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
94 template <class Context>
97 OPM_THROW(std::logic_error,
98 "Not implemented: Problem::ergunCoefficient()");
112 template <
class Context>
116 return 1.0 / context.intensiveQuantities(spaceIdx, timeIdx).fluidState().viscosity(
125 template <
class TypeTag>
126 class ForchheimerIntensiveQuantities
128 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
129 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
143 return ergunCoefficient_;
151 return mobilityPassabilityRatio_[phaseIdx];
155 void update_(
const ElementContext &elemCtx,
int dofIdx,
int timeIdx)
157 const auto &problem = elemCtx.problem();
158 ergunCoefficient_ = problem.ergunCoefficient(elemCtx, dofIdx, timeIdx);
160 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
161 mobilityPassabilityRatio_[phaseIdx] =
162 problem.mobilityPassabilityRatio(elemCtx,
169 Scalar ergunCoefficient_;
170 Scalar mobilityPassabilityRatio_[numPhases];
215 template <
class TypeTag>
216 class ForchheimerExtensiveQuantities
217 :
public DarcyExtensiveQuantities<TypeTag>
219 typedef DarcyExtensiveQuantities<TypeTag> DarcyExtQuants;
220 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
221 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
222 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
223 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
224 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
225 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
228 enum { dimWorld = GridView::dimensionworld };
231 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
232 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
257 const auto &intQuantsIn = elemCtx.intensiveQuantities(this->
interiorDofIdx_, timeIdx);
258 const auto &intQuantsEx = elemCtx.intensiveQuantities(this->
exteriorDofIdx_, timeIdx);
263 for (
int i = 0; i < dimWorld; ++i)
264 sqrtK_[i][i] = std::sqrt(this->
K_[i][i]);
267 ergunCoefficient_ = (intQuantsIn.ergunCoefficient() + intQuantsEx.ergunCoefficient())/2;
270 for (
int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
271 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
274 const auto &up = elemCtx.intensiveQuantities(this->
upstreamIndex_(phaseIdx), timeIdx);
276 density_[phaseIdx] = up.fluidState().density(phaseIdx);
281 template <
class Flu
idState>
285 const FluidState& fluidState,
286 const typename FluidSystem::ParameterCache& paramCache)
294 const auto &intQuantsIn = elemCtx.intensiveQuantities(this->
interiorDofIdx_, timeIdx);
303 for (
int i = 0; i < dimWorld; ++i)
304 sqrtK_[i][i] = std::sqrt(this->
K_[i][i]);
306 for (
int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
307 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
310 density_[phaseIdx] = intQuantsIn.fluidState().density(phaseIdx);
323 const auto &intQuantsI = elemCtx.intensiveQuantities(asImp_().interiorIndex(), timeIdx);
324 const auto &intQuantsJ = elemCtx.intensiveQuantities(asImp_().exteriorIndex(), timeIdx);
326 const auto &scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
327 const auto &normal = scvf.normal();
328 Valgrind::CheckDefined(normal);
332 ergunCoefficient_ = (intQuantsI.ergunCoefficient() + intQuantsJ.ergunCoefficient()) / 2;
337 for (
int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
338 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
356 const auto &boundaryFace = elemCtx.stencil(timeIdx).boundaryFace(bfIdx);
357 const auto &normal = boundaryFace.normal();
362 for (
int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
363 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
381 DimVector deltaV(1e5);
390 while (deltaV.two_norm() > 1e-11) {
391 if (newtonIter >= 50)
392 OPM_THROW(Opm::NumericalProblem,
393 "Could not determine Forchheimer velocity within "
394 << newtonIter <<
" iterations");
401 gradResid.solve(deltaV, residual);
411 Scalar mobility = this->
mobility_[phaseIdx];
412 Scalar density =
density_[phaseIdx];
422 this->
K_.usmv(mobility, pGrad, residual);
435 Valgrind::CheckDefined(residual);
446 for (
int i = 0; i < dimWorld; ++i) {
447 Scalar coordEps = std::max(eps, velocity[i] * (1 + eps));
448 velocity[i] += coordEps;
453 velocity[i] -= coordEps;
466 for (
int i = 0; i < dimWorld; i++) {
467 for (
int j = 0; j < dimWorld; j++) {
471 if (std::abs(K[i][j]) > 1e-25)
479 Implementation &asImp_()
480 {
return *
static_cast<Implementation *
>(
this); }
482 const Implementation &asImp_()
const
483 {
return *
static_cast<const Implementation *
>(
this); }
Specifies a flux module which uses the Forchheimer relation.
Definition: forchheimerfluxmodule.hh:61
Scalar ergunCoefficient(Context &context, int spaceIdx, int timeIdx) const
Returns the Ergun coefficient.
Definition: forchheimerfluxmodule.hh:95
void calculateForchheimerFlux_(int phaseIdx)
Definition: forchheimerfluxmodule.hh:374
short interiorDofIdx_
Definition: darcyfluxmodule.hh:521
Scalar mobilityPassabilityRatio(int phaseIdx) const
Returns the passability of a phase.
Definition: forchheimerfluxmodule.hh:149
void forchheimerResid_(DimVector &residual, int phaseIdx) const
Definition: forchheimerfluxmodule.hh:406
Declare the properties used by the infrastructure code of the finite volume discretizations.
DimMatrix sqrtK_
Definition: forchheimerfluxmodule.hh:487
Scalar ergunCoefficient() const
Returns the Ergun coefficient.
Definition: forchheimerfluxmodule.hh:141
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: forchheimerfluxmodule.hh:70
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
ForchheimerBaseProblem< TypeTag > FluxBaseProblem
Definition: forchheimerfluxmodule.hh:65
ForchheimerIntensiveQuantities< TypeTag > FluxIntensiveQuantities
Definition: forchheimerfluxmodule.hh:63
Scalar density_[numPhases]
Definition: forchheimerfluxmodule.hh:496
Provides the intensive quantities for the Forchheimer module.
Definition: forchheimerfluxmodule.hh:48
void calculateBoundaryGradients_(const ElementContext &elemCtx, int boundaryFaceIdx, int timeIdx, const FluidState &fluidState, const typename FluidSystem::ParameterCache ¶mCache)
Definition: forchheimerfluxmodule.hh:282
void calculateFluxes_(const ElementContext &elemCtx, int scvfIdx, int timeIdx)
Calculate the volumetric fluxes of all phases.
Definition: forchheimerfluxmodule.hh:321
NEW_PROP_TAG(Grid)
The type of the DUNE grid.
short upstreamIndex_(int phaseIdx) const
Definition: darcyfluxmodule.hh:182
Provides the Forchheimer flux module.
Definition: forchheimerfluxmodule.hh:51
Scalar ergunCoefficient_
Definition: forchheimerfluxmodule.hh:490
void gradForchheimerResid_(DimVector &residual, DimMatrix &gradResid, int phaseIdx)
Definition: forchheimerfluxmodule.hh:438
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
Scalar mobilityPassabilityRatio(int phaseIdx) const
Return the ratio of the mobility divided by the passability at the face's integration point for a giv...
Definition: forchheimerfluxmodule.hh:247
Scalar mobilityPassabilityRatio(Context &context, int spaceIdx, int timeIdx, int phaseIdx) const
Returns the ratio between the phase mobility and its passability for a given fluid phase ...
Definition: forchheimerfluxmodule.hh:113
Evaluation volumeFlux_[numPhases]
Definition: darcyfluxmodule.hh:534
This file contains the necessary classes to calculate the volumetric fluxes out of a pressure potenti...
DimMatrix K_
Definition: darcyfluxmodule.hh:518
Scalar ergunCoefficient() const
Return the Ergun coefficent at the face's integration point.
Definition: forchheimerfluxmodule.hh:238
bool isDiagonal_(const DimMatrix &K) const
Check whether all off-diagonal entries of a tensor are zero.
Definition: forchheimerfluxmodule.hh:464
void calculateBoundaryFluxes_(const ElementContext &elemCtx, int bfIdx, int timeIdx)
Calculate the volumetric flux rates of all phases at the domain boundary.
Definition: forchheimerfluxmodule.hh:352
EvalDimVector filterVelocity_[numPhases]
Definition: darcyfluxmodule.hh:530
Provides the defaults for the parameters required by the Forchheimer velocity approach.
Definition: forchheimerfluxmodule.hh:54
void calculateGradients_(const ElementContext &elemCtx, int faceIdx, int timeIdx)
Definition: forchheimerfluxmodule.hh:251
void calculateGradients_(const ElementContext &elemCtx, int faceIdx, int timeIdx)
Calculate the gradients which are required to determine the volumetric fluxes.
Definition: darcyfluxmodule.hh:193
EvalDimVector potentialGrad_[numPhases]
Definition: darcyfluxmodule.hh:537
Scalar mobilityPassabilityRatio_[numPhases]
Definition: forchheimerfluxmodule.hh:493
void update_(const ElementContext &elemCtx, int dofIdx, int timeIdx)
Definition: forchheimerfluxmodule.hh:155
ForchheimerExtensiveQuantities< TypeTag > FluxExtensiveQuantities
Definition: forchheimerfluxmodule.hh:64
short exteriorDofIdx_
Definition: darcyfluxmodule.hh:522