All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
darcyfluxmodule.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_DARCY_FLUX_MODULE_HH
31 #define EWOMS_DARCY_FLUX_MODULE_HH
32 
35 
36 #include <opm/common/Valgrind.hpp>
37 #include <opm/common/Unused.hpp>
38 #include <opm/common/ErrorMacros.hpp>
39 #include <opm/common/Exceptions.hpp>
40 
41 #include <dune/common/fvector.hh>
42 #include <dune/common/fmatrix.hh>
43 
44 #include <cmath>
45 
46 namespace Ewoms {
47 namespace Properties {
48 NEW_PROP_TAG(MaterialLaw);
49 }
50 
51 template <class TypeTag>
53 
54 template <class TypeTag>
56 
57 template <class TypeTag>
59 
64 template <class TypeTag>
66 {
70 
74  static void registerParameters()
75  { }
76 };
77 
83 template <class TypeTag>
84 class DarcyBaseProblem
85 { };
86 
91 template <class TypeTag>
92 class DarcyIntensiveQuantities
93 {
94  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
95 protected:
96  void update_(const ElementContext& elemCtx OPM_UNUSED,
97  unsigned dofIdx OPM_UNUSED,
98  unsigned timeIdx OPM_UNUSED)
99  { }
100 };
101 
118 template <class TypeTag>
119 class DarcyExtensiveQuantities
120 {
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;
128 
129  enum { dimWorld = GridView::dimensionworld };
130  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
131 
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;
137 
138 public:
143  const DimMatrix& intrinsicPermability() const
144  { return K_; }
145 
152  const EvalDimVector& potentialGrad(unsigned phaseIdx) const
153  { return potentialGrad_[phaseIdx]; }
154 
161  const EvalDimVector& filterVelocity(unsigned phaseIdx) const
162  { return filterVelocity_[phaseIdx]; }
163 
173  const Evaluation& volumeFlux(unsigned phaseIdx) const
174  { return volumeFlux_[phaseIdx]; }
175 
176 protected:
177  short upstreamIndex_(unsigned phaseIdx) const
178  { return upstreamDofIdx_[phaseIdx]; }
179 
180  short downstreamIndex_(unsigned phaseIdx) const
181  { return downstreamDofIdx_[phaseIdx]; }
182 
188  void calculateGradients_(const ElementContext& elemCtx,
189  unsigned faceIdx,
190  unsigned timeIdx)
191  {
192  const auto& gradCalc = elemCtx.gradientCalculator();
193  Ewoms::PressureCallback<TypeTag> pressureCallback(elemCtx);
194 
195  const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
196  const auto& faceNormal = scvf.normal();
197 
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();
203 
204  // calculate the "raw" pressure gradient
205  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
206  if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
207  Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
208  continue;
209  }
210 
211  pressureCallback.setPhaseIndex(phaseIdx);
212  gradCalc.calculateGradient(potentialGrad_[phaseIdx],
213  elemCtx,
214  faceIdx,
215  pressureCallback);
216  Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
217  }
218 
219  // correct the pressure gradients by the gravitational acceleration
220  if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) {
221  // estimate the gravitational acceleration at a given SCV face
222  // using the arithmetic mean
223  const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
224  const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
225 
226  const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
227  const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
228 
229  const auto& posIn = elemCtx.pos(i, timeIdx);
230  const auto& posEx = elemCtx.pos(j, timeIdx);
231  const auto& posFace = scvf.integrationPos();
232 
233  // the distance between the centers of the control volumes
234  DimVector distVecIn(posIn);
235  DimVector distVecEx(posEx);
236  DimVector distVecTotal(posEx);
237 
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))
244  continue;
245 
246  // calculate the hydrostatic pressure at the integration point of the face
247  Evaluation pStatIn;
248 
249  if (std::is_same<Scalar, Evaluation>::value ||
250  interiorDofIdx_ == static_cast<int>(focusDofIdx))
251  {
252  const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
253  pStatIn = - rhoIn*(gIn*distVecIn);
254  }
255  else {
256  Scalar rhoIn = Toolbox::value(intQuantsIn.fluidState().density(phaseIdx));
257  pStatIn = - rhoIn*(gIn*distVecIn);
258  }
259 
260  // the quantities on the exterior side of the face do not influence the
261  // result for the TPFA scheme, so they can be treated as scalar values.
262  Evaluation pStatEx;
263 
264  if (std::is_same<Scalar, Evaluation>::value ||
265  exteriorDofIdx_ == static_cast<int>(focusDofIdx))
266  {
267  const Evaluation& rhoEx = intQuantsEx.fluidState().density(phaseIdx);
268  pStatEx = - rhoEx*(gEx*distVecEx);
269  }
270  else {
271  Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
272  pStatEx = - rhoEx*(gEx*distVecEx);
273  }
274 
275  // compute the hydrostatic gradient between the two control volumes (this
276  // gradient exhibitis the same direction as the vector between the two
277  // control volume centers and the length (pStaticExterior -
278  // pStaticInterior)/distanceInteriorToExterior
279  Dune::FieldVector<Evaluation, dimWorld> f(distVecTotal);
280  f *= (pStatEx - pStatIn)/absDistTotalSquared;
281 
282  // calculate the final potential gradient
283  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
284  potentialGrad_[phaseIdx][dimIdx] += f[dimIdx];
285 
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) << "'");
291  }
292  }
293  }
294  }
295 
296  Opm::Valgrind::SetUndefined(K_);
297  elemCtx.problem().intersectionIntrinsicPermeability(K_, elemCtx, faceIdx, timeIdx);
298  Opm::Valgrind::CheckDefined(K_);
299 
300  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
301  if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
302  Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
303  continue;
304  }
305 
306  // determine the upstream and downstream DOFs
307  Evaluation tmp = 0.0;
308  for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
309  tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx];
310 
311  if (tmp > 0) {
312  upstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
313  downstreamDofIdx_[phaseIdx] = interiorDofIdx_;
314  }
315  else {
316  upstreamDofIdx_[phaseIdx] = interiorDofIdx_;
317  downstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
318  }
319 
320  // we only carry the derivatives along if the upstream DOF is the one which
321  // we currently focus on
322  const auto& up = elemCtx.intensiveQuantities(upstreamDofIdx_[phaseIdx], timeIdx);
323  if (upstreamDofIdx_[phaseIdx] == static_cast<int>(focusDofIdx))
324  mobility_[phaseIdx] = up.mobility(phaseIdx);
325  else
326  mobility_[phaseIdx] = Toolbox::value(up.mobility(phaseIdx));
327  }
328  }
329 
336  template <class FluidState>
337  void calculateBoundaryGradients_(const ElementContext& elemCtx,
338  unsigned boundaryFaceIdx,
339  unsigned timeIdx,
340  const FluidState& fluidState,
341  const typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache)
342  {
343  const auto& gradCalc = elemCtx.gradientCalculator();
344  Ewoms::BoundaryPressureCallback<TypeTag, FluidState> pressureCallback(elemCtx, fluidState);
345 
346  // calculate the pressure gradient
347  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
348  if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
349  Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
350  continue;
351  }
352 
353  pressureCallback.setPhaseIndex(phaseIdx);
354  gradCalc.calculateBoundaryGradient(potentialGrad_[phaseIdx],
355  elemCtx,
356  boundaryFaceIdx,
357  pressureCallback);
358  Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
359  }
360 
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();
366 
367  // calculate the intrinsic permeability
368  const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
369  K_ = intQuantsIn.intrinsicPermeability();
370 
371  // correct the pressure gradients by the gravitational acceleration
372  if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) {
373  // estimate the gravitational acceleration at a given SCV face
374  // using the arithmetic mean
375  const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
376  const auto& posIn = elemCtx.pos(i, timeIdx);
377  const auto& posFace = scvf.integrationPos();
378 
379  // the distance between the face center and the center of the control volume
380  DimVector distVecIn(posIn);
381  distVecIn -= posFace;
382  Scalar absDist = distVecIn.two_norm();
383  Scalar gTimesDist = gIn*distVecIn;
384 
385  for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
386  if (!elemCtx.model().phaseIsConsidered(phaseIdx))
387  continue;
388 
389  // calculate the hydrostatic pressure at the integration point of the face
390  Evaluation rhoIn = intQuantsIn.fluidState().density(phaseIdx);
391  Evaluation pStatIn = - gTimesDist*rhoIn;
392 
393  Opm::Valgrind::CheckDefined(pStatIn);
394 
395  // compute the hydrostatic gradient between the two control volumes (this
396  // gradient exhibitis the same direction as the vector between the two
397  // control volume centers and the length (pStaticExterior -
398  // pStaticInterior)/distanceInteriorToExterior. Note that for the
399  // boundary, 'pStaticExterior' is zero as the boundary pressure is
400  // defined on boundary face's integration point...
401  EvalDimVector f(distVecIn);
402  f *= - pStatIn/absDist;
403 
404  // calculate the final potential gradient
405  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
406  potentialGrad_[phaseIdx][dimIdx] += f[dimIdx];
407 
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) << "'");
414  }
415  }
416  }
417  }
418 
419  // determine the upstream and downstream DOFs
420  const auto& faceNormal = scvf.normal();
421 
422  const auto& matParams = elemCtx.problem().materialLawParams(elemCtx, i, timeIdx);
423 
424  Scalar kr[numPhases];
425  MaterialLaw::relativePermeabilities(kr, matParams, fluidState);
426  Opm::Valgrind::CheckDefined(kr);
427 
428  for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
429  if (!elemCtx.model().phaseIsConsidered(phaseIdx))
430  continue;
431 
432  Evaluation tmp = 0.0;
433  for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
434  tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx];
435 
436  if (tmp > 0) {
437  upstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
438  downstreamDofIdx_[phaseIdx] = interiorDofIdx_;
439  }
440  else {
441  upstreamDofIdx_[phaseIdx] = interiorDofIdx_;
442  downstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
443  }
444 
445  // take the phase mobility from the DOF in upstream direction
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));
451  else
452  mobility_[phaseIdx] = intQuantsIn.mobility(phaseIdx);
453  Opm::Valgrind::CheckDefined(mobility_[phaseIdx]);
454  }
455  }
456 
463  void calculateFluxes_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
464  {
465  const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
466  const DimVector& normal = scvf.normal();
467  Opm::Valgrind::CheckDefined(normal);
468 
469  for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
470  filterVelocity_[phaseIdx] = 0.0;
471  volumeFlux_[phaseIdx] = 0.0;
472  if (!elemCtx.model().phaseIsConsidered(phaseIdx))
473  continue;
474 
475  asImp_().calculateFilterVelocity_(phaseIdx);
476  Opm::Valgrind::CheckDefined(filterVelocity_[phaseIdx]);
477 
478  volumeFlux_[phaseIdx] = 0.0;
479  for (unsigned i = 0; i < normal.size(); ++i)
480  volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i];
481  }
482  }
483 
490  void calculateBoundaryFluxes_(const ElementContext& elemCtx,
491  unsigned boundaryFaceIdx,
492  unsigned timeIdx)
493  {
494  const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
495  const DimVector& normal = scvf.normal();
496  Opm::Valgrind::CheckDefined(normal);
497 
498  for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
499  if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
500  filterVelocity_[phaseIdx] = 0.0;
501  volumeFlux_[phaseIdx] = 0.0;
502  continue;
503  }
504 
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];
510  }
511  }
512 
513  void calculateFilterVelocity_(unsigned phaseIdx)
514  {
515 #ifndef NDEBUG
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]));
520 #endif
521 
522  K_.mv(potentialGrad_[phaseIdx], filterVelocity_[phaseIdx]);
523  filterVelocity_[phaseIdx] *= - mobility_[phaseIdx];
524 
525 #ifndef NDEBUG
526  for (unsigned i = 0; i < filterVelocity_[phaseIdx].size(); ++ i)
527  assert(std::isfinite(Toolbox::value(filterVelocity_[phaseIdx][i])));
528 #endif
529  }
530 
531 private:
532  Implementation& asImp_()
533  { return *static_cast<Implementation*>(this); }
534 
535  const Implementation& asImp_() const
536  { return *static_cast<const Implementation*>(this); }
537 
538 protected:
539  // intrinsic permeability tensor and its square root
540  DimMatrix K_;
541 
542  // mobilities of all fluid phases [1 / (Pa s)]
543  Evaluation mobility_[numPhases];
544 
545  // filter velocities of all phases [m/s]
546  EvalDimVector filterVelocity_[numPhases];
547 
548  // the volumetric flux of all fluid phases over the control
549  // volume's face [m^3/s / m^2]
550  Evaluation volumeFlux_[numPhases];
551 
552  // pressure potential gradients of all phases [Pa / m]
553  EvalDimVector potentialGrad_[numPhases];
554 
555  // upstream, downstream, interior and exterior DOFs
556  short upstreamDofIdx_[numPhases];
557  short downstreamDofIdx_[numPhases];
558  short interiorDofIdx_;
559  short exteriorDofIdx_;
560 };
561 
562 } // namespace Ewoms
563 
564 #endif
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&#39;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&#39;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 > &paramCache)
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&#39;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