3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
flux/cctpfa/darcyslaw.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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program 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 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH
25#define DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH
26
27#include <dumux/common/math.hh>
30
33
34namespace Dumux {
35
36// forward declarations
37template<class TypeTag, DiscretizationMethod discMethod>
38class DarcysLawImplementation;
39
48template<class Scalar, class GridGeometry, bool isNetwork>
50
56template <class TypeTag>
58: public CCTpfaDarcysLaw<GetPropType<TypeTag, Properties::Scalar>,
59 GetPropType<TypeTag, Properties::GridGeometry>,
60 (GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension < GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld)>
61{};
62
67template<class GridGeometry>
68class TpfaDarcysLawCacheFiller
69{
70 using FVElementGeometry = typename GridGeometry::LocalView;
71 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
72 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
73
74public:
78 template<class FluxVariablesCache, class Problem, class ElementVolumeVariables, class FluxVariablesCacheFiller>
79 static void fill(FluxVariablesCache& scvfFluxVarsCache,
80 const Problem& problem,
81 const Element& element,
82 const FVElementGeometry& fvGeometry,
83 const ElementVolumeVariables& elemVolVars,
84 const SubControlVolumeFace& scvf,
85 const FluxVariablesCacheFiller& fluxVarsCacheFiller)
86 {
87 scvfFluxVarsCache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf);
88 }
89};
90
95template<class AdvectionType, class GridGeometry>
96class TpfaDarcysLawCache
97{
98 using Scalar = typename AdvectionType::Scalar;
99 using FVElementGeometry = typename GridGeometry::LocalView;
100 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
101 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
102
103public:
104 using Filler = TpfaDarcysLawCacheFiller<GridGeometry>;
105
106 template<class Problem, class ElementVolumeVariables>
107 void updateAdvection(const Problem& problem,
108 const Element& element,
109 const FVElementGeometry& fvGeometry,
110 const ElementVolumeVariables& elemVolVars,
111 const SubControlVolumeFace &scvf)
112 {
113 tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf);
114 }
115
116 const Scalar& advectionTij() const
117 { return tij_; }
118
119private:
120 Scalar tij_;
121};
122
127template<class ScalarType, class GridGeometry>
128class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
129{
130 using ThisType = CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>;
131 using FVElementGeometry = typename GridGeometry::LocalView;
132 using SubControlVolume = typename GridGeometry::SubControlVolume;
133 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
134 using GridView = typename GridGeometry::GridView;
135 using Element = typename GridView::template Codim<0>::Entity;
136
137 static constexpr int dim = GridView::dimension;
138 static constexpr int dimWorld = GridView::dimensionworld;
139
140 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
141
142 public:
144 using Scalar = ScalarType;
145
147 static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
148
150 using Cache = TpfaDarcysLawCache<ThisType, GridGeometry>;
151
153 template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
154 static Scalar flux(const Problem& problem,
155 const Element& element,
156 const FVElementGeometry& fvGeometry,
157 const ElementVolumeVariables& elemVolVars,
158 const SubControlVolumeFace& scvf,
159 int phaseIdx,
160 const ElementFluxVarsCache& elemFluxVarsCache)
161 {
162 static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
163
164 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
165
166 // Get the inside and outside volume variables
167 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
168 const auto& insideVolVars = elemVolVars[insideScv];
169 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
170
171 if (enableGravity)
172 {
173 // do averaging for the density over all neighboring elements
174 const auto rho = scvf.boundary() ? outsideVolVars.density(phaseIdx)
175 : (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5;
176
177 // Obtain inside and outside pressures
178 const auto pInside = insideVolVars.pressure(phaseIdx);
179 const auto pOutside = outsideVolVars.pressure(phaseIdx);
180
181 const auto& tij = fluxVarsCache.advectionTij();
182 const auto& g = problem.spatialParams().gravity(scvf.ipGlobal());
183
185 const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
186
187 Scalar flux = tij*(pInside - pOutside) + rho*scvf.area()*alpha_inside;
188
190 if (!scvf.boundary())
191 {
192 const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
193 const auto outsideK = outsideVolVars.permeability();
194 const auto outsideTi = fvGeometry.gridGeometry().isPeriodic()
195 ? computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, outsideK, outsideVolVars.extrusionFactor())
196 : -1.0*computeTpfaTransmissibility(scvf, outsideScv, outsideK, outsideVolVars.extrusionFactor());
197 const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideK, g)*outsideVolVars.extrusionFactor();
198
199 flux -= rho*tij/outsideTi*(alpha_inside - alpha_outside);
200 }
201
202 return flux;
203 }
204 else
205 {
206 // Obtain inside and outside pressures
207 const auto pInside = insideVolVars.pressure(phaseIdx);
208 const auto pOutside = outsideVolVars.pressure(phaseIdx);
209
210 // return flux
211 return fluxVarsCache.advectionTij()*(pInside - pOutside);
212 }
213 }
214
215 // The flux variables cache has to be bound to an element prior to flux calculations
216 // During the binding, the transmissibility will be computed and stored using the method below.
217 template<class Problem, class ElementVolumeVariables>
218 static Scalar calculateTransmissibility(const Problem& problem,
219 const Element& element,
220 const FVElementGeometry& fvGeometry,
221 const ElementVolumeVariables& elemVolVars,
222 const SubControlVolumeFace& scvf)
223 {
224 Scalar tij;
225
226 const auto insideScvIdx = scvf.insideScvIdx();
227 const auto& insideScv = fvGeometry.scv(insideScvIdx);
228 const auto& insideVolVars = elemVolVars[insideScvIdx];
229
230 const Scalar ti = computeTpfaTransmissibility(scvf, insideScv,
231 getPermeability_(problem, insideVolVars, scvf.ipGlobal()),
232 insideVolVars.extrusionFactor());
233
234 // on the boundary (dirichlet) we only need ti
235 if (scvf.boundary())
236 tij = scvf.area()*ti;
237
238 // otherwise we compute a tpfa harmonic mean
239 else
240 {
241 const auto outsideScvIdx = scvf.outsideScvIdx();
242 // as we assemble fluxes from the neighbor to our element the outside index
243 // refers to the scv of our element, so we use the scv method
244 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
245 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
246 const Scalar tj = fvGeometry.gridGeometry().isPeriodic()
247 ? computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor())
248 : -1.0*computeTpfaTransmissibility(scvf, outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor());
249
250 // harmonic mean (check for division by zero!)
251 // TODO: This could lead to problems!? Is there a better way to do this?
252 if (ti*tj <= 0.0)
253 tij = 0;
254 else
255 tij = scvf.area()*(ti * tj)/(ti + tj);
256 }
257
258 return tij;
259 }
260
261private:
262 template<class Problem, class VolumeVariables,
263 std::enable_if_t<!Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
264 static decltype(auto) getPermeability_(const Problem& problem,
265 const VolumeVariables& volVars,
266 const GlobalPosition& scvfIpGlobal)
267 { return volVars.permeability(); }
268
269 template<class Problem, class VolumeVariables,
270 std::enable_if_t<Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
271 static decltype(auto) getPermeability_(const Problem& problem,
272 const VolumeVariables& volVars,
273 const GlobalPosition& scvfIpGlobal)
274 { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); }
275};
276
281template<class ScalarType, class GridGeometry>
282class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ true>
283{
284 using ThisType = CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ true>;
285 using FVElementGeometry = typename GridGeometry::LocalView;
286 using SubControlVolume = typename GridGeometry::SubControlVolume;
287 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
288 using GridView = typename GridGeometry::GridView;
289 using Element = typename GridView::template Codim<0>::Entity;
290
291 static constexpr int dim = GridView::dimension;
292 static constexpr int dimWorld = GridView::dimensionworld;
293
294 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
295
296public:
298 using Scalar = ScalarType;
299
301 static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
302
304 using Cache = TpfaDarcysLawCache<ThisType, GridGeometry>;
305
307 template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
308 static Scalar flux(const Problem& problem,
309 const Element& element,
310 const FVElementGeometry& fvGeometry,
311 const ElementVolumeVariables& elemVolVars,
312 const SubControlVolumeFace& scvf,
313 int phaseIdx,
314 const ElementFluxVarsCache& elemFluxVarsCache)
315 {
316 static const bool gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
317
318 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
319
320 // Get the inside and outside volume variables
321 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
322 const auto& insideVolVars = elemVolVars[insideScv];
323 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
324
325 if (gravity)
326 {
327 // do averaging for the density over all neighboring elements
328 const auto rho = [&]()
329 {
330 // boundaries
331 if (scvf.boundary())
332 return outsideVolVars.density(phaseIdx);
333
334 // inner faces with two neighboring elements
335 else if (scvf.numOutsideScvs() == 1)
336 return (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5;
337
338 // inner faces in networks (general case)
339 else
340 {
341 Scalar rho(insideVolVars.density(phaseIdx));
342 for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
343 {
344 const auto outsideScvIdx = scvf.outsideScvIdx(i);
345 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
346 rho += outsideVolVars.density(phaseIdx);
347 }
348 return rho/(scvf.numOutsideScvs()+1);
349 }
350 }();
351
352 const auto& tij = fluxVarsCache.advectionTij();
353 const auto& g = problem.spatialParams().gravity(scvf.ipGlobal());
354
355 // Obtain inside and outside pressures
356 const auto pInside = insideVolVars.pressure(phaseIdx);
357 const auto pOutside = [&]()
358 {
359 // Dirichlet boundaries and inner faces with one neighbor
360 if (scvf.numOutsideScvs() == 1)
361 return outsideVolVars.pressure(phaseIdx);
362
363 // inner faces in networks (general case)
364 else
365 {
366 Scalar sumTi(tij);
367 Scalar sumPTi(tij*pInside);
368
369 // add inside gravitational contribution
370 sumPTi += rho*scvf.area()
371 *insideVolVars.extrusionFactor()
372 *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g);
373
374 for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
375 {
376 const auto outsideScvIdx = scvf.outsideScvIdx(i);
377 const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
378 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
379 const auto& outsideFluxVarsCache = elemFluxVarsCache[flippedScvf];
380 sumTi += outsideFluxVarsCache.advectionTij();
381 sumPTi += outsideFluxVarsCache.advectionTij()*outsideVolVars.pressure(phaseIdx);
382
383 // add outside gravitational contribution
384 sumPTi += rho*scvf.area()
385 *outsideVolVars.extrusionFactor()
386 *vtmv(flippedScvf.unitOuterNormal(), outsideVolVars.permeability(), g);
387 }
388 return sumPTi/sumTi;
389 }
390 }();
391
393 const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
394
395 Scalar flux = tij*(pInside - pOutside) + scvf.area()*rho*alpha_inside;
396
398 if (!scvf.boundary() && scvf.numOutsideScvs() == 1)
399 {
400 const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
401 const auto& outsideScvf = fvGeometry.flipScvf(scvf.index());
402 const auto outsideK = outsideVolVars.permeability();
403 const auto outsideTi = computeTpfaTransmissibility(outsideScvf, outsideScv, outsideK, outsideVolVars.extrusionFactor());
404 const auto alpha_outside = vtmv(outsideScvf.unitOuterNormal(), outsideK, g)*outsideVolVars.extrusionFactor();
405
406 flux -= rho*tij/outsideTi*(alpha_inside + alpha_outside);
407 }
408
409 return flux;
410 }
411 else
412 {
413 // Obtain inside and outside pressures
414 const auto pInside = insideVolVars.pressure(phaseIdx);
415 const auto pOutside = [&]()
416 {
417 // Dirichlet boundaries and inner faces with two neighboring elements
418 if (scvf.numOutsideScvs() <= 1)
419 return outsideVolVars.pressure(phaseIdx);
420
421 // inner faces in networks (general case)
422 else
423 {
424 const auto& insideFluxVarsCache = elemFluxVarsCache[scvf];
425 Scalar sumTi(insideFluxVarsCache.advectionTij());
426 Scalar sumPTi(insideFluxVarsCache.advectionTij()*pInside);
427
428 for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
429 {
430 const auto outsideScvIdx = scvf.outsideScvIdx(i);
431 const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
432 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
433 const auto& outsideFluxVarsCache = elemFluxVarsCache[flippedScvf];
434 sumTi += outsideFluxVarsCache.advectionTij();
435 sumPTi += outsideFluxVarsCache.advectionTij()*outsideVolVars.pressure(phaseIdx);
436 }
437 return sumPTi/sumTi;
438 }
439 }();
440
441 // return flux
442 return fluxVarsCache.advectionTij()*(pInside - pOutside);
443 }
444 }
445
446 // The flux variables cache has to be bound to an element prior to flux calculations
447 // During the binding, the transmissibility will be computed and stored using the method below.
448 template<class Problem, class ElementVolumeVariables>
449 static Scalar calculateTransmissibility(const Problem& problem,
450 const Element& element,
451 const FVElementGeometry& fvGeometry,
452 const ElementVolumeVariables& elemVolVars,
453 const SubControlVolumeFace& scvf)
454 {
455 Scalar tij;
456
457 const auto insideScvIdx = scvf.insideScvIdx();
458 const auto& insideScv = fvGeometry.scv(insideScvIdx);
459 const auto& insideVolVars = elemVolVars[insideScvIdx];
460
461 const Scalar ti = computeTpfaTransmissibility(scvf, insideScv,
462 getPermeability_(problem, insideVolVars, scvf.ipGlobal()),
463 insideVolVars.extrusionFactor());
464
465 // for the boundary (dirichlet) or at branching points we only need ti
466 if (scvf.boundary() || scvf.numOutsideScvs() > 1)
467 tij = scvf.area()*ti;
468
469 // otherwise we compute a tpfa harmonic mean
470 else
471 {
472 const auto outsideScvIdx = scvf.outsideScvIdx();
473 // as we assemble fluxes from the neighbor to our element the outside index
474 // refers to the scv of our element, so we use the scv method
475 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
476 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
477 const Scalar tj = computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv,
478 getPermeability_(problem, outsideVolVars, scvf.ipGlobal()),
479 outsideVolVars.extrusionFactor());
480
481 // harmonic mean (check for division by zero!)
482 // TODO: This could lead to problems!? Is there a better way to do this?
483 if (ti*tj <= 0.0)
484 tij = 0;
485 else
486 tij = scvf.area()*(ti * tj)/(ti + tj);
487 }
488
489 return tij;
490 }
491
492private:
493 template<class Problem, class VolumeVariables,
494 std::enable_if_t<!Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
495 static decltype(auto) getPermeability_(const Problem& problem,
496 const VolumeVariables& volVars,
497 const GlobalPosition& scvfIpGlobal)
498 { return volVars.permeability(); }
499
500 template<class Problem, class VolumeVariables,
501 std::enable_if_t<Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
502 static decltype(auto) getPermeability_(const Problem& problem,
503 const VolumeVariables& volVars,
504 const GlobalPosition& scvfIpGlobal)
505 { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); }
506};
507
508} // end namespace Dumux
509
510#endif
Define some often used mathematical functions.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
forward declaration of the method-specific implementation
Definition: flux/darcyslaw.hh:38
Darcy's law for cell-centered finite volume schemes with two-point flux approximation.
Definition: flux/cctpfa/darcyslaw.hh:49
Declares all properties used in Dumux.
Free functions to evaluate the transmissibilities associated with flux evaluations across sub-control...