forchheimerfluxmodule.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
30#ifndef EWOMS_FORCHHEIMER_FLUX_MODULE_HH
31#define EWOMS_FORCHHEIMER_FLUX_MODULE_HH
32
33#include "darcyfluxmodule.hh"
34
35#include <opm/common/Exceptions.hpp>
36
38
39#include <opm/material/common/Valgrind.hpp>
40
41#include <dune/common/fvector.hh>
42#include <dune/common/fmatrix.hh>
43
44#include <cmath>
45
46namespace Opm {
47template <class TypeTag>
48class ForchheimerIntensiveQuantities;
49
50template <class TypeTag>
51class ForchheimerExtensiveQuantities;
52
53template <class TypeTag>
54class ForchheimerBaseProblem;
55
60template <class TypeTag>
62{
66
70 static void registerParameters()
71 {}
72};
73
79template <class TypeTag>
81{
84
85public:
95 template <class Context>
96 Scalar ergunCoefficient(const Context&,
97 unsigned,
98 unsigned) const
99 {
100 throw std::logic_error("Not implemented: Problem::ergunCoefficient()");
101 }
102
114 template <class Context>
115 Evaluation mobilityPassabilityRatio(Context& context,
116 unsigned spaceIdx,
117 unsigned timeIdx,
118 unsigned phaseIdx) const
119 {
120 return 1.0 / context.intensiveQuantities(spaceIdx, timeIdx).fluidState().viscosity(phaseIdx);
121 }
122};
123
128template <class TypeTag>
130{
134
135 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
136
137public:
145 const Evaluation& ergunCoefficient() const
146 { return ergunCoefficient_; }
147
151 const Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const
152 { return mobilityPassabilityRatio_[phaseIdx]; }
153
154protected:
155 void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
156 {
157 const auto& problem = elemCtx.problem();
158 ergunCoefficient_ = problem.ergunCoefficient(elemCtx, dofIdx, timeIdx);
159
160 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
161 mobilityPassabilityRatio_[phaseIdx] =
162 problem.mobilityPassabilityRatio(elemCtx,
163 dofIdx,
164 timeIdx,
165 phaseIdx);
166 }
167
168private:
169 Evaluation ergunCoefficient_;
170 Evaluation mobilityPassabilityRatio_[numPhases];
171};
172
212template <class TypeTag>
214 : public DarcyExtensiveQuantities<TypeTag>
215{
224
225 enum { dimWorld = GridView::dimensionworld };
226 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
227
228 using Toolbox = MathToolbox<Evaluation>;
229
230 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
231 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
232 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
233 using DimEvalMatrix = Dune::FieldMatrix<Evaluation, dimWorld, dimWorld>;
234
235public:
239 const Evaluation& ergunCoefficient() const
240 { return ergunCoefficient_; }
241
248 Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const
249 { return mobilityPassabilityRatio_[phaseIdx]; }
250
251protected:
252 void calculateGradients_(const ElementContext& elemCtx,
253 unsigned faceIdx,
254 unsigned timeIdx)
255 {
256 DarcyExtQuants::calculateGradients_(elemCtx, faceIdx, timeIdx);
257
258 auto focusDofIdx = elemCtx.focusDofIndex();
259 unsigned i = static_cast<unsigned>(this->interiorDofIdx_);
260 unsigned j = static_cast<unsigned>(this->exteriorDofIdx_);
261 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
262 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
263
264 // calculate the square root of the intrinsic permeability
265 assert(isDiagonal_(this->K_));
266 sqrtK_ = 0.0;
267 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
268 sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
269
270 // obtain the Ergun coefficient. Lacking better ideas, we use its the arithmetic mean.
271 if (focusDofIdx == i) {
273 (intQuantsIn.ergunCoefficient() +
274 getValue(intQuantsEx.ergunCoefficient()))/2;
275 }
276 else if (focusDofIdx == j)
278 (getValue(intQuantsIn.ergunCoefficient()) +
279 intQuantsEx.ergunCoefficient())/2;
280 else
282 (getValue(intQuantsIn.ergunCoefficient()) +
283 getValue(intQuantsEx.ergunCoefficient()))/2;
284
285 // obtain the mobility to passability ratio for each phase.
286 for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
287 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
288 continue;
289
290 unsigned upIdx = static_cast<unsigned>(this->upstreamIndex_(phaseIdx));
291 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
292
293 if (focusDofIdx == upIdx) {
294 density_[phaseIdx] =
295 up.fluidState().density(phaseIdx);
296 mobilityPassabilityRatio_[phaseIdx] =
297 up.mobilityPassabilityRatio(phaseIdx);
298 }
299 else {
300 density_[phaseIdx] =
301 getValue(up.fluidState().density(phaseIdx));
302 mobilityPassabilityRatio_[phaseIdx] =
303 getValue(up.mobilityPassabilityRatio(phaseIdx));
304 }
305 }
306 }
307
308 template <class FluidState>
309 void calculateBoundaryGradients_(const ElementContext& elemCtx,
310 unsigned boundaryFaceIdx,
311 unsigned timeIdx,
312 const FluidState& fluidState)
313 {
315 boundaryFaceIdx,
316 timeIdx,
317 fluidState);
318
319 auto focusDofIdx = elemCtx.focusDofIndex();
320 unsigned i = static_cast<unsigned>(this->interiorDofIdx_);
321 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
322
323 // obtain the Ergun coefficient. Because we are on the boundary here, we will
324 // take the Ergun coefficient of the interior
325 if (focusDofIdx == i)
326 ergunCoefficient_ = intQuantsIn.ergunCoefficient();
327 else
328 ergunCoefficient_ = getValue(intQuantsIn.ergunCoefficient());
329
330 // calculate the square root of the intrinsic permeability
331 assert(isDiagonal_(this->K_));
332 sqrtK_ = 0.0;
333 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
334 sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
335
336 for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
337 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
338 continue;
339
340 if (focusDofIdx == i) {
341 density_[phaseIdx] = intQuantsIn.fluidState().density(phaseIdx);
342 mobilityPassabilityRatio_[phaseIdx] = intQuantsIn.mobilityPassabilityRatio(phaseIdx);
343 }
344 else {
345 density_[phaseIdx] =
346 getValue(intQuantsIn.fluidState().density(phaseIdx));
347 mobilityPassabilityRatio_[phaseIdx] =
348 getValue(intQuantsIn.mobilityPassabilityRatio(phaseIdx));
349 }
350 }
351 }
352
359 void calculateFluxes_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
360 {
361 auto focusDofIdx = elemCtx.focusDofIndex();
362 auto i = asImp_().interiorIndex();
363 auto j = asImp_().exteriorIndex();
364 const auto& intQuantsI = elemCtx.intensiveQuantities(i, timeIdx);
365 const auto& intQuantsJ = elemCtx.intensiveQuantities(j, timeIdx);
366
367 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
368 const auto& normal = scvf.normal();
369 Valgrind::CheckDefined(normal);
370
371 // obtain the Ergun coefficient from the intensive quantity object. Until a
372 // better method comes along, we use arithmetic averaging.
373 if (focusDofIdx == i)
375 (intQuantsI.ergunCoefficient() +
376 getValue(intQuantsJ.ergunCoefficient())) / 2;
377 else if (focusDofIdx == j)
379 (getValue(intQuantsI.ergunCoefficient()) +
380 intQuantsJ.ergunCoefficient()) / 2;
381 else
383 (getValue(intQuantsI.ergunCoefficient()) +
384 getValue(intQuantsJ.ergunCoefficient())) / 2;
385
387 // calculate the weights of the upstream and the downstream control volumes
389 for (unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
390 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
391 this->filterVelocity_[phaseIdx] = 0.0;
392 this->volumeFlux_[phaseIdx] = 0.0;
393 continue;
394 }
395
397
398 this->volumeFlux_[phaseIdx] = 0.0;
399 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
400 this->volumeFlux_[phaseIdx] +=
401 this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx];
402 }
403 }
404
408 void calculateBoundaryFluxes_(const ElementContext& elemCtx,
409 unsigned bfIdx,
410 unsigned timeIdx)
411 {
412 const auto& boundaryFace = elemCtx.stencil(timeIdx).boundaryFace(bfIdx);
413 const auto& normal = boundaryFace.normal();
414
416 // calculate the weights of the upstream and the downstream degrees of freedom
418 for (unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
419 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
420 this->filterVelocity_[phaseIdx] = 0.0;
421 this->volumeFlux_[phaseIdx] = 0.0;
422 continue;
423 }
424
426
427 this->volumeFlux_[phaseIdx] = 0.0;
428 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
429 this->volumeFlux_[phaseIdx] +=
430 this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx];
431 }
432 }
433
434 void calculateForchheimerFlux_(unsigned phaseIdx)
435 {
436 // initial guess: filter velocity is zero
437 DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
438 velocity = 0.0;
439
440 // the change of velocity between two consecutive Newton iterations
441 DimEvalVector deltaV(1e5);
442 // the function value that is to be minimized of the equation that is to be
443 // fulfilled
444 DimEvalVector residual;
445 // derivative of equation that is to be solved
446 DimEvalMatrix gradResid;
447
448 // search by means of the Newton method for a root of Forchheimer equation
449 unsigned newtonIter = 0;
450 while (deltaV.one_norm() > 1e-11) {
451 if (newtonIter >= 50)
452 throw NumericalProblem("Could not determine Forchheimer velocity within "
453 + std::to_string(newtonIter)+" iterations");
454 ++newtonIter;
455
456 // calculate the residual and its Jacobian matrix
457 gradForchheimerResid_(residual, gradResid, phaseIdx);
458
459 // newton method
460 gradResid.solve(deltaV, residual);
461 velocity -= deltaV;
462 }
463 }
464
465 void forchheimerResid_(DimEvalVector& residual, unsigned phaseIdx) const
466 {
467 const DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
468
469 // Obtaining the upstreamed quantities
470 const auto& mobility = this->mobility_[phaseIdx];
471 const auto& density = density_[phaseIdx];
473
474 // optain the quantites for the integration point
475 const auto& pGrad = this->potentialGrad_[phaseIdx];
476
477 // residual = v_\alpha
478 residual = velocity;
479
480 // residual += mobility_\alpha K(\grad p_\alpha - \rho_\alpha g)
481 // -> this->K_.usmv(mobility, pGrad, residual);
482 assert(isDiagonal_(this->K_));
483 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
484 residual[dimIdx] += mobility*pGrad[dimIdx]*this->K_[dimIdx][dimIdx];
485
486 // Forchheimer turbulence correction:
487 //
488 // residual +=
489 // \rho_\alpha
490 // * mobility_\alpha
491 // * C_E / \eta_{r,\alpha}
492 // * abs(v_\alpha) * sqrt(K)*v_\alpha
493 //
494 // -> sqrtK_.usmv(density*mobilityPassabilityRatio*ergunCoefficient_*velocity.two_norm(),
495 // velocity,
496 // residual);
497 Evaluation absVel = 0.0;
498 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
499 absVel += velocity[dimIdx]*velocity[dimIdx];
500 // the derivatives of the square root of 0 are undefined, so we must guard
501 // against this case
502 if (absVel <= 0.0)
503 absVel = 0.0;
504 else
505 absVel = Toolbox::sqrt(absVel);
506 const auto& alpha = density*mobilityPassabilityRatio*ergunCoefficient_*absVel;
507 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
508 residual[dimIdx] += sqrtK_[dimIdx]*alpha*velocity[dimIdx];
509 Valgrind::CheckDefined(residual);
510 }
511
512 void gradForchheimerResid_(DimEvalVector& residual,
513 DimEvalMatrix& gradResid,
514 unsigned phaseIdx)
515 {
516 // TODO (?) use AD for this.
517 DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
518 forchheimerResid_(residual, phaseIdx);
519
520 Scalar eps = 1e-11;
521 DimEvalVector tmp;
522 for (unsigned i = 0; i < dimWorld; ++i) {
523 Scalar coordEps = std::max(eps, Toolbox::scalarValue(velocity[i]) * (1 + eps));
524 velocity[i] += coordEps;
525 forchheimerResid_(tmp, phaseIdx);
526 tmp -= residual;
527 tmp /= coordEps;
528 gradResid[i] = tmp;
529 velocity[i] -= coordEps;
530 }
531 }
532
540 bool isDiagonal_(const DimMatrix& K) const
541 {
542 for (unsigned i = 0; i < dimWorld; i++) {
543 for (unsigned j = 0; j < dimWorld; j++) {
544 if (i == j)
545 continue;
546
547 if (std::abs(K[i][j]) > 1e-25)
548 return false;
549 }
550 }
551 return true;
552 }
553
554private:
555 Implementation& asImp_()
556 { return *static_cast<Implementation *>(this); }
557
558 const Implementation& asImp_() const
559 { return *static_cast<const Implementation *>(this); }
560
561protected:
562 // intrinsic permeability tensor and its square root
563 DimVector sqrtK_;
564
565 // Ergun coefficient of all phases at the integration point
567
568 // Passability of all phases at the integration point
569 Evaluation mobilityPassabilityRatio_[numPhases];
570
571 // Density of all phases at the integration point
572 Evaluation density_[numPhases];
573};
574
575} // namespace Opm
576
577#endif
Provides the Darcy flux module.
Definition: darcyfluxmodule.hh:115
Evaluation mobility_[numPhases]
Definition: darcyfluxmodule.hh:539
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:331
short interiorDofIdx_
Definition: darcyfluxmodule.hh:554
short upstreamIndex_(unsigned phaseIdx) const
Definition: darcyfluxmodule.hh:172
EvalDimVector potentialGrad_[numPhases]
Definition: darcyfluxmodule.hh:549
EvalDimVector filterVelocity_[numPhases]
Definition: darcyfluxmodule.hh:542
Evaluation volumeFlux_[numPhases]
Definition: darcyfluxmodule.hh:546
void calculateGradients_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Calculate the gradients which are required to determine the volumetric fluxes.
Definition: darcyfluxmodule.hh:183
short exteriorDofIdx_
Definition: darcyfluxmodule.hh:555
DimMatrix K_
Definition: darcyfluxmodule.hh:536
Provides the defaults for the parameters required by the Forchheimer velocity approach.
Definition: forchheimerfluxmodule.hh:81
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:115
Scalar ergunCoefficient(const Context &, unsigned, unsigned) const
Returns the Ergun coefficient.
Definition: forchheimerfluxmodule.hh:96
Provides the Forchheimer flux module.
Definition: forchheimerfluxmodule.hh:215
void calculateGradients_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Definition: forchheimerfluxmodule.hh:252
void calculateForchheimerFlux_(unsigned phaseIdx)
Definition: forchheimerfluxmodule.hh:434
const Evaluation & ergunCoefficient() const
Return the Ergun coefficent at the face's integration point.
Definition: forchheimerfluxmodule.hh:239
void calculateBoundaryGradients_(const ElementContext &elemCtx, unsigned boundaryFaceIdx, unsigned timeIdx, const FluidState &fluidState)
Definition: forchheimerfluxmodule.hh:309
void forchheimerResid_(DimEvalVector &residual, unsigned phaseIdx) const
Definition: forchheimerfluxmodule.hh:465
Evaluation density_[numPhases]
Definition: forchheimerfluxmodule.hh:572
Evaluation mobilityPassabilityRatio_[numPhases]
Definition: forchheimerfluxmodule.hh:569
DimVector sqrtK_
Definition: forchheimerfluxmodule.hh:563
void gradForchheimerResid_(DimEvalVector &residual, DimEvalMatrix &gradResid, unsigned phaseIdx)
Definition: forchheimerfluxmodule.hh:512
void calculateBoundaryFluxes_(const ElementContext &elemCtx, unsigned bfIdx, unsigned timeIdx)
Calculate the volumetric flux rates of all phases at the domain boundary.
Definition: forchheimerfluxmodule.hh:408
Evaluation ergunCoefficient_
Definition: forchheimerfluxmodule.hh:566
void calculateFluxes_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Calculate the volumetric fluxes of all phases.
Definition: forchheimerfluxmodule.hh:359
bool isDiagonal_(const DimMatrix &K) const
Check whether all off-diagonal entries of a tensor are zero.
Definition: forchheimerfluxmodule.hh:540
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:248
Provides the intensive quantities for the Forchheimer module.
Definition: forchheimerfluxmodule.hh:130
const Evaluation & ergunCoefficient() const
Returns the Ergun coefficient.
Definition: forchheimerfluxmodule.hh:145
void update_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: forchheimerfluxmodule.hh:155
const Evaluation & mobilityPassabilityRatio(unsigned phaseIdx) const
Returns the passability of a phase.
Definition: forchheimerfluxmodule.hh:151
This file contains the necessary classes to calculate the volumetric fluxes out of a pressure potenti...
Declare the properties used by the infrastructure code of the finite volume discretizations.
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:235
Specifies a flux module which uses the Forchheimer relation.
Definition: forchheimerfluxmodule.hh:62
static void registerParameters()
Register all run-time parameters for the flux module.
Definition: forchheimerfluxmodule.hh:70