3.2-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
velocityadaptive.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_FVVELOCITY2P_ADAPTIVE_HH
25#define DUMUX_FVVELOCITY2P_ADAPTIVE_HH
26
27#include <dune/common/float_cmp.hh>
28#include <dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh>
29
30namespace Dumux {
31
38template<class TypeTag>
39class FVVelocity2PAdaptive: public FVVelocity2P<TypeTag>
40{
41 using ParentType = FVVelocity2P<TypeTag>;
42
43 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
44 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
45 using Problem = GetPropType<TypeTag, Properties::Problem>;
46
47
48 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
49
50 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
51 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
52
53 using CellData = GetPropType<TypeTag, Properties::CellData>;
54
55 using Intersection = typename GridView::Intersection;
56
57 enum
58 {
59 dim = GridView::dimension, dimWorld = GridView::dimensionworld
60 };
61
62 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
63
64 enum
65 {
66 pw = Indices::pressureW,
67 pn = Indices::pressureNw,
68 pGlobal = Indices::pressureGlobal,
69 vw = Indices::velocityW,
70 vn = Indices::velocityNw,
71 vt = Indices::velocityTotal
72 };
73 enum
74 {
75 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, numPhases = getPropValue<TypeTag, Properties::NumPhases>()
76 };
77
78 using Element = typename GridView::template Codim<0>::Entity;
79
80 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
81 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
82 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
83
84public:
89 FVVelocity2PAdaptive(Problem& problem)
90 : ParentType(problem), problem_(problem), gravity_(problem.gravity())
91 {
92 if (getPropValue<TypeTag, Properties::EnableCompressibility>() && velocityType_ == vt)
93 {
94 DUNE_THROW(Dune::NotImplemented, "Total velocity - global pressure - model cannot be used with compressible fluids!");
95 }
96 if (velocityType_ != vw && velocityType_ != vn && velocityType_ != vt)
97 {
98 DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
99 }
100
101 density_[wPhaseIdx] = 0.;
102 density_[nPhaseIdx] = 0.;
103 viscosity_[wPhaseIdx] = 0.;
104 viscosity_[nPhaseIdx] = 0.;
105 }
106
108 void initialize()
109 {
110 ParentType::initialize();
111
112 if (!compressibility_)
113 {
114 const auto element = *problem_.gridView().template begin<0>();
115 FluidState fluidState;
116 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
117 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
118 fluidState.setTemperature(problem_.temperature(element));
119 fluidState.setSaturation(wPhaseIdx, 1.);
120 fluidState.setSaturation(nPhaseIdx, 0.);
121 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
122 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
123 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
124 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
125 }
126 }
127
129 void calculateVelocity(const Intersection& intersection, CellData& cellData);
130
137 bool calculateVelocityInTransport()
138 {
139 return false;
140 }
141
142private:
143 Problem& problem_;
144 const GravityVector& gravity_;
145 Scalar density_[numPhases];
146 Scalar viscosity_[numPhases];
147
149 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
150 static const bool compressibility_ = getPropValue<TypeTag, Properties::EnableCompressibility>();
152 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
154 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
155};
156
164template<class TypeTag>
165void FVVelocity2PAdaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
166{
167 auto elementI = intersection.inside();
168 auto elementJ = intersection.outside();
169
170 if (elementI.level() == elementJ.level())
171 {
172 ParentType::calculateVelocity(intersection, cellData);
173 }
174 else if (elementJ.level() == elementI.level() + 1 && dim == 2)
175 {
176 CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(elementJ));
177
178 // get global coordinates of cell centers
179 const GlobalPosition& globalPosI = elementI.geometry().center();
180 const GlobalPosition& globalPosJ = elementJ.geometry().center();
181
182 // get mobilities and fractional flow factors
183 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
184 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
185 Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx);
186 Scalar fractionalNwI = cellData.fracFlowFunc(nPhaseIdx);
187 Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx);
188 Scalar lambdaNwJ = cellDataJ.mobility(nPhaseIdx);
189 Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx);
190 Scalar fractionalNwJ = cellDataJ.fracFlowFunc(nPhaseIdx);
191
192 // get capillary pressure
193 Scalar pcI = cellData.capillaryPressure();
194 Scalar pcJ = cellDataJ.capillaryPressure();
195
196 // get face index
197 int isIndexI = intersection.indexInInside();
198
199 Scalar faceArea = intersection.geometry().volume();
200 Scalar faceAreaSum = faceArea;
201
202 // get face normal
203 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
204
205 // Count number of hanging nodes
206 // not really necessary
207 int isIndexJ = intersection.indexInOutside();
208
209 int globalIdxK = 0;
210 auto elementK = intersection.outside();
211 // We are looking for two things:
212 // IsIndexJ, the index of the interface from the neighbor-cell point of view
213 // GlobalIdxK, the index of the third cell
214 // for efficiency this is done in one intersection loop
215
216 // Intersectioniterator around cell I
217 for (const auto& intersectionI : intersections(problem_.gridView(), elementI))
218 {
219 if (intersectionI.neighbor())
220 {
221 auto neighbor2 = intersectionI.outside();
222
223 // make sure we do not choose elemntI as third element
224 // -> faces with hanging node have more than one intersection but only one face index!
225 if (neighbor2 != elementJ && intersectionI.indexInInside() == isIndexI)
226 {
227 globalIdxK = problem_.variables().index(neighbor2);
228 elementK = neighbor2;
229 faceAreaSum += intersectionI.geometry().volume();
230
231 break;
232 }
233 }
234 }
235
236 CellData& cellDataK = problem_.variables().cellData(globalIdxK);
237
238 // face global coordinates
239 const GlobalPosition& globalPosInterface = intersection.geometry().center();
240
241 Dune::FieldVector < Scalar, dimWorld > distVec = globalPosInterface - globalPosI;
242 Scalar lI = distVec * unitOuterNormal;
243 distVec = globalPosJ - globalPosInterface;
244 Scalar lJ = distVec * unitOuterNormal;
245 Scalar l = lI + lJ;
246
247 DimMatrix permeabilityI(0);
248 DimMatrix permeabilityJ(0);
249 DimMatrix permeabilityK(0);
250
251 problem_.spatialParams().meanK(permeabilityI,
252 problem_.spatialParams().intrinsicPermeability(elementI));
253 problem_.spatialParams().meanK(permeabilityJ,
254 problem_.spatialParams().intrinsicPermeability(elementJ));
255 problem_.spatialParams().meanK(permeabilityK,
256 problem_.spatialParams().intrinsicPermeability(elementK));
257
258 // Calculate permeablity component normal to interface
259 Scalar kI, kJ, kK, ng, kMean; //, gI, gJ, gK;
260 Dune::FieldVector < Scalar, dim > permI(0);
261 Dune::FieldVector < Scalar, dim > permJ(0);
262 Dune::FieldVector < Scalar, dim > permK(0);
263
264 permeabilityI.mv(unitOuterNormal, permI);
265 permeabilityJ.mv(unitOuterNormal, permJ);
266 permeabilityK.mv(unitOuterNormal, permK);
267
268 // kI,kJ,kK = (n^T)Kn
269 kI = unitOuterNormal * permI;
270 kJ = unitOuterNormal * permJ;
271 kK = unitOuterNormal * permK;
272 kMean = kI * kJ * kK * l / (kJ * kK * lI + kI * (kJ + kK) / 2 * lJ);
273
274 ng = gravity_ * unitOuterNormal;
275
276 Scalar fractionalWK = cellDataK.fracFlowFunc(wPhaseIdx);
277 Scalar fractionalNwK = cellDataK.fracFlowFunc(nPhaseIdx);
278
279 Scalar pcK = cellDataK.capillaryPressure();
280 Scalar pcJK = (pcJ + pcK) / 2;
281
282 // calculate potential gradients
283 // reuse potentials from fvpressure2padaptive
284 Scalar potentialDiffWIJ = cellData.fluxData().upwindPotential(wPhaseIdx, isIndexI);
285 Scalar potentialDiffNwIJ = cellData.fluxData().upwindPotential(nPhaseIdx, isIndexI);
286 Scalar potentialDiffWIK = potentialDiffWIJ;
287 Scalar potentialDiffNwIK = potentialDiffNwIJ;
288 // preliminary potential. The "real" ones are found below
289
290 // Comment: reversed weighting is plausible, too (swap lJ and lI)
291 Scalar rhoMeanWIJ = density_[wPhaseIdx];
292 Scalar rhoMeanNwIJ = density_[nPhaseIdx];
293 Scalar rhoMeanWIK = density_[wPhaseIdx];
294 Scalar rhoMeanNwIK = density_[nPhaseIdx];
295
296 if (compressibility_)
297 {
298 rhoMeanWIJ = (lJ * cellData.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l;
299 rhoMeanNwIJ = (lJ * cellData.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l;
300 rhoMeanWIK = (lJ * cellData.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l;
301 rhoMeanNwIK = (lJ * cellData.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l;
302 }
303
304 Scalar fMeanWIJ = (lJ * fractionalWI + lI * fractionalWJ) / l;
305 Scalar fMeanNwIJ = (lJ * fractionalNwI + lI * fractionalNwJ) / l;
306 Scalar fMeanWIK = (lJ * fractionalWI + lI * fractionalWK) / l;
307 Scalar fMeanNwIK = (lJ * fractionalNwI + lI * fractionalNwK) / l;
308
309 Scalar densityWIJ = density_[wPhaseIdx];
310 Scalar densityNwIJ = density_[nPhaseIdx];
311 Scalar densityWIK = density_[wPhaseIdx];
312 Scalar densityNwIK = density_[nPhaseIdx];
313
314 if (compressibility_)
315 {
316 // Upwinding for finding the upwinding direction
317 densityWIJ = (potentialDiffWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
318 densityNwIJ = (potentialDiffNwIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
319 densityWIK = (potentialDiffWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
320 densityNwIK = (potentialDiffNwIK > 0.) ? cellData.density(nPhaseIdx) : cellDataK.density(nPhaseIdx);
321
322 densityWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? rhoMeanWIJ : densityWIJ;
323 densityNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? rhoMeanNwIJ : densityNwIJ;
324 densityWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? rhoMeanWIK : densityWIK;
325 densityNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? rhoMeanNwIK : densityNwIK;
326 }
327
328 Scalar fractionalWIJ = (potentialDiffWIJ > 0.) ? fractionalWI : fractionalWJ;
329 Scalar fractionalNwIJ = (potentialDiffNwIJ > 0.) ? fractionalNwI : fractionalNwJ;
330 Scalar fractionalWIK = (potentialDiffWIK > 0.) ? fractionalWI : fractionalWK;
331 Scalar fractionalNwIK = (potentialDiffNwIK > 0.) ? fractionalNwI : fractionalNwK;
332
333 fractionalWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? fMeanWIJ : fractionalWIJ;
334 fractionalNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? fMeanNwIJ : fractionalNwIJ;
335 fractionalWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? fMeanWIK : fractionalWIK;
336 fractionalNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? fMeanNwIK : fractionalNwIK;
337
338 switch (pressureType_)
339 {
340 case pGlobal:
341 {
342 Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2;
343
344 potentialDiffWIJ = (cellData.globalPressure() - fMeanNwIJ * pcI - (pressJK - fMeanNwIJ * pcJK)) / l
345 + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
346 potentialDiffNwIJ = (cellData.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l
347 + (densityNwIJ - lJ / l * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2) * ng;
348 potentialDiffWIK = (cellData.globalPressure() - fMeanNwIK * pcI - (pressJK - fMeanNwIK * pcJK)) / l
349 + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
350 potentialDiffNwIK = (cellData.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l
351 + (densityNwIK - lJ / l * (kI + kK) / kI * (densityNwIJ - densityNwIK) / 2) * ng;
352 break;
353 }
354 default:
355 {
356 potentialDiffWIJ = (cellData.pressure(wPhaseIdx)
357 - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
358 + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
359 potentialDiffNwIJ = (cellData.pressure(nPhaseIdx)
360 - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
361 + (densityNwIJ - lJ / l * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2) * ng;
362 potentialDiffWIK = (cellData.pressure(wPhaseIdx)
363 - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
364 + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
365 potentialDiffNwIK = (cellData.pressure(nPhaseIdx)
366 - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
367 + (densityNwIK - lJ / l * (kI + kK) / kI * (densityNwIJ - densityNwIK) / 2) * ng;
368 break;
369 }
370 }
371
372 // store potentials for further calculations (velocity, saturation, ...)
373 // these quantities only have correct sign (needed for upwinding)
374 // potentials are defined slightly different for adaptive scheme
375 cellData.fluxData().addUpwindPotential(wPhaseIdx, isIndexI, potentialDiffWIJ);
376 cellData.fluxData().addUpwindPotential(nPhaseIdx, isIndexI, potentialDiffNwIJ);
377 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, isIndexJ, -potentialDiffWIJ);
378 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, isIndexJ, -potentialDiffNwIJ);
379
380 //do the upwinding of the mobility depending on the phase potentials
381 Scalar lambdaWIJ = (potentialDiffWIJ > 0.) ? lambdaWI : lambdaWJ;
382 lambdaWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
383 Scalar lambdaNwIJ = (potentialDiffNwIJ > 0.) ? lambdaNwI : lambdaNwJ;
384 lambdaNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNwIJ;
385
386 if (compressibility_)
387 {
388 densityWIJ = (potentialDiffWIJ > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
389 densityNwIJ = (potentialDiffNwIJ > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
390 densityWIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIJ, 0.0, 1.0e-30)) ? rhoMeanWIJ : densityWIJ;
391 densityNwIJ = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIJ, 0.0, 1.0e-30)) ? rhoMeanNwIJ : densityNwIJ;
392 densityWIK = (potentialDiffWIK > 0.) ? cellData.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
393 densityNwIK = (potentialDiffNwIK > 0.) ? cellData.density(nPhaseIdx) : cellDataK.density(nPhaseIdx);
394 densityWIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffWIK, 0.0, 1.0e-30)) ? rhoMeanWIK : densityWIK;
395 densityNwIK = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNwIK, 0.0, 1.0e-30)) ? rhoMeanNwIK : densityNwIK;
396 }
397
398 // calculate velocities and the gravity term
399 Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
400 Dune::FieldVector < Scalar, dimWorld > velocityNw(unitOuterNormal);
401 Dune::FieldVector < Scalar, dimWorld > gravityTermW(unitOuterNormal);
402 Dune::FieldVector < Scalar, dimWorld > gravityTermNw(unitOuterNormal);
403
404 gravityTermW *= lambdaWIJ * kMean * ng;
405 gravityTermW *= densityWIJ - (lJ / l) * (kI + kK) / kI * (densityWIK - densityWIJ) / 2;
406 gravityTermNw *= lambdaNwIJ * kMean * ng;
407 gravityTermNw *= densityNwIJ - (lJ / l) * (kI + kK) / kI * (densityNwIK - densityNwIJ) / 2;
408
409 switch (pressureType_)
410 {
411 case pw:
412 {
413 velocityW *= lambdaWIJ * kMean * (cellData.pressure(wPhaseIdx) -
414 (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l;
415 velocityNw *= lambdaNwIJ * kMean * (cellData.pressure(nPhaseIdx) -
416 (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l;
417
418 velocityW += gravityTermW;
419 velocityNw += gravityTermNw;
420 break;
421 }
422 case pn:
423 {
424 velocityNw *= lambdaNwIJ * kMean * (cellData.pressure(nPhaseIdx) -
425 (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l;
426 velocityW *= lambdaWIJ * kMean * (cellData.pressure(wPhaseIdx) -
427 (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l;
428
429 velocityW += gravityTermW;
430 velocityNw += gravityTermNw;
431 break;
432 }
433 case pGlobal:
434 {
435 Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2;
436
437 velocityW *= lambdaWIJ * kMean * (cellData.globalPressure() - pressJK) / l;
438 velocityW += gravityTermW;
439 velocityW += gravityTermNw;
440 velocityNw = 0;
441 break;
442 }
443 }
444
445 cellDataJ.fluxData().setVelocity(wPhaseIdx, isIndexJ, velocityW);
446 cellDataJ.fluxData().setVelocity(nPhaseIdx, isIndexJ, velocityNw);
447 cellDataJ.fluxData().setVelocityMarker(isIndexJ);
448
449 //times 0.5 because cell face with hanging node is called twice! Do not set marker because it should be called twice!
450 velocityW *= faceArea/faceAreaSum;
451 velocityNw *= faceArea/faceAreaSum;
452 cellData.fluxData().addVelocity(wPhaseIdx, isIndexI, velocityW);
453 cellData.fluxData().addVelocity(nPhaseIdx, isIndexI, velocityNw);
454 }
455 else if (elementI.level() > elementJ.level() && dim == 3)
456 {
457 int globalIdxJ = problem_.variables().index(elementJ);
458
459 CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
460
461 // get global coordinates of cell centers
462 const GlobalPosition& globalPosI = elementI.geometry().center();
463 const GlobalPosition& globalPosJ = elementJ.geometry().center();
464
465 // get mobilities and fractional flow factors
466 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
467 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
468 Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx);
469 Scalar lambdaNwJ = cellDataJ.mobility(nPhaseIdx);
470
471 // get capillary pressure
472 Scalar pcI = cellData.capillaryPressure();
473 Scalar pcJ = cellDataJ.capillaryPressure();
474
475 // get face index
476 int isIndexI = intersection.indexInInside();
477 int isIndexJ = intersection.indexInOutside();
478
479 // get face normal
480 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
481
482 // distance vector between barycenters
483 GlobalPosition distVec = globalPosJ - globalPosI;
484
485 // compute distance between cell centers
486 Scalar dist = distVec.two_norm();
487
488 // compute vectorized permeabilities
489 DimMatrix meanPermeability(0);
490
491 problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
492 problem_.spatialParams().intrinsicPermeability(elementJ));
493
494 Dune::FieldVector<Scalar, dim> permeability(0);
495 meanPermeability.mv(unitOuterNormal, permeability);
496
497 // calculate potential gradients
498 Scalar potentialDiffW = cellData.potential(wPhaseIdx) - cellDataJ.potential(wPhaseIdx);
499 Scalar potentialDiffNw = cellData.potential(nPhaseIdx) - cellDataJ.potential(nPhaseIdx);
500
501 if (compressibility_)
502 {
503 density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
504 density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
505
506 density_[wPhaseIdx] =
507 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
508 density_[wPhaseIdx];
509 density_[nPhaseIdx] =
510 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
511 density_[nPhaseIdx];
512
513 potentialDiffW = (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx));
514 potentialDiffNw = (cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx));
515
516
517 potentialDiffW += density_[wPhaseIdx] * (distVec * gravity_); //delta z/delta x in unitOuterNormal[z]
518 potentialDiffNw += density_[nPhaseIdx] * (distVec * gravity_);
519 }
520
521 // store potentials for further calculations (velocity, saturation, ...)
522 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndexI, potentialDiffW);
523 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndexI, potentialDiffNw);
524
525 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, isIndexJ, -potentialDiffW);
526 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, isIndexJ, -potentialDiffNw);
527
528 // do the upwinding of the mobility depending on the phase potentials
529 Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWJ;
530 lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30 )) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
531 Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwJ;
532 lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
533
534 if (compressibility_)
535 {
536 density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
537 density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
538
539 density_[wPhaseIdx] =
540 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
541 density_[wPhaseIdx];
542 density_[nPhaseIdx] =
543 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
544 density_[nPhaseIdx];
545 }
546
547 Scalar scalarPerm = permeability.two_norm();
548
549 // calculate the gravity term
550 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
551 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
552
553 // calculate unit distVec
554 distVec /= dist;
555 Scalar areaScaling = (unitOuterNormal * distVec);
556 // this treatment of g allows to account for gravity flux through faces
557 // where the face normal has no z component (e.g. parallelepiped grids)
558 Scalar gravityTermW = (gravity_ * distVec) * density_[wPhaseIdx] * areaScaling;
559 Scalar gravityTermNw = (gravity_ * distVec) * density_[nPhaseIdx] * areaScaling;
560
561 // calculate velocity depending on the pressure used -> use pc = pn - pw
562 switch (pressureType_)
563 {
564 case pw:
565 {
566 velocityW *= lambdaW * scalarPerm
567 * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermW);
568 velocityNw *= lambdaNw * scalarPerm
569 * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermNw)
570 + 0.5 * (lambdaNwI + lambdaNwJ) * scalarPerm * (pcI - pcJ) / dist;
571 break;
572 }
573 case pn:
574 {
575 velocityW *= lambdaW * scalarPerm
576 * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermW)
577 - 0.5 * (lambdaWI + lambdaWJ) * scalarPerm * (pcI - pcJ) / dist;
578 velocityNw *= lambdaNw * scalarPerm
579 * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermNw);
580 break;
581 }
582 case pGlobal:
583 {
584 velocityW *= (lambdaW + lambdaNw) * scalarPerm * (cellData.globalPressure() - cellDataJ.globalPressure()) / dist
585 + scalarPerm * (lambdaW * gravityTermW + lambdaNw * gravityTermNw);
586 velocityNw = 0;
587 break;
588 }
589 }
590
591 // store velocities
592 cellData.fluxData().setVelocity(wPhaseIdx, isIndexI, velocityW);
593 cellData.fluxData().setVelocity(nPhaseIdx, isIndexI, velocityNw);
594 cellData.fluxData().setVelocityMarker(isIndexI);
595
596 using std::pow;
597 Scalar weightingFactor = pow(0.5, (dim - 1)*(elementI.level() - elementJ.level()));
598
599 velocityW *= weightingFactor;
600 velocityNw *= weightingFactor;
601
602 cellDataJ.fluxData().addVelocity(wPhaseIdx, isIndexJ, velocityW);
603 cellDataJ.fluxData().addVelocity(nPhaseIdx, isIndexJ, velocityNw);
604 cellDataJ.fluxData().setVelocityMarker(isIndexJ);
605 }
606
607 return;
608}
609} // end namespace Dumux
610#endif
Dumux
Definition: adapt.hh:29
Dumux::GetPropType
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:149
Dumux::IOName::viscosity
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
Dumux::IOName::permeability
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
Dumux::IOName::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Dumux::FVVelocity2P
Determines the velocity from a finite volume solution of the pressure equation of a sequential model ...
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:61
Dumux::FVVelocity2P::initialize
void initialize()
For initialization.
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:143
Dumux::FVVelocity2PAdaptive
Determines the velocity from a finite volume solution of the pressure equation of a sequential model ...
Definition: velocityadaptive.hh:40
Dumux::FVVelocity2PAdaptive::initialize
void initialize()
For initialization.
Definition: velocityadaptive.hh:108
Dumux::FVVelocity2PAdaptive::calculateVelocityInTransport
bool calculateVelocityInTransport()
Indicates if velocity is reconstructed in the pressure step or in the transport step.
Definition: velocityadaptive.hh:137
Dumux::FVVelocity2PAdaptive::FVVelocity2PAdaptive
FVVelocity2PAdaptive(Problem &problem)
Constructs a FVVelocity2PAdaptive object.
Definition: velocityadaptive.hh:89
Dumux::FVVelocity2PAdaptive::calculateVelocity
void calculateVelocity(const Intersection &intersection, CellData &cellData)
Calculates the velocity at a cell-cell interface.
Definition: velocityadaptive.hh:165
velocity.hh
Velocity field from a finite volume solution of a pressure equation.
Generated by doxygen 1.9.3