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 OPM_DIFFUSION_MODULE_HH
29#define OPM_DIFFUSION_MODULE_HH
30
31#include <dune/common/fvector.hh>
32
33#include <opm/material/common/Valgrind.hpp>
34
38
39#include <array>
40#include <stdexcept>
41
42namespace Opm {
43
50template <class TypeTag, bool enableDiffusion>
52
56template <class TypeTag>
57class DiffusionModule<TypeTag, /*enableDiffusion=*/false>
58{
62
63public:
67 static void registerParameters()
68 {}
69
74 template <class Context>
75 static void addDiffusiveFlux(RateVector&,
76 const Context&,
77 unsigned,
78 unsigned)
79 {}
80};
81
85template <class TypeTag>
86class DiffusionModule<TypeTag, /*enableDiffusion=*/true>
87{
93
94 enum { numPhases = FluidSystem::numPhases };
95 enum { numComponents = FluidSystem::numComponents };
96 enum { conti0EqIdx = Indices::conti0EqIdx };
97
98 using Toolbox = Opm::MathToolbox<Evaluation>;
99
100public:
104 static void registerParameters()
105 {}
106
111 template <class Context>
112 static void addDiffusiveFlux(RateVector& flux, const Context& context,
113 unsigned spaceIdx, unsigned timeIdx)
114 {
115 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
116
117 const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
118 const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
119
120 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
121 // arithmetic mean of the phase's molar density
122 Evaluation rhoMolar = fluidStateI.molarDensity(phaseIdx);
123 rhoMolar += Toolbox::value(fluidStateJ.molarDensity(phaseIdx));
124 rhoMolar /= 2;
125
126 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
127 // mass flux due to molecular diffusion
128 flux[conti0EqIdx + compIdx] +=
129 -rhoMolar
130 * extQuants.moleFractionGradientNormal(phaseIdx, compIdx)
131 * extQuants.effectiveDiffusionCoefficient(phaseIdx, compIdx);
132 }
133 }
134 }
135};
136
144template <class TypeTag, bool enableDiffusion>
146
150template <class TypeTag>
151class DiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/false>
152{
156
157public:
162 Scalar tortuosity(unsigned) const
163 {
164 throw std::logic_error("Method tortuosity() does not make sense "
165 "if diffusion is disabled");
166 }
167
172 Scalar diffusionCoefficient(unsigned, unsigned) const
173 {
174 throw std::logic_error("Method diffusionCoefficient() does not "
175 "make sense if diffusion is disabled");
176 }
177
182 Scalar effectiveDiffusionCoefficient(unsigned, unsigned) const
183 {
184 throw std::logic_error("Method effectiveDiffusionCoefficient() "
185 "does not make sense if diffusion is disabled");
186 }
187
188protected:
193 template <class FluidState>
194 void update_(FluidState&,
195 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>&,
196 const ElementContext&,
197 unsigned,
198 unsigned)
199 { }
200};
201
205template <class TypeTag>
206class DiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
207{
212
213 enum { numPhases = FluidSystem::numPhases };
214 enum { numComponents = FluidSystem::numComponents };
215
216public:
221 Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
222 { return diffusionCoefficient_[phaseIdx][compIdx]; }
223
228 Evaluation tortuosity(unsigned phaseIdx) const
229 { return tortuosity_[phaseIdx]; }
230
235 Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
236 { return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx]; }
237
238protected:
243 template <class FluidState>
244 void update_(FluidState& fluidState,
245 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
246 const ElementContext& elemCtx,
247 unsigned dofIdx,
248 unsigned timeIdx)
249 {
250 using Toolbox = Opm::MathToolbox<Evaluation>;
251
252 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
253 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
254 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
255 continue;
256 }
257
258 // TODO: let the problem do this (this is a constitutive
259 // relation of which the model should be free of from the
260 // abstraction POV!)
261 // Based on Millington, R. J., & Quirk, J. P. (1961).
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, 10.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
280private:
281 std::array<Evaluation, numPhases> tortuosity_;
282 std::array<std::array<Evaluation, numComponents>, numPhases> diffusionCoefficient_;
283};
284
291template <class TypeTag, bool enableDiffusion>
293
297template <class TypeTag>
298class DiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/false>
299{
303
304protected:
309 void update_(const ElementContext&,
310 unsigned,
311 unsigned)
312 {}
313
314 template <class Context, class FluidState>
315 void updateBoundary_(const Context&,
316 unsigned,
317 unsigned,
318 const FluidState&)
319 {}
320
321public:
328 const Evaluation& moleFractionGradientNormal(unsigned,
329 unsigned) const
330 {
331 throw std::logic_error("The method moleFractionGradient() does not "
332 "make sense if diffusion is disabled.");
333 }
334
342 const Evaluation& effectiveDiffusionCoefficient(unsigned,
343 unsigned) const
344 {
345 throw std::logic_error("The method effectiveDiffusionCoefficient() "
346 "does not make sense if diffusion is disabled.");
347 }
348};
349
353template <class TypeTag>
354class DiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
355{
360
361 enum { dimWorld = GridView::dimensionworld };
362 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
363 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
364
365 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
366 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
367
368protected:
373 void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx)
374 {
375 const auto& gradCalc = elemCtx.gradientCalculator();
376 Opm::MoleFractionCallback<TypeTag> moleFractionCallback(elemCtx);
377
378 const auto& face = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
379 const auto& normal = face.normal();
380 const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
381
382 const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
383 const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
384
385 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
386 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
387 continue;
388 }
389
390 moleFractionCallback.setPhaseIndex(phaseIdx);
391 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
392 moleFractionCallback.setComponentIndex(compIdx);
393
394 DimEvalVector moleFractionGradient(0.0);
395 gradCalc.calculateGradient(moleFractionGradient,
396 elemCtx,
397 faceIdx,
398 moleFractionCallback);
399
400 moleFractionGradientNormal_[phaseIdx][compIdx] = 0.0;
401 for (unsigned i = 0; i < normal.size(); ++i) {
402 moleFractionGradientNormal_[phaseIdx][compIdx] +=
403 normal[i]*moleFractionGradient[i];
404 }
405 Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
406
407 // use the arithmetic average for the effective
408 // diffusion coefficients.
409 effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
410 (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx)
411 +
412 intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx))
413 / 2;
414
415 Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
416 }
417 }
418 }
419
420 template <class Context, class FluidState>
421 void updateBoundary_(const Context& context,
422 unsigned bfIdx,
423 unsigned timeIdx,
424 const FluidState& fluidState)
425 {
426 const auto& stencil = context.stencil(timeIdx);
427 const auto& face = stencil.boundaryFace(bfIdx);
428
429 const auto& elemCtx = context.elementContext();
430 const unsigned insideScvIdx = face.interiorIndex();
431 const auto& insideScv = stencil.subControlVolume(insideScvIdx);
432
433 const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
434 const auto& fluidStateInside = intQuantsInside.fluidState();
435
436 // distance between the center of the SCV and center of the boundary face
437 DimVector distVec = face.integrationPos();
438 distVec -= context.element().geometry().global(insideScv.localGeometry().center());
439
440 const Scalar dist = distVec * face.normal();
441
442 // if the following assertation triggers, the center of the
443 // center of the interior SCV was not inside the element!
444 assert(dist > 0);
445
446 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
447 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
448 continue;
449 }
450
451 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
452 // calculate mole fraction gradient using two-point
453 // gradients
454 moleFractionGradientNormal_[phaseIdx][compIdx] =
455 (fluidState.moleFraction(phaseIdx, compIdx)
456 -
457 fluidStateInside.moleFraction(phaseIdx, compIdx))
458 / dist;
459 Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
460
461 // use effective diffusion coefficients of the interior finite
462 // volume.
463 effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
464 intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx);
465
466 Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
467 }
468 }
469 }
470
471public:
478 const Evaluation& moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
479 { return moleFractionGradientNormal_[phaseIdx][compIdx]; }
480
488 const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
489 { return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
490
491private:
492 std::array<std::array<Evaluation, numComponents>, numPhases> moleFractionGradientNormal_;
493 std::array<std::array<Evaluation, numComponents>, numPhases> effectiveDiffusionCoefficient_;
494};
495
496} // namespace Opm
497
498#endif
void update_(const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:309
const Evaluation & effectiveDiffusionCoefficient(unsigned, unsigned) const
The effective diffusion coeffcient of a component in a fluid phase at the face's integration point.
Definition: diffusionmodule.hh:342
const Evaluation & moleFractionGradientNormal(unsigned, unsigned) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:328
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: diffusionmodule.hh:315
void updateBoundary_(const Context &context, unsigned bfIdx, unsigned timeIdx, const FluidState &fluidState)
Definition: diffusionmodule.hh:421
const Evaluation & moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:478
const Evaluation & effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
The effective diffusion coeffcient of a component in a fluid phase at the face's integration point.
Definition: diffusionmodule.hh:488
void update_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:373
Provides the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:292
Scalar diffusionCoefficient(unsigned, unsigned) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:172
Scalar effectiveDiffusionCoefficient(unsigned, unsigned) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: diffusionmodule.hh:182
void update_(FluidState &, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &, const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:194
Scalar tortuosity(unsigned) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:162
Evaluation tortuosity(unsigned phaseIdx) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:228
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:244
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:235
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:221
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes.
Definition: diffusionmodule.hh:145
static void addDiffusiveFlux(RateVector &, const Context &, unsigned, unsigned)
Adds the diffusive mass flux flux to the flux vector over a flux integration point.
Definition: diffusionmodule.hh:75
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: diffusionmodule.hh:67
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:112
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: diffusionmodule.hh:104
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:51
Callback class for a mole fraction of a component in a phase.
Definition: quantitycallbacks.hh:426
void setComponentIndex(unsigned compIdx)
Set the index of the component for which the mole fraction should be returned.
Definition: quantitycallbacks.hh:461
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which a mole fraction should be returned.
Definition: quantitycallbacks.hh:454
Declare the properties used by the infrastructure code of the finite volume discretizations.
Defines the common properties required by the porous medium multi-phase models.
Definition: blackoilboundaryratevector.hh:39
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:233
This method contains all callback classes for quantities that are required by some extensive quantiti...