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
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes...
Definition: diffusionmodule.hh:146
Definition: baseauxiliarymodule.hh:37
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
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
Evaluation tortuosity(unsigned phaseIdx) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:230
const Evaluation & moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:479
#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 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
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
void update_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:375
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
Declare the properties used by the infrastructure code of the finite volume discretizations.
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:53
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 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...
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
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
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
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
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
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
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:223