3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
cctpfa/forchheimerslaw.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_FORCHHEIMERS_LAW_HH
25#define DUMUX_DISCRETIZATION_CC_TPFA_FORCHHEIMERS_LAW_HH
26
27#include <dune/common/fvector.hh>
28#include <dune/common/fmatrix.hh>
29
30#include <dumux/common/math.hh>
34
37
38namespace Dumux {
39
40// forward declarations
41template<class TypeTag, DiscretizationMethod discMethod>
42class ForchheimersLawImplementation;
43
52template<class Scalar, class GridGeometry, bool isNetwork>
54
60template <class TypeTag>
62: public CCTpfaForchheimersLaw<GetPropType<TypeTag, Properties::Scalar>,
63 GetPropType<TypeTag, Properties::GridGeometry>,
64 (GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension < GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld)>
65{};
66
71template<class GridGeometry>
72class TpfaForchheimersLawCacheFiller
73{
74 using FVElementGeometry = typename GridGeometry::LocalView;
75 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
76 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
77
78public:
82 template<class FluxVariablesCache, class Problem, class ElementVolumeVariables, class FluxVariablesCacheFiller>
83 static void fill(FluxVariablesCache& scvfFluxVarsCache,
84 const Problem& problem,
85 const Element& element,
86 const FVElementGeometry& fvGeometry,
87 const ElementVolumeVariables& elemVolVars,
88 const SubControlVolumeFace& scvf,
89 const FluxVariablesCacheFiller& fluxVarsCacheFiller)
90 {
91 scvfFluxVarsCache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf);
92 }
93};
94
99template<class AdvectionType, class GridGeometry>
100class TpfaForchheimersLawCache
101{
102 using Scalar = typename AdvectionType::Scalar;
103 using FVElementGeometry = typename GridGeometry::LocalView;
104 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
105 using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
106 static constexpr int dimWorld = GridGeometry::GridView::dimensionworld;
107 using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
108
109public:
110 using Filler = TpfaForchheimersLawCacheFiller<GridGeometry>;
111
112 template<class Problem, class ElementVolumeVariables>
113 void updateAdvection(const Problem& problem,
114 const Element& element,
115 const FVElementGeometry& fvGeometry,
116 const ElementVolumeVariables& elemVolVars,
117 const SubControlVolumeFace &scvf)
118 {
119 tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf);
120 harmonicMeanSqrtK_ = AdvectionType::calculateHarmonicMeanSqrtPermeability(problem, elemVolVars, scvf);
121 }
122
123 const Scalar& advectionTij() const
124 { return tij_; }
125
126 const DimWorldMatrix& harmonicMeanSqrtPermeability() const
127 { return harmonicMeanSqrtK_; }
128
129private:
130 Scalar tij_;
131 DimWorldMatrix harmonicMeanSqrtK_;
132};
133
138template<class ScalarType, class GridGeometry>
139class CCTpfaForchheimersLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
140{
141 using ThisType = CCTpfaForchheimersLaw<ScalarType, GridGeometry, /*isNetwork*/ false>;
142 using FVElementGeometry = typename GridGeometry::LocalView;
143 using SubControlVolume = typename GridGeometry::SubControlVolume;
144 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
145 using GridView = typename GridGeometry::GridView;
146 using Element = typename GridView::template Codim<0>::Entity;
147
148 static constexpr int dim = GridView::dimension;
149 static constexpr int dimWorld = GridView::dimensionworld;
150
151 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
152 using DimWorldVector = Dune::FieldVector<ScalarType, dimWorld>;
153 using DimWorldMatrix = Dune::FieldMatrix<ScalarType, dimWorld, dimWorld>;
154
155 using DarcysLaw = CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>;
156
157 public:
159 using Scalar = ScalarType;
160
162 static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
163
165 using Cache = TpfaForchheimersLawCache<ThisType, GridGeometry>;
166
224 template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
225 static Scalar flux(const Problem& problem,
226 const Element& element,
227 const FVElementGeometry& fvGeometry,
228 const ElementVolumeVariables& elemVolVars,
229 const SubControlVolumeFace& scvf,
230 int phaseIdx,
231 const ElementFluxVarsCache& elemFluxVarsCache)
232 {
233 const auto velocity = velocity_(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
234 Scalar flux = velocity * scvf.unitOuterNormal();
235 flux *= scvf.area();
236 return flux;
237 }
238
241 template<class Problem, class ElementVolumeVariables>
242 static Scalar calculateTransmissibility(const Problem& problem,
243 const Element& element,
244 const FVElementGeometry& fvGeometry,
245 const ElementVolumeVariables& elemVolVars,
246 const SubControlVolumeFace& scvf)
247 {
248 return DarcysLaw::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf);
249 }
250
255 template<class Problem, class ElementVolumeVariables,
256 bool scalarPerm = std::is_same<typename Problem::SpatialParams::PermeabilityType, Scalar>::value,
257 std::enable_if_t<scalarPerm, int> = 0>
258 static DimWorldMatrix calculateHarmonicMeanSqrtPermeability(const Problem& problem,
259 const ElementVolumeVariables& elemVolVars,
260 const SubControlVolumeFace& scvf)
261 {
262 using std::sqrt;
263 Scalar harmonicMeanSqrtK(0.0);
264
265 const auto insideScvIdx = scvf.insideScvIdx();
266 const auto& insideVolVars = elemVolVars[insideScvIdx];
267 const Scalar Ki = getPermeability_(problem, insideVolVars, scvf.ipGlobal());
268 const Scalar sqrtKi = sqrt(Ki);
269
270 if (!scvf.boundary())
271 {
272 const auto outsideScvIdx = scvf.outsideScvIdx();
273 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
274 const Scalar Kj = getPermeability_(problem, outsideVolVars, scvf.ipGlobal());
275 const Scalar sqrtKj = sqrt(Kj);
276 harmonicMeanSqrtK = problem.spatialParams().harmonicMean(sqrtKi, sqrtKj, scvf.unitOuterNormal());
277 }
278 else
279 harmonicMeanSqrtK = sqrtKi;
280
281 // Convert to tensor
282 DimWorldMatrix result(0.0);
283 for (int i = 0; i < dimWorld; ++i)
284 result[i][i] = harmonicMeanSqrtK;
285
286 return result;
287 }
288
293 template<class Problem, class ElementVolumeVariables,
294 bool scalarPerm = std::is_same<typename Problem::SpatialParams::PermeabilityType, Scalar>::value,
295 std::enable_if_t<!scalarPerm, int> = 0>
296 static DimWorldMatrix calculateHarmonicMeanSqrtPermeability(const Problem& problem,
297 const ElementVolumeVariables& elemVolVars,
298 const SubControlVolumeFace& scvf)
299 {
300 using std::sqrt;
301 DimWorldMatrix sqrtKi(0.0);
302 DimWorldMatrix sqrtKj(0.0);
303 DimWorldMatrix harmonicMeanSqrtK(0.0);
304
305 const auto insideScvIdx = scvf.insideScvIdx();
306 const auto& insideVolVars = elemVolVars[insideScvIdx];
307 const auto& Ki = getPermeability_(problem, insideVolVars, scvf.ipGlobal());
308 // Make sure the permeability matrix does not have off-diagonal entries
309 // TODO implement method using eigenvalues and eigenvectors for general tensors
310 if (!isDiagonal_(Ki))
311 DUNE_THROW(Dune::NotImplemented, "Only diagonal permeability tensors are supported.");
312
313 for (int i = 0; i < dim; ++i)
314 sqrtKi[i][i] = sqrt(Ki[i][i]);
315
316 if (!scvf.boundary())
317 {
318 const auto outsideScvIdx = scvf.outsideScvIdx();
319 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
320 const auto& Kj = getPermeability_(problem, outsideVolVars, scvf.ipGlobal());
321 // Make sure the permeability matrix does not have off-diagonal entries
322 if (!isDiagonal_(Kj))
323 DUNE_THROW(Dune::NotImplemented, "Only diagonal permeability tensors are supported.");
324
325 for (int i = 0; i < dim; ++i)
326 sqrtKj[i][i] = sqrt(Kj[i][i]);
327
328 harmonicMeanSqrtK = problem.spatialParams().harmonicMean(sqrtKi, sqrtKj, scvf.unitOuterNormal());
329 }
330 else
331 harmonicMeanSqrtK = sqrtKi;
332
333 return harmonicMeanSqrtK;
334 }
335
336private:
337
339 template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
340 static DimWorldVector velocity_(const Problem& problem,
341 const Element& element,
342 const FVElementGeometry& fvGeometry,
343 const ElementVolumeVariables& elemVolVars,
344 const SubControlVolumeFace& scvf,
345 int phaseIdx,
346 const ElementFluxVarsCache& elemFluxVarsCache)
347 {
348 // Get the volume flux based on Darcy's law. The value returned by this method needs to be multiplied with the
349 // mobility (upwinding).
350 Scalar darcyFlux = DarcysLaw::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
351 auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.mobility(phaseIdx); };
352 const Scalar darcyUpwindMobility = doUpwind_(scvf, elemVolVars, upwindTerm, !std::signbit(darcyFlux) /*insideIsUpstream*/);
353 darcyFlux *= darcyUpwindMobility;
354
355 // Convert to Darcy velocity
356 DimWorldVector darcyVelocity = scvf.unitOuterNormal();
357 darcyVelocity *= darcyFlux / scvf.area();
358
359 // Get the harmonic mean of the square root of permeability
360 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
361 const auto& sqrtK = fluxVarsCache.harmonicMeanSqrtPermeability();
362
363 // Obtain the Forchheimer coefficient from the spatial parameters
364 const Scalar forchCoeff = problem.spatialParams().forchCoeff(scvf);
365
366 // Initial guess of velocity: Darcy relation
367 DimWorldVector velocity = darcyVelocity;
368
369 DimWorldVector deltaV(0.0); // the change in velocity between Newton iterations
370 DimWorldVector residual(10e10); // the residual (function value that is to be minimized)
371 DimWorldMatrix gradF(0.0); // slope of equation that is to be solved
372
373 // Search by means of the Newton method for a root of Forchheimer equation
374 static const Scalar epsilon = getParamFromGroup<Scalar>(problem.paramGroup(), "Forchheimer.NewtonTolerance", 1e-12);
375 static const std::size_t maxNumIter = getParamFromGroup<std::size_t>(problem.paramGroup(), "Forchheimer.MaxIterations", 30);
376 for (int k = 0; residual.two_norm() > epsilon ; ++k)
377 {
378 if (k >= maxNumIter)
379 DUNE_THROW(NumericalProblem, "could not determine forchheimer velocity within "
380 << k << " iterations");
381
382 // calculate current value of Forchheimer equation
383 forchheimerResidual_(residual,
384 forchCoeff,
385 sqrtK,
386 darcyVelocity,
387 velocity,
388 elemVolVars,
389 scvf,
390 phaseIdx);
391
392 // newton's method: slope (gradF) and current value (residual) of function is known,
393 forchheimerDerivative_(gradF,
394 forchCoeff,
395 sqrtK,
396 velocity,
397 elemVolVars,
398 scvf,
399 phaseIdx);
400
401 // solve for change in velocity ("x-Axis intercept")
402 gradF.solve(deltaV, residual);
403 velocity -= deltaV;
404 }
405
406 // The fluxvariables expect a value on which an upwinding of the mobility is performed.
407 // We therefore divide by the mobility here.
408 const Scalar forchheimerUpwindMobility = doUpwind_(scvf, elemVolVars, upwindTerm,
409 !std::signbit(velocity * scvf.unitOuterNormal()) /*insideIsUpstream*/);
410
411 // Do not divide by zero. If the mobility is zero, the resulting flux with upwinding will be zero anyway.
412 if (forchheimerUpwindMobility > 0.0)
413 velocity /= forchheimerUpwindMobility;
414
415 return velocity;
416 }
417
419 template<class ElementVolumeVariables>
420 static void forchheimerResidual_(DimWorldVector& residual,
421 const Scalar forchCoeff,
422 const DimWorldMatrix& sqrtK,
423 const DimWorldVector& darcyVelocity,
424 const DimWorldVector& oldVelocity,
425 const ElementVolumeVariables& elemVolVars,
426 const SubControlVolumeFace& scvf,
427 const int phaseIdx)
428 {
429 residual = oldVelocity;
430
431 // residual += k_r/mu K grad p
432 residual -= darcyVelocity;
433
434 // residual += c_F rho / mu abs(v) sqrt(K) v
435 auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx) / volVars.viscosity(phaseIdx); };
436 bool insideIsUpstream = !std::signbit(oldVelocity * scvf.unitOuterNormal());
437 const Scalar rhoOverMu = doUpwind_(scvf, elemVolVars, upwindTerm, insideIsUpstream);
438 sqrtK.usmv(forchCoeff * rhoOverMu * oldVelocity.two_norm(), oldVelocity, residual);
439 }
440
442 template<class ElementVolumeVariables>
443 static void forchheimerDerivative_(DimWorldMatrix& derivative,
444 const Scalar forchCoeff,
445 const DimWorldMatrix& sqrtK,
446 const DimWorldVector& velocity,
447 const ElementVolumeVariables& elemVolVars,
448 const SubControlVolumeFace& scvf,
449 const int phaseIdx)
450 {
451
452
453 // Initialize because for low velocities, we add and do not set the entries.
454 derivative = 0.0;
455
456 auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx) / volVars.viscosity(phaseIdx); };
457 bool insideIsUpstream = !std::signbit(velocity * scvf.unitOuterNormal());
458
459 const Scalar absV = velocity.two_norm() ;
460 // This part of the derivative is only calculated if the velocity is sufficiently small:
461 // do not divide by zero!
462 // This should not be very bad because the derivative is only intended to give an
463 // approximation of the gradient of the
464 // function that goes into the Newton scheme.
465 // This is important in the case of a one-phase region in two-phase flow. The non-existing
466 // phase has a velocity of zero (kr=0).
467 // f = sqrtK * vTimesV (this is a matrix)
468 // f *= forchCoeff density / |velocity| / viscosity (multiply all entries with scalars)
469 const Scalar rhoOverMu = doUpwind_(scvf, elemVolVars, upwindTerm, insideIsUpstream);
470 if (absV > 1e-20)
471 {
472 for (int i = 0; i < dim; i++)
473 {
474 for (int k = 0; k <= i; k++)
475 {
476 derivative[i][k] = sqrtK[i][i] * velocity[i] * velocity[k] * forchCoeff
477 / absV * rhoOverMu;
478
479 if (k < i)
480 derivative[k][i] = derivative[i][k];
481 }
482 }
483 }
484
485 // add on the main diagonal:
486 // 1 + sqrtK_i forchCoeff density |v| / viscosity
487 for (int i = 0; i < dim; i++)
488 derivative[i][i] += (1.0 + (sqrtK[i][i]*forchCoeff * absV * rhoOverMu));
489 }
490
491 template<class ElementVolumeVariables, class UpwindTermFunction>
492 static Scalar doUpwind_(const SubControlVolumeFace& scvf,
493 const ElementVolumeVariables& elemVolVars,
494 const UpwindTermFunction& upwindTerm,
495 const bool insideIsUpstream)
496 {
497 static const Scalar upwindWeight = getParam<Scalar>("Flux.UpwindWeight");
498
499 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
500 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
501
502 if (!insideIsUpstream) // if sign of flux is negative
503 return (upwindWeight*upwindTerm(outsideVolVars)
504 + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
505 else
506 return (upwindWeight*upwindTerm(insideVolVars)
507 + (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
508 }
509
518 static bool isDiagonal_(const DimWorldMatrix & K)
519 {
520 using std::abs;
521 for (int i = 0; i < dim; i++) {
522 for (int k = 0; k < dim; k++) {
523 if ((i != k) && (abs(K[i][k]) >= 1e-25)) {
524 return false;
525 }
526 }
527 }
528 return true;
529 }
530
531 template<class Problem, class VolumeVariables,
532 std::enable_if_t<!Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
533 static decltype(auto) getPermeability_(const Problem& problem,
534 const VolumeVariables& volVars,
535 const GlobalPosition& scvfIpGlobal)
536 { return volVars.permeability(); }
537
538 template<class Problem, class VolumeVariables,
539 std::enable_if_t<Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
540 static decltype(auto) getPermeability_(const Problem& problem,
541 const VolumeVariables& volVars,
542 const GlobalPosition& scvfIpGlobal)
543 { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); }
544};
545
550template<class ScalarType, class GridGeometry>
551class CCTpfaForchheimersLaw<ScalarType, GridGeometry, /*isNetwork*/ true>
552{
553 static_assert(AlwaysFalse<ScalarType>::value, "Forchheimer not implemented for network grids");
554};
555
556} // end namespace Dumux
557
558#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 declare
Definition: forchheimerslaw.hh:38
Forchheimer's law for cell-centered finite volume schemes with two-point flux approximation.
Definition: cctpfa/forchheimerslaw.hh:53
Declares all properties used in Dumux.
Darcy's law for cell-centered finite volume schemes with two-point flux approximation.