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