All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
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 
36 
37 #include <opm/common/Valgrind.hpp>
38 #include <opm/common/Unused.hpp>
39 #include <opm/common/ErrorMacros.hpp>
40 #include <opm/common/Exceptions.hpp>
41 
42 #include <dune/common/fvector.hh>
43 #include <dune/common/fmatrix.hh>
44 
45 #include <cmath>
46 
47 namespace Ewoms {
48 namespace Properties {
49 NEW_PROP_TAG(MaterialLaw);
50 }}
51 
52 namespace Ewoms {
53 template <class TypeTag>
55 
56 template <class TypeTag>
58 
59 template <class TypeTag>
61 
66 template <class TypeTag>
68 {
72 
76  static void registerParameters()
77  {}
78 };
79 
85 template <class TypeTag>
86 class ForchheimerBaseProblem
87 {
88  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
89  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
90 
91 public:
101  template <class Context>
102  Scalar ergunCoefficient(const Context& context OPM_UNUSED,
103  unsigned spaceIdx OPM_UNUSED,
104  unsigned timeIdx OPM_UNUSED) const
105  {
106  OPM_THROW(std::logic_error,
107  "Not implemented: Problem::ergunCoefficient()");
108  }
109 
121  template <class Context>
122  Evaluation mobilityPassabilityRatio(Context& context,
123  unsigned spaceIdx,
124  unsigned timeIdx,
125  unsigned phaseIdx) const
126  {
127  return 1.0 / context.intensiveQuantities(spaceIdx, timeIdx).fluidState().viscosity(phaseIdx);
128  }
129 };
130 
135 template <class TypeTag>
136 class ForchheimerIntensiveQuantities
137 {
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;
141 
142  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
143 
144 public:
152  const Evaluation& ergunCoefficient() const
153  { return ergunCoefficient_; }
154 
158  const Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const
159  { return mobilityPassabilityRatio_[phaseIdx]; }
160 
161 protected:
162  void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
163  {
164  const auto& problem = elemCtx.problem();
165  ergunCoefficient_ = problem.ergunCoefficient(elemCtx, dofIdx, timeIdx);
166 
167  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
168  mobilityPassabilityRatio_[phaseIdx] =
169  problem.mobilityPassabilityRatio(elemCtx,
170  dofIdx,
171  timeIdx,
172  phaseIdx);
173  }
174 
175 private:
176  Evaluation ergunCoefficient_;
177  Evaluation mobilityPassabilityRatio_[numPhases];
178 };
179 
219 template <class TypeTag>
220 class ForchheimerExtensiveQuantities
221  : public DarcyExtensiveQuantities<TypeTag>
222 {
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;
231 
232  enum { dimWorld = GridView::dimensionworld };
233  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
234 
235  typedef Opm::MathToolbox<Evaluation> Toolbox;
236 
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;
241 
242 public:
246  const Evaluation& ergunCoefficient() const
247  { return ergunCoefficient_; }
248 
255  Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const
256  { return mobilityPassabilityRatio_[phaseIdx]; }
257 
258 protected:
259  void calculateGradients_(const ElementContext& elemCtx,
260  unsigned faceIdx,
261  unsigned timeIdx)
262  {
263  DarcyExtQuants::calculateGradients_(elemCtx, faceIdx, timeIdx);
264 
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);
270 
271  // calculate the square root of the intrinsic permeability
272  assert(isDiagonal_(this->K_));
273  sqrtK_ = 0.0;
274  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
275  sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
276 
277  // obtain the Ergun coefficient. Lacking better ideas, we use its the arithmetic mean.
278  if (focusDofIdx == i) {
279  ergunCoefficient_ =
280  (intQuantsIn.ergunCoefficient() +
281  Toolbox::value(intQuantsEx.ergunCoefficient()))/2;
282  }
283  else if (focusDofIdx == j)
284  ergunCoefficient_ =
285  (Toolbox::value(intQuantsIn.ergunCoefficient()) +
286  intQuantsEx.ergunCoefficient())/2;
287  else
288  ergunCoefficient_ =
289  (Toolbox::value(intQuantsIn.ergunCoefficient()) +
290  Toolbox::value(intQuantsEx.ergunCoefficient()))/2;
291 
292  // obtain the mobility to passability ratio for each phase.
293  for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
294  if (!elemCtx.model().phaseIsConsidered(phaseIdx))
295  continue;
296 
297  unsigned upIdx = static_cast<unsigned>(this->upstreamIndex_(phaseIdx));
298  const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
299 
300  if (focusDofIdx == upIdx) {
301  density_[phaseIdx] =
302  up.fluidState().density(phaseIdx);
303  mobilityPassabilityRatio_[phaseIdx] =
304  up.mobilityPassabilityRatio(phaseIdx);
305  }
306  else {
307  density_[phaseIdx] =
308  Toolbox::value(up.fluidState().density(phaseIdx));
309  mobilityPassabilityRatio_[phaseIdx] =
310  Toolbox::value(up.mobilityPassabilityRatio(phaseIdx));
311  }
312  }
313  }
314 
315  template <class FluidState>
316  void calculateBoundaryGradients_(const ElementContext& elemCtx,
317  unsigned boundaryFaceIdx,
318  unsigned timeIdx,
319  const FluidState& fluidState,
320  const typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache)
321  {
323  boundaryFaceIdx,
324  timeIdx,
325  fluidState,
326  paramCache);
327 
328  auto focusDofIdx = elemCtx.focusDofIndex();
329  unsigned i = static_cast<unsigned>(this->interiorDofIdx_);
330  const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
331 
332  // obtain the Ergun coefficient. Because we are on the boundary here, we will
333  // take the Ergun coefficient of the interior
334  if (focusDofIdx == i)
335  ergunCoefficient_ = intQuantsIn.ergunCoefficient();
336  else
337  ergunCoefficient_ = Toolbox::value(intQuantsIn.ergunCoefficient());
338 
339  // calculate the square root of the intrinsic permeability
340  assert(isDiagonal_(this->K_));
341  sqrtK_ = 0.0;
342  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
343  sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
344 
345  for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
346  if (!elemCtx.model().phaseIsConsidered(phaseIdx))
347  continue;
348 
349  if (focusDofIdx == i) {
350  density_[phaseIdx] = intQuantsIn.fluidState().density(phaseIdx);
351  mobilityPassabilityRatio_[phaseIdx] = intQuantsIn.mobilityPassabilityRatio(phaseIdx);
352  }
353  else {
354  density_[phaseIdx] =
355  Toolbox::value(intQuantsIn.fluidState().density(phaseIdx));
356  mobilityPassabilityRatio_[phaseIdx] =
357  Toolbox::value(intQuantsIn.mobilityPassabilityRatio(phaseIdx));
358  }
359  }
360  }
361 
368  void calculateFluxes_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
369  {
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);
375 
376  const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
377  const auto& normal = scvf.normal();
378  Opm::Valgrind::CheckDefined(normal);
379 
380  // obtain the Ergun coefficient from the intensive quantity object. Until a
381  // better method comes along, we use arithmetic averaging.
382  if (focusDofIdx == i)
383  ergunCoefficient_ =
384  (intQuantsI.ergunCoefficient() +
385  Toolbox::value(intQuantsJ.ergunCoefficient())) / 2;
386  else if (focusDofIdx == j)
387  ergunCoefficient_ =
388  (Toolbox::value(intQuantsI.ergunCoefficient()) +
389  intQuantsJ.ergunCoefficient()) / 2;
390  else
391  ergunCoefficient_ =
392  (Toolbox::value(intQuantsI.ergunCoefficient()) +
393  Toolbox::value(intQuantsJ.ergunCoefficient())) / 2;
394 
396  // calculate the weights of the upstream and the downstream control volumes
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;
402  continue;
403  }
404 
405  calculateForchheimerFlux_(phaseIdx);
406 
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];
411  }
412  }
413 
417  void calculateBoundaryFluxes_(const ElementContext& elemCtx,
418  unsigned bfIdx,
419  unsigned timeIdx)
420  {
421  const auto& boundaryFace = elemCtx.stencil(timeIdx).boundaryFace(bfIdx);
422  const auto& normal = boundaryFace.normal();
423 
425  // calculate the weights of the upstream and the downstream degrees of freedom
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;
431  continue;
432  }
433 
434  calculateForchheimerFlux_(phaseIdx);
435 
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];
440  }
441  }
442 
443  void calculateForchheimerFlux_(unsigned phaseIdx)
444  {
445  // initial guess: filter velocity is zero
446  DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
447  velocity = 0.0;
448 
449  // the change of velocity between two consecutive Newton iterations
450  DimEvalVector deltaV(1e5);
451  // the function value that is to be minimized of the equation that is to be
452  // fulfilled
453  DimEvalVector residual;
454  // derivative of equation that is to be solved
455  DimEvalMatrix gradResid;
456 
457  // search by means of the Newton method for a root of Forchheimer equation
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");
464  ++newtonIter;
465 
466  // calculate the residual and its Jacobian matrix
467  gradForchheimerResid_(residual, gradResid, phaseIdx);
468 
469  // newton method
470  gradResid.solve(deltaV, residual);
471  velocity -= deltaV;
472  }
473  }
474 
475  void forchheimerResid_(DimEvalVector& residual, unsigned phaseIdx) const
476  {
477  const DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
478 
479  // Obtaining the upstreamed quantities
480  const auto& mobility = this->mobility_[phaseIdx];
481  const auto& density = density_[phaseIdx];
482  const auto& mobilityPassabilityRatio = mobilityPassabilityRatio_[phaseIdx];
483 
484  // optain the quantites for the integration point
485  const auto& pGrad = this->potentialGrad_[phaseIdx];
486 
487  // residual = v_\alpha
488  residual = velocity;
489 
490  // residual += mobility_\alpha K(\grad p_\alpha - \rho_\alpha g)
491  // -> this->K_.usmv(mobility, pGrad, residual);
492  assert(isDiagonal_(this->K_));
493  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
494  residual[dimIdx] += mobility*pGrad[dimIdx]*this->K_[dimIdx][dimIdx];
495 
496  // Forchheimer turbulence correction:
497  //
498  // residual +=
499  // \rho_\alpha
500  // * mobility_\alpha
501  // * C_E / \eta_{r,\alpha}
502  // * abs(v_\alpha) * sqrt(K)*v_\alpha
503  //
504  // -> sqrtK_.usmv(density*mobilityPassabilityRatio*ergunCoefficient_*velocity.two_norm(),
505  // velocity,
506  // residual);
507  Evaluation absVel = 0.0;
508  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
509  absVel += velocity[dimIdx]*velocity[dimIdx];
510  // the derivatives of the square root of 0 are undefined, so we must guard
511  // against this case
512  if (absVel <= 0.0)
513  absVel = 0.0;
514  else
515  absVel = Toolbox::sqrt(absVel);
516  const auto& alpha = density*mobilityPassabilityRatio*ergunCoefficient_*absVel;
517  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
518  residual[dimIdx] += sqrtK_[dimIdx]*alpha*velocity[dimIdx];
519  Opm::Valgrind::CheckDefined(residual);
520  }
521 
522  void gradForchheimerResid_(DimEvalVector& residual,
523  DimEvalMatrix& gradResid,
524  unsigned phaseIdx)
525  {
526  // TODO (?) use AD for this.
527  DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
528  forchheimerResid_(residual, phaseIdx);
529 
530  Scalar eps = 1e-11;
531  DimEvalVector tmp;
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);
536  tmp -= residual;
537  tmp /= coordEps;
538  gradResid[i] = tmp;
539  velocity[i] -= coordEps;
540  }
541  }
542 
550  bool isDiagonal_(const DimMatrix& K) const
551  {
552  for (unsigned i = 0; i < dimWorld; i++) {
553  for (unsigned j = 0; j < dimWorld; j++) {
554  if (i == j)
555  continue;
556 
557  if (std::abs(K[i][j]) > 1e-25)
558  return false;
559  }
560  }
561  return true;
562  }
563 
564 private:
565  Implementation& asImp_()
566  { return *static_cast<Implementation *>(this); }
567 
568  const Implementation& asImp_() const
569  { return *static_cast<const Implementation *>(this); }
570 
571 protected:
572  // intrinsic permeability tensor and its square root
573  DimVector sqrtK_;
574 
575  // Ergun coefficient of all phases at the integration point
576  Evaluation ergunCoefficient_;
577 
578  // Passability of all phases at the integration point
579  Evaluation mobilityPassabilityRatio_[numPhases];
580 
581  // Density of all phases at the integration point
582  Evaluation density_[numPhases];
583 };
584 
585 } // namespace Ewoms
586 
587 #endif
Provides the intensive quantities for the Forchheimer module.
Definition: forchheimerfluxmodule.hh:54
Provides the defaults for the parameters required by the Forchheimer velocity approach.
Definition: forchheimerfluxmodule.hh:60
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
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
Evaluation & mobilityPassabilityRatio(unsigned phaseIdx) const
Return the ratio of the mobility divided by the passability at the face&#39;s integration point for a giv...
Definition: forchheimerfluxmodule.hh:255
Scalar ergunCoefficient(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Returns the Ergun coefficient.
Definition: forchheimerfluxmodule.hh:102
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
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 & ergunCoefficient() const
Return the Ergun coefficent at the face&#39;s integration point.
Definition: forchheimerfluxmodule.hh:246
void calculateBoundaryGradients_(const ElementContext &elemCtx, unsigned boundaryFaceIdx, unsigned timeIdx, const FluidState &fluidState, const typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache)
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
const Evaluation & mobilityPassabilityRatio(unsigned phaseIdx) const
Returns the passability of a phase.
Definition: forchheimerfluxmodule.hh:158
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