All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
diffusionmodule.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 */
28 #ifndef EWOMS_DIFFUSION_MODULE_HH
29 #define EWOMS_DIFFUSION_MODULE_HH
30 
33 
34 #include <opm/common/Valgrind.hpp>
35 #include <opm/common/Unused.hpp>
36 #include <opm/common/ErrorMacros.hpp>
37 #include <opm/common/Exceptions.hpp>
38 
39 #include <dune/common/fvector.hh>
40 
41 namespace Ewoms {
42 namespace Properties {
43 NEW_PROP_TAG(Indices);
44 }
45 
52 template <class TypeTag, bool enableDiffusion>
54 
58 template <class TypeTag>
59 class DiffusionModule<TypeTag, /*enableDiffusion=*/false>
60 {
61  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
62  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
63  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
64 
65 public:
69  static void registerParameters()
70  {}
71 
76  template <class Context>
77  static void addDiffusiveFlux(RateVector& flux OPM_UNUSED,
78  const Context& context OPM_UNUSED,
79  unsigned spaceIdx OPM_UNUSED,
80  unsigned timeIdx OPM_UNUSED)
81  {}
82 };
83 
87 template <class TypeTag>
88 class DiffusionModule<TypeTag, /*enableDiffusion=*/true>
89 {
90  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
91  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
92  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
93  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
94  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
95 
96  enum { numPhases = FluidSystem::numPhases };
97  enum { numComponents = FluidSystem::numComponents };
98  enum { conti0EqIdx = Indices::conti0EqIdx };
99 
100  typedef Opm::MathToolbox<Evaluation> Toolbox;
101 
102 public:
106  static void registerParameters()
107  {}
108 
113  template <class Context>
114  static void addDiffusiveFlux(RateVector& flux, const Context& context,
115  unsigned spaceIdx, unsigned timeIdx)
116  {
117  const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
118 
119  const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
120  const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
121 
122  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
123  // arithmetic mean of the phase's molar density
124  Evaluation rhoMolar = fluidStateI.molarDensity(phaseIdx);
125  rhoMolar += Toolbox::value(fluidStateJ.molarDensity(phaseIdx));
126  rhoMolar /= 2;
127 
128  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
129  // mass flux due to molecular diffusion
130  flux[conti0EqIdx + compIdx] +=
131  -rhoMolar
132  * extQuants.moleFractionGradientNormal(phaseIdx, compIdx)
133  * extQuants.effectiveDiffusionCoefficient(phaseIdx, compIdx);
134  }
135  }
136 };
137 
145 template <class TypeTag, bool enableDiffusion>
147 
151 template <class TypeTag>
152 class DiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/false>
153 {
154  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
155  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
156  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
157 
158 public:
163  Scalar tortuosity(unsigned phaseIdx OPM_UNUSED) const
164  {
165  OPM_THROW(std::logic_error, "Method tortuosity() does not make sense "
166  "if diffusion is disabled");
167  }
168 
173  Scalar diffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
174  {
175  OPM_THROW(std::logic_error, "Method diffusionCoefficient() does not "
176  "make sense if diffusion is disabled");
177  }
178 
183  Scalar effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
184  {
185  OPM_THROW(std::logic_error, "Method effectiveDiffusionCoefficient() "
186  "does not make sense if diffusion is "
187  "disabled");
188  }
189 
190 protected:
195  template <class FluidState>
196  void update_(FluidState& fs OPM_UNUSED,
197  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache OPM_UNUSED,
198  const ElementContext& elemCtx OPM_UNUSED,
199  unsigned dofIdx OPM_UNUSED,
200  unsigned timeIdx OPM_UNUSED)
201  { }
202 };
203 
207 template <class TypeTag>
208 class DiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
209 {
210  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
211  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
212  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
213  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
214 
215  enum { numPhases = FluidSystem::numPhases };
216  enum { numComponents = FluidSystem::numComponents };
217 
218 public:
223  Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
224  { return diffusionCoefficient_[phaseIdx][compIdx]; }
225 
230  Evaluation tortuosity(unsigned phaseIdx) const
231  { return tortuosity_[phaseIdx]; }
232 
237  Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
238  { return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx]; }
239 
240 protected:
245  template <class FluidState>
246  void update_(FluidState& fluidState,
247  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
248  const ElementContext& elemCtx,
249  unsigned dofIdx,
250  unsigned timeIdx)
251  {
252  typedef Opm::MathToolbox<Evaluation> Toolbox;
253 
254  const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
255  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
256  if (!elemCtx.model().phaseIsConsidered(phaseIdx))
257  continue;
258 
259  // TODO: let the problem do this (this is a constitutive
260  // relation of which the model should be free of from the
261  // abstraction POV!)
262  const Evaluation& base =
263  Toolbox::max(0.0001,
264  intQuants.porosity()
265  * intQuants.fluidState().saturation(phaseIdx));
266  tortuosity_[phaseIdx] =
267  1.0 / (intQuants.porosity() * intQuants.porosity())
268  * Toolbox::pow(base, 7.0/3.0);
269 
270  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
271  diffusionCoefficient_[phaseIdx][compIdx] =
272  FluidSystem::diffusionCoefficient(fluidState,
273  paramCache,
274  phaseIdx,
275  compIdx);
276  }
277  }
278  }
279 
280 private:
281  Evaluation tortuosity_[numPhases];
282  Evaluation diffusionCoefficient_[numPhases][numComponents];
283 };
284 
291 template <class TypeTag, bool enableDiffusion>
293 
297 template <class TypeTag>
298 class DiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/false>
299 {
300  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
301  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
302  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
303 
304 protected:
309  void update_(const ElementContext& elemCtx OPM_UNUSED,
310  unsigned faceIdx OPM_UNUSED,
311  unsigned timeIdx OPM_UNUSED)
312  {}
313 
314  template <class Context, class FluidState>
315  void updateBoundary_(const Context& context OPM_UNUSED,
316  unsigned bfIdx OPM_UNUSED,
317  unsigned timeIdx OPM_UNUSED,
318  const FluidState& fluidState OPM_UNUSED)
319  {}
320 
321 public:
328  const Evaluation& moleFractionGradientNormal(unsigned phaseIdx OPM_UNUSED,
329  unsigned compIdx OPM_UNUSED) const
330  {
331  OPM_THROW(std::logic_error,
332  "The method moleFractionGradient() does not "
333  "make sense if diffusion is disabled.");
334  }
335 
343  const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED,
344  unsigned compIdx OPM_UNUSED) const
345  {
346  OPM_THROW(std::logic_error,
347  "The method effectiveDiffusionCoefficient() "
348  "does not make sense if diffusion is disabled.");
349  }
350 };
351 
355 template <class TypeTag>
356 class DiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
357 {
358  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
359  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
360  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
361  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
362 
363  enum { dimWorld = GridView::dimensionworld };
364  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
365  enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
366 
367  typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
368  typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
369 
370 protected:
375  void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx)
376  {
377  const auto& gradCalc = elemCtx.gradientCalculator();
378  Ewoms::MoleFractionCallback<TypeTag> moleFractionCallback(elemCtx);
379 
380  DimEvalVector temperatureGrad;
381 
382  const auto& face = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
383  const auto& normal = face.normal();
384  const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
385 
386  const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
387  const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
388 
389  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
390  if (!elemCtx.model().phaseIsConsidered(phaseIdx))
391  continue;
392 
393  moleFractionCallback.setPhaseIndex(phaseIdx);
394  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
395  moleFractionCallback.setComponentIndex(compIdx);
396 
397  DimEvalVector moleFractionGradient(0.0);
398  gradCalc.calculateGradient(moleFractionGradient,
399  elemCtx,
400  faceIdx,
401  moleFractionCallback);
402 
403  moleFractionGradientNormal_[phaseIdx][compIdx] = 0.0;
404  for (unsigned i = 0; i < normal.size(); ++i)
405  moleFractionGradientNormal_[phaseIdx][compIdx] +=
406  normal[i]*moleFractionGradient[i];
407  Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
408 
409  // use the arithmetic average for the effective
410  // diffusion coefficients.
411  effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
412  (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx)
413  +
414  intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx))
415  / 2;
416 
417  Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
418  }
419  }
420  }
421 
422  template <class Context, class FluidState>
423  void updateBoundary_(const Context& context,
424  unsigned bfIdx,
425  unsigned timeIdx,
426  const FluidState& fluidState)
427  {
428  const auto& stencil = context.stencil(timeIdx);
429  const auto& face = stencil.boundaryFace(bfIdx);
430 
431  const auto& elemCtx = context.elementContext();
432  unsigned insideScvIdx = face.interiorIndex();
433  const auto& insideScv = stencil.subControlVolume(insideScvIdx);
434 
435  const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
436  const auto& fluidStateInside = intQuantsInside.fluidState();
437 
438  // distance between the center of the SCV and center of the boundary face
439  DimVector distVec = face.integrationPos();
440  distVec -= context.element().geometry().global(insideScv.localGeometry().center());
441 
442  Scalar dist = distVec * face.normal();
443 
444  // if the following assertation triggers, the center of the
445  // center of the interior SCV was not inside the element!
446  assert(dist > 0);
447 
448  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
449  if (!elemCtx.model().phaseIsConsidered(phaseIdx))
450  continue;
451 
452  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
453  // calculate mole fraction gradient using two-point
454  // gradients
455  moleFractionGradientNormal_[phaseIdx][compIdx] =
456  (fluidState.moleFraction(phaseIdx, compIdx)
457  -
458  fluidStateInside.moleFraction(phaseIdx, compIdx))
459  / dist;
460  Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
461 
462  // use effective diffusion coefficients of the interior finite
463  // volume.
464  effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
465  intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx);
466 
467  Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
468  }
469  }
470  }
471 
472 public:
479  const Evaluation& moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
480  { return moleFractionGradientNormal_[phaseIdx][compIdx]; }
481 
489  const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
490  { return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
491 
492 private:
493  Evaluation moleFractionGradientNormal_[numPhases][numComponents];
494  Evaluation effectiveDiffusionCoefficient_[numPhases][numComponents];
495 };
496 
497 } // namespace Ewoms
498 
499 #endif
Scalar tortuosity(unsigned phaseIdx OPM_UNUSED) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:163
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes...
Definition: diffusionmodule.hh:146
Scalar diffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:173
static void addDiffusiveFlux(RateVector &flux, const Context &context, unsigned spaceIdx, unsigned timeIdx)
Adds the mass flux due to molecular diffusion to the flux vector over the flux integration point...
Definition: diffusionmodule.hh:114
Scalar effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: diffusionmodule.hh:183
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const Evaluation & effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
The effective diffusion coeffcient of a component in a fluid phase at the face&#39;s integration point...
Definition: diffusionmodule.hh:489
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
const Evaluation & moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:479
void update_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:375
Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: diffusionmodule.hh:237
Declare the properties used by the infrastructure code of the finite volume discretizations.
Evaluation tortuosity(unsigned phaseIdx) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:230
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:53
const Evaluation & moleFractionGradientNormal(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:328
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: diffusionmodule.hh:106
This method contains all callback classes for quantities that are required by some extensive quantiti...
const Evaluation & effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
The effective diffusion coeffcient of a component in a fluid phase at the face&#39;s integration point...
Definition: diffusionmodule.hh:343
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:223
void setComponentIndex(unsigned compIdx)
Set the index of the component for which the mole fraction should be returned.
Definition: quantitycallbacks.hh:460
void update_(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:246
static void addDiffusiveFlux(RateVector &flux OPM_UNUSED, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Adds the diffusive mass flux flux to the flux vector over a flux integration point.
Definition: diffusionmodule.hh:77
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which a mole fraction should be returned.
Definition: quantitycallbacks.hh:453
Callback class for a mole fraction of a component in a phase.
Definition: quantitycallbacks.hh:424
void update_(const ElementContext &elemCtx OPM_UNUSED, unsigned faceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:309
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: diffusionmodule.hh:69
Provides the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:292
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
void update_(FluidState &fs OPM_UNUSED, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache OPM_UNUSED, const ElementContext &elemCtx OPM_UNUSED, unsigned dofIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Update the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:196