30 #ifndef EWOMS_FORCHHEIMER_FLUX_MODULE_HH 31 #define EWOMS_FORCHHEIMER_FLUX_MODULE_HH 37 #include <opm/common/Valgrind.hpp> 38 #include <opm/common/Unused.hpp> 39 #include <opm/common/ErrorMacros.hpp> 40 #include <opm/common/Exceptions.hpp> 42 #include <dune/common/fvector.hh> 43 #include <dune/common/fmatrix.hh> 48 namespace Properties {
53 template <
class TypeTag>
56 template <
class TypeTag>
59 template <
class TypeTag>
66 template <
class TypeTag>
85 template <
class TypeTag>
86 class ForchheimerBaseProblem
88 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
89 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
101 template <
class Context>
103 unsigned spaceIdx OPM_UNUSED,
104 unsigned timeIdx OPM_UNUSED)
const 106 OPM_THROW(std::logic_error,
107 "Not implemented: Problem::ergunCoefficient()");
121 template <
class Context>
125 unsigned phaseIdx)
const 127 return 1.0 / context.intensiveQuantities(spaceIdx, timeIdx).fluidState().viscosity(phaseIdx);
135 template <
class TypeTag>
136 class ForchheimerIntensiveQuantities
138 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
139 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
140 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
153 {
return ergunCoefficient_; }
159 {
return mobilityPassabilityRatio_[phaseIdx]; }
162 void update_(
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
164 const auto& problem = elemCtx.problem();
165 ergunCoefficient_ = problem.ergunCoefficient(elemCtx, dofIdx, timeIdx);
167 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
168 mobilityPassabilityRatio_[phaseIdx] =
169 problem.mobilityPassabilityRatio(elemCtx,
176 Evaluation ergunCoefficient_;
177 Evaluation mobilityPassabilityRatio_[numPhases];
219 template <
class TypeTag>
220 class ForchheimerExtensiveQuantities
221 :
public DarcyExtensiveQuantities<TypeTag>
223 typedef DarcyExtensiveQuantities<TypeTag> DarcyExtQuants;
224 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
225 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
226 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
227 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
228 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
229 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
230 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
232 enum { dimWorld = GridView::dimensionworld };
235 typedef Opm::MathToolbox<Evaluation> Toolbox;
237 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
238 typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
239 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
240 typedef Dune::FieldMatrix<Evaluation, dimWorld, dimWorld> DimEvalMatrix;
247 {
return ergunCoefficient_; }
256 {
return mobilityPassabilityRatio_[phaseIdx]; }
259 void calculateGradients_(
const ElementContext& elemCtx,
265 auto focusDofIdx = elemCtx.focusDofIndex();
266 unsigned i =
static_cast<unsigned>(this->interiorDofIdx_);
267 unsigned j =
static_cast<unsigned>(this->exteriorDofIdx_);
268 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
269 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
274 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
275 sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
278 if (focusDofIdx == i) {
280 (intQuantsIn.ergunCoefficient() +
281 Toolbox::value(intQuantsEx.ergunCoefficient()))/2;
283 else if (focusDofIdx == j)
285 (Toolbox::value(intQuantsIn.ergunCoefficient()) +
286 intQuantsEx.ergunCoefficient())/2;
289 (Toolbox::value(intQuantsIn.ergunCoefficient()) +
290 Toolbox::value(intQuantsEx.ergunCoefficient()))/2;
293 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
294 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
297 unsigned upIdx =
static_cast<unsigned>(this->upstreamIndex_(phaseIdx));
298 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
300 if (focusDofIdx == upIdx) {
302 up.fluidState().density(phaseIdx);
303 mobilityPassabilityRatio_[phaseIdx] =
304 up.mobilityPassabilityRatio(phaseIdx);
308 Toolbox::value(up.fluidState().density(phaseIdx));
309 mobilityPassabilityRatio_[phaseIdx] =
310 Toolbox::value(up.mobilityPassabilityRatio(phaseIdx));
315 template <
class Flu
idState>
316 void calculateBoundaryGradients_(
const ElementContext& elemCtx,
317 unsigned boundaryFaceIdx,
319 const FluidState& fluidState,
320 const typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache)
328 auto focusDofIdx = elemCtx.focusDofIndex();
329 unsigned i =
static_cast<unsigned>(this->interiorDofIdx_);
330 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
334 if (focusDofIdx == i)
335 ergunCoefficient_ = intQuantsIn.ergunCoefficient();
337 ergunCoefficient_ = Toolbox::value(intQuantsIn.ergunCoefficient());
342 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
343 sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
345 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
346 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
349 if (focusDofIdx == i) {
350 density_[phaseIdx] = intQuantsIn.fluidState().density(phaseIdx);
351 mobilityPassabilityRatio_[phaseIdx] = intQuantsIn.mobilityPassabilityRatio(phaseIdx);
355 Toolbox::value(intQuantsIn.fluidState().density(phaseIdx));
356 mobilityPassabilityRatio_[phaseIdx] =
357 Toolbox::value(intQuantsIn.mobilityPassabilityRatio(phaseIdx));
370 auto focusDofIdx = elemCtx.focusDofIndex();
371 auto i = asImp_().interiorIndex();
372 auto j = asImp_().exteriorIndex();
373 const auto& intQuantsI = elemCtx.intensiveQuantities(i, timeIdx);
374 const auto& intQuantsJ = elemCtx.intensiveQuantities(j, timeIdx);
376 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
377 const auto& normal = scvf.normal();
378 Opm::Valgrind::CheckDefined(normal);
382 if (focusDofIdx == i)
384 (intQuantsI.ergunCoefficient() +
385 Toolbox::value(intQuantsJ.ergunCoefficient())) / 2;
386 else if (focusDofIdx == j)
388 (Toolbox::value(intQuantsI.ergunCoefficient()) +
389 intQuantsJ.ergunCoefficient()) / 2;
392 (Toolbox::value(intQuantsI.ergunCoefficient()) +
393 Toolbox::value(intQuantsJ.ergunCoefficient())) / 2;
398 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
399 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
400 this->filterVelocity_[phaseIdx] = 0.0;
401 this->volumeFlux_[phaseIdx] = 0.0;
405 calculateForchheimerFlux_(phaseIdx);
407 this->volumeFlux_[phaseIdx] = 0.0;
408 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
409 this->volumeFlux_[phaseIdx] +=
410 this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx];
421 const auto& boundaryFace = elemCtx.stencil(timeIdx).boundaryFace(bfIdx);
422 const auto& normal = boundaryFace.normal();
427 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
428 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
429 this->filterVelocity_[phaseIdx] = 0.0;
430 this->volumeFlux_[phaseIdx] = 0.0;
434 calculateForchheimerFlux_(phaseIdx);
436 this->volumeFlux_[phaseIdx] = 0.0;
437 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
438 this->volumeFlux_[phaseIdx] +=
439 this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx];
443 void calculateForchheimerFlux_(
unsigned phaseIdx)
446 DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
450 DimEvalVector deltaV(1e5);
453 DimEvalVector residual;
455 DimEvalMatrix gradResid;
458 unsigned newtonIter = 0;
459 while (deltaV.one_norm() > 1e-11) {
460 if (newtonIter >= 50)
461 OPM_THROW(Opm::NumericalProblem,
462 "Could not determine Forchheimer velocity within " 463 << newtonIter <<
" iterations");
467 gradForchheimerResid_(residual, gradResid, phaseIdx);
470 gradResid.solve(deltaV, residual);
475 void forchheimerResid_(DimEvalVector& residual,
unsigned phaseIdx)
const 477 const DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
480 const auto& mobility = this->mobility_[phaseIdx];
481 const auto& density = density_[phaseIdx];
485 const auto& pGrad = this->potentialGrad_[phaseIdx];
493 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
494 residual[dimIdx] += mobility*pGrad[dimIdx]*this->K_[dimIdx][dimIdx];
507 Evaluation absVel = 0.0;
508 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
509 absVel += velocity[dimIdx]*velocity[dimIdx];
515 absVel = Toolbox::sqrt(absVel);
517 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
518 residual[dimIdx] += sqrtK_[dimIdx]*alpha*velocity[dimIdx];
519 Opm::Valgrind::CheckDefined(residual);
522 void gradForchheimerResid_(DimEvalVector& residual,
523 DimEvalMatrix& gradResid,
527 DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
528 forchheimerResid_(residual, phaseIdx);
532 for (
unsigned i = 0; i < dimWorld; ++i) {
533 Scalar coordEps = std::max(eps, Toolbox::scalarValue(velocity[i]) * (1 + eps));
534 velocity[i] += coordEps;
535 forchheimerResid_(tmp, phaseIdx);
539 velocity[i] -= coordEps;
552 for (
unsigned i = 0; i < dimWorld; i++) {
553 for (
unsigned j = 0; j < dimWorld; j++) {
557 if (std::abs(K[i][j]) > 1e-25)
565 Implementation& asImp_()
566 {
return *
static_cast<Implementation *
>(
this); }
568 const Implementation& asImp_()
const 569 {
return *
static_cast<const Implementation *
>(
this); }
576 Evaluation ergunCoefficient_;
579 Evaluation mobilityPassabilityRatio_[numPhases];
582 Evaluation density_[numPhases];
Provides the intensive quantities for the Forchheimer module.
Definition: forchheimerfluxmodule.hh:54
Scalar ergunCoefficient(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Returns the Ergun coefficient.
Definition: forchheimerfluxmodule.hh:102
Provides the defaults for the parameters required by the Forchheimer velocity approach.
Definition: forchheimerfluxmodule.hh:60
Evaluation & mobilityPassabilityRatio(unsigned phaseIdx) const
Return the ratio of the mobility divided by the passability at the face's integration point for a giv...
Definition: forchheimerfluxmodule.hh:255
Definition: baseauxiliarymodule.hh:37
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const Evaluation & ergunCoefficient() const
Return the Ergun coefficent at the face's integration point.
Definition: forchheimerfluxmodule.hh:246
Evaluation mobilityPassabilityRatio(Context &context, unsigned spaceIdx, unsigned timeIdx, unsigned phaseIdx) const
Returns the ratio between the phase mobility and its passability for a given fluid phase ...
Definition: forchheimerfluxmodule.hh:122
const Evaluation & mobilityPassabilityRatio(unsigned phaseIdx) const
Returns the passability of a phase.
Definition: forchheimerfluxmodule.hh:158
Declare the properties used by the infrastructure code of the finite volume discretizations.
const Evaluation & ergunCoefficient() const
Returns the Ergun coefficient.
Definition: forchheimerfluxmodule.hh:152
This file contains the necessary classes to calculate the volumetric fluxes out of a pressure potenti...
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: forchheimerfluxmodule.hh:76
void calculateFluxes_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Calculate the volumetric fluxes of all phases.
Definition: forchheimerfluxmodule.hh:368
Specifies a flux module which uses the Forchheimer relation.
Definition: forchheimerfluxmodule.hh:67
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
bool isDiagonal_(const DimMatrix &K) const
Check whether all off-diagonal entries of a tensor are zero.
Definition: forchheimerfluxmodule.hh:550
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
void calculateBoundaryFluxes_(const ElementContext &elemCtx, unsigned bfIdx, unsigned timeIdx)
Calculate the volumetric flux rates of all phases at the domain boundary.
Definition: forchheimerfluxmodule.hh:417
Provides the Forchheimer flux module
Definition: forchheimerfluxmodule.hh:57