3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
2p/sequential/diffusion/cellcentered/velocity.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_HH
25#define DUMUX_FVVELOCITY2P_HH
26
27#include <dune/common/float_cmp.hh>
28#include <dune/grid/common/gridenums.hh>
30
31namespace Dumux {
32
59template<class TypeTag>
61{
62 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
63 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
64 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
65
66 using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
67 using MaterialLaw = typename SpatialParams::MaterialLaw;
68
69 using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
70
71 using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
72 using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
73
74 using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
75 using SolutionTypes = typename GET_PROP(TypeTag, SolutionTypes);
76 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
77 using CellData = typename GET_PROP_TYPE(TypeTag, CellData);
78
79 using Element = typename GridView::Traits::template Codim<0>::Entity;
80 using Intersection = typename GridView::Intersection;
81
82 using Geometry = typename Element::Geometry;
83 using JacobianTransposed = typename Geometry::JacobianTransposed;
84
85 enum
86 {
87 dim = GridView::dimension, dimWorld = GridView::dimensionworld
88 };
89
90 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
91
92 enum
93 {
94 pw = Indices::pressureW,
95 pn = Indices::pressureNw,
96 pGlobal = Indices::pressureGlobal,
97 vw = Indices::velocityW,
98 vn = Indices::velocityNw,
99 vt = Indices::velocityTotal,
100 sw = Indices::saturationW,
101 sn = Indices::saturationNw,
102 pressureIdx = Indices::pressureIdx,
103 saturationIdx = Indices::saturationIdx,
104 eqIdxPress = Indices::pressureEqIdx,
105 eqIdxSat = Indices::satEqIdx
106 };
107 enum
108 {
109 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
110 };
111
112 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
113 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
114 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
115
116public:
121 FVVelocity2P(Problem& problem) :
122 problem_(problem), gravity_(problem.gravity())
123 {
124 if (GET_PROP_VALUE(TypeTag, EnableCompressibility) && velocityType_ == vt)
125 {
126 DUNE_THROW(Dune::NotImplemented,
127 "Total velocity - global pressure - model cannot be used with compressible fluids!");
128 }
129 if (velocityType_ != vw && velocityType_ != vn && velocityType_ != vt)
130 {
131 DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
132 }
133
134 density_[wPhaseIdx] = 0.;
135 density_[nPhaseIdx] = 0.;
136 viscosity_[wPhaseIdx] = 0.;
137 viscosity_[nPhaseIdx] = 0.;
138
139 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
140 }
141
144 {
145 if (!compressibility_)
146 {
147 const auto element = *problem_.gridView().template begin<0> ();
148 FluidState fluidState;
149 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
150 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
151 fluidState.setTemperature(problem_.temperature(element));
152 fluidState.setSaturation(wPhaseIdx, 1.);
153 fluidState.setSaturation(nPhaseIdx, 0.);
154 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
155 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
156 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
157 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
158 }
159 }
160
162 void calculateVelocity(const Intersection&, CellData&);
163
165 void calculateVelocityOnBoundary(const Intersection&, CellData&);
166
173 {
174 return true;
175 }
176
185 template<class MultiWriter>
186 void addOutputVtkFields(MultiWriter &writer)
187 {
188 if (vtkOutputLevel_ > 0)
189 {
190 Dune::BlockVector < Dune::FieldVector<Scalar, dim> > &velocity = *(writer.template allocateManagedBuffer<Scalar,
191 dim>(problem_.gridView().size(0)));
192 Dune::BlockVector < Dune::FieldVector<Scalar, dim> > &velocitySecondPhase =
193 *(writer.template allocateManagedBuffer<Scalar, dim>(problem_.gridView().size(0)));
194
195 // compute update vector
196 for (const auto& element : elements(problem_.gridView()))
197 {
198 // cell index
199 int eIdxGlobal = problem_.variables().index(element);
200
201 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
202
203 const typename Element::Geometry& geometry = element.geometry();
204
205 // get corresponding reference element
206 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
207 const auto refElement = ReferenceElements::general(geometry.type());
208
209 const int numberOfFaces=refElement.size(1);
210 std::vector<Scalar> fluxW(numberOfFaces,0);
211 std::vector<Scalar> fluxNw(numberOfFaces,0);
212
213 // run through all intersections with neighbors and boundary
214 for (const auto& intersection : intersections(problem_.gridView(), element))
215 {
216 int isIndex = intersection.indexInInside();
217
218 fluxW[isIndex] += intersection.geometry().volume()
219 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
220 fluxNw[isIndex] += intersection.geometry().volume()
221 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
222 }
223
224 // calculate velocity on reference element as the Raviart-Thomas-0
225 // interpolant of the fluxes
226 Dune::FieldVector<Scalar, dim> refVelocity;
227 // simplices
228 if (refElement.type().isSimplex()) {
229 for (int dimIdx = 0; dimIdx < dim; dimIdx++)
230 {
231 refVelocity[dimIdx] = -fluxW[dim - 1 - dimIdx];
232 for (int fIdx = 0; fIdx < dim + 1; fIdx++)
233 {
234 refVelocity[dimIdx] += fluxW[fIdx]/(dim + 1);
235 }
236 }
237 }
238 // cubes
239 else if (refElement.type().isCube()){
240 for (int i = 0; i < dim; i++)
241 refVelocity[i] = 0.5 * (fluxW[2*i + 1] - fluxW[2*i]);
242 }
243 // 3D prism and pyramids
244 else {
245 DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
246 }
247
248 const Dune::FieldVector<Scalar, dim> localPos =
249 refElement.position(0, 0);
250
251 // get the transposed Jacobian of the element mapping
252 const JacobianTransposed jacobianT =
253 geometry.jacobianTransposed(localPos);
254
255 // calculate the element velocity by the Piola transformation
256 Dune::FieldVector < Scalar, dim > elementVelocity(0);
257 jacobianT.umtv(refVelocity, elementVelocity);
258 elementVelocity /= geometry.integrationElement(localPos);
259
260 velocity[eIdxGlobal] = elementVelocity;
261
262 // calculate velocity on reference element as the Raviart-Thomas-0
263 // interpolant of the fluxes
264 // simplices
265 if (refElement.type().isSimplex()) {
266 for (int dimIdx = 0; dimIdx < dim; dimIdx++)
267 {
268 refVelocity[dimIdx] = -fluxNw[dim - 1 - dimIdx];
269 for (int fIdx = 0; fIdx < dim + 1; fIdx++)
270 {
271 refVelocity[dimIdx] += fluxNw[fIdx]/(dim + 1);
272 }
273 }
274 }
275 // cubes
276 else if (refElement.type().isCube()){
277 for (int i = 0; i < dim; i++)
278 refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
279 }
280 // 3D prism and pyramids
281 else {
282 DUNE_THROW(Dune::NotImplemented, "velocity output for prism/pyramid not implemented");
283 }
284
285 // calculate the element velocity by the Piola transformation
286 elementVelocity = 0;
287 jacobianT.umtv(refVelocity, elementVelocity);
288 elementVelocity /= geometry.integrationElement(localPos);
289
290 velocitySecondPhase[eIdxGlobal] = elementVelocity;
291 }
292
293 // switch velocities
294 if (velocityType_ == vt)
295 {
296 writer.attachCellData(velocity, "total velocity", dim);
297 }
298 else
299 {
300 writer.attachCellData(velocity, "wetting-velocity", dim);
301 writer.attachCellData(velocitySecondPhase, "non-wetting-velocity", dim);
302 }
303 }
304
305 return;
306 }
307
308private:
309 Problem& problem_;
310 const GravityVector& gravity_;
311 Scalar density_[numPhases];
312 Scalar viscosity_[numPhases];
313
314 int vtkOutputLevel_;
315
317 static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation);
318 static const bool compressibility_ = GET_PROP_VALUE(TypeTag, EnableCompressibility);
320 static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation);
322 static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
323};
324
333template<class TypeTag>
334void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
335{
336 auto elementI = intersection.inside();
337 auto elementJ = intersection.outside();
338
339 int eIdxGlobalJ = problem_.variables().index(elementJ);
340
341 CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
342
343 // get global coordinates of cell centers
344 const GlobalPosition& globalPosI = (elementI).geometry().center();
345 const GlobalPosition& globalPosJ = (elementJ).geometry().center();
346
347 // get mobilities and fractional flow factors
348 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
349 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
350 Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx);
351 Scalar lambdaNwJ = cellDataJ.mobility(nPhaseIdx);
352
353 // get capillary pressure
354 Scalar pcI = cellData.capillaryPressure();
355 Scalar pcJ = cellDataJ.capillaryPressure();
356
357 // get face index
358 int isIndexI = intersection.indexInInside();
359 int isIndexJ = intersection.indexInOutside();
360
361 // get face normal
362 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
363
364 // distance vector between barycenters
365 GlobalPosition distVec = globalPosJ - globalPosI;
366
367 // compute distance between cell centers
368 Scalar dist = distVec.two_norm();
369
370 // compute vectorized permeabilities
371 DimMatrix meanPermeability(0);
372
373 problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(elementI),
374 problem_.spatialParams().intrinsicPermeability(elementJ));
375
376 Dune::FieldVector<Scalar, dim> permeability(0);
377 meanPermeability.mv(unitOuterNormal, permeability);
378
379 // calculate potential gradients
380 Scalar potentialDiffW = cellData.potential(wPhaseIdx) - cellDataJ.potential(wPhaseIdx);
381 Scalar potentialDiffNw = cellData.potential(nPhaseIdx) - cellDataJ.potential(nPhaseIdx);
382
383 if (compressibility_)
384 {
385 density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
386 density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
387
388 density_[wPhaseIdx] =
389 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
390 density_[wPhaseIdx];
391 density_[nPhaseIdx] =
392 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
393 density_[nPhaseIdx];
394
395 potentialDiffW = (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx));
396 potentialDiffNw = (cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx));
397
398 potentialDiffW += density_[wPhaseIdx] * (distVec * gravity_); //delta z/delta x in unitOuterNormal[z]
399 potentialDiffNw += density_[nPhaseIdx] * (distVec * gravity_);
400 }
401
402 // store potentials for further calculations (velocity, saturation, ...)
403 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndexI, potentialDiffW);
404 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndexI, potentialDiffNw);
405
406 cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, isIndexJ, -potentialDiffW);
407 cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, isIndexJ, -potentialDiffNw);
408
409 // do the upwinding of the mobility depending on the phase potentials
410 Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWJ;
411 lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
412 Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwJ;
413 lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwJ) : lambdaNw;
414
415 if (compressibility_)
416 {
417 density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
418 density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
419
420 density_[wPhaseIdx] =
421 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + cellDataJ.density(wPhaseIdx)) :
422 density_[wPhaseIdx];
423 density_[nPhaseIdx] =
424 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + cellDataJ.density(nPhaseIdx)) :
425 density_[nPhaseIdx];
426 }
427
428 Scalar scalarPerm = permeability.two_norm();
429
430 // calculate the gravity term
431 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
432 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
433
434 // calculate unit distVec
435 distVec /= dist;
436 Scalar areaScaling = (unitOuterNormal * distVec);
437 // this treatment of g allows to account for gravity flux through faces where the face normal
438 // has no z component (e.g. parallelepiped grids)
439 Scalar gravityTermW = (gravity_ * distVec) * density_[wPhaseIdx] * areaScaling;
440 Scalar gravityTermNw = (gravity_ * distVec) * density_[nPhaseIdx] * areaScaling;
441
442 // calculate velocity depending on the pressure used -> use pc = pn - pw
443 switch (pressureType_)
444 {
445 case pw:
446 {
447 velocityW *= lambdaW * scalarPerm
448 * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermW);
449 velocityNw *= lambdaNw * scalarPerm
450 * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermNw)
451 + 0.5 * (lambdaNwI + lambdaNwJ) * scalarPerm * (pcI - pcJ) / dist;
452 break;
453 }
454 case pn:
455 {
456 velocityW *= lambdaW * scalarPerm
457 * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermW)
458 - 0.5 * (lambdaWI + lambdaWJ) * scalarPerm * (pcI - pcJ) / dist;
459 velocityNw *= lambdaNw * scalarPerm
460 * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermNw);
461 break;
462 }
463 case pGlobal:
464 {
465 velocityW *= (lambdaW + lambdaNw) * scalarPerm * (cellData.globalPressure() - cellDataJ.globalPressure()) / dist
466 + scalarPerm * (lambdaW * gravityTermW + lambdaNw * gravityTermNw);
467 velocityNw = 0;
468 break;
469 }
470 }
471
472 // store velocities
473 cellData.fluxData().setVelocity(wPhaseIdx, isIndexI, velocityW);
474 cellData.fluxData().setVelocity(nPhaseIdx, isIndexI, velocityNw);
475 cellData.fluxData().setVelocityMarker(isIndexI);
476
477 cellDataJ.fluxData().setVelocity(wPhaseIdx, isIndexJ, velocityW);
478 cellDataJ.fluxData().setVelocity(nPhaseIdx, isIndexJ, velocityNw);
479 cellDataJ.fluxData().setVelocityMarker(isIndexJ);
480
481// printvector(std::cout, cellData.fluxData().velocity(), "velocity", "row", 4, 1, 3);
482 return;
483}
484
493template<class TypeTag>
494void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
495{
496 auto element = intersection.inside();
497
498 // get face index
499 int isIndex = intersection.indexInInside();
500
501 // get face normal
502 const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
503
504 BoundaryTypes bcType;
505 //get boundary type
506 problem_.boundaryTypes(bcType, intersection);
507 PrimaryVariables boundValues(0.0);
508
509 if (bcType.isDirichlet(eqIdxPress))
510 {
511 problem_.dirichlet(boundValues, intersection);
512
513 // get global coordinates of cell centers
514 const GlobalPosition& globalPosI = element.geometry().center();
515
516 // center of face in global coordinates
517 const GlobalPosition& globalPosJ = intersection.geometry().center();
518
519 // get mobilities and fractional flow factors
520 Scalar lambdaWI = cellData.mobility(wPhaseIdx);
521 Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
522 Scalar fractionalWI = cellData.fracFlowFunc(wPhaseIdx);
523 Scalar fractionalNwI = cellData.fracFlowFunc(nPhaseIdx);
524
525 // get capillary pressure
526 Scalar pcI = cellData.capillaryPressure();
527
528 // distance vector between barycenters
529 GlobalPosition distVec = globalPosJ - globalPosI;
530
531 // compute distance between cell centers
532 Scalar dist = distVec.two_norm();
533
534 // permeability vector at boundary
535 // compute vectorized permeabilities
536 DimMatrix meanPermeability(0);
537
538 problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(element));
539
540 Dune::FieldVector<Scalar, dim> permeability(0);
541 meanPermeability.mv(unitOuterNormal, permeability);
542
543 // determine saturation at the boundary -> if no saturation is known directly at the boundary use the cell saturation
544 Scalar satW = 0;
545 Scalar satNw = 0;
546 if (bcType.isDirichlet(eqIdxSat))
547 {
548 switch (saturationType_)
549 {
550 case sw:
551 {
552 satW = boundValues[saturationIdx];
553 satNw = 1 - boundValues[saturationIdx];
554 break;
555 }
556 case sn:
557 {
558 satW = 1 - boundValues[saturationIdx];
559 satNw = boundValues[saturationIdx];
560 break;
561 }
562 }
563 }
564 else
565 {
566 satW = cellData.saturation(wPhaseIdx);
567 satNw = cellData.saturation(nPhaseIdx);
568 }
569
570 Scalar pressBound = boundValues[pressureIdx];
571 Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
572
573 // determine phase pressures from primary pressure variable
574 Scalar pressWBound = 0;
575 Scalar pressNwBound = 0;
576 if (pressureType_ == pw)
577 {
578 pressWBound = pressBound;
579 pressNwBound = pressBound + pcBound;
580 }
581 else if (pressureType_ == pn)
582 {
583 pressWBound = pressBound - pcBound;
584 pressNwBound = pressBound;
585 }
586
587 // get temperature at current position
588 Scalar temperature = problem_.temperature(element);
589
590 Scalar densityWBound = density_[wPhaseIdx];
591 Scalar densityNwBound = density_[nPhaseIdx];
592 Scalar viscosityWBound = viscosity_[wPhaseIdx];
593 Scalar viscosityNwBound = viscosity_[nPhaseIdx];
594
595 if (compressibility_)
596 {
597 FluidState fluidState;
598 fluidState.setSaturation(wPhaseIdx, satW);
599 fluidState.setSaturation(nPhaseIdx, satNw);
600 fluidState.setTemperature(temperature);
601 fluidState.setPressure(wPhaseIdx, pressWBound);
602 fluidState.setPressure(nPhaseIdx, pressNwBound);
603
604 densityWBound = FluidSystem::density(fluidState, wPhaseIdx);
605 densityNwBound = FluidSystem::density(fluidState, nPhaseIdx);
606 viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx) / densityWBound;
607 viscosityNwBound = FluidSystem::viscosity(fluidState, nPhaseIdx) / densityNwBound;
608 }
609
610 Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
611 / viscosityWBound;
612 Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
613 / viscosityNwBound;
614
615 Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
616 Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
617
618 if (compressibility_)
619 {
620 density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : densityWBound;
621 density_[wPhaseIdx] =
622 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
623 density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : densityNwBound;
624 density_[nPhaseIdx] =
625 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + densityNwBound) : density_[nPhaseIdx];
626 }
627
628 // calculate potential gradient
629 if (pressureType_ == pGlobal)
630 {
631 potentialDiffW = (cellData.globalPressure() - pressBound - fractionalNwI * (pcI - pcBound));
632 potentialDiffNw = (cellData.globalPressure() - pressBound + fractionalWI * (pcI - pcBound));
633 }
634 else
635 {
636 potentialDiffW = (cellData.pressure(wPhaseIdx) - pressWBound);
637 potentialDiffNw = (cellData.pressure(nPhaseIdx) - pressNwBound);
638 }
639
640 potentialDiffW += density_[wPhaseIdx] * (distVec * gravity_);
641 potentialDiffNw += density_[nPhaseIdx] * (distVec * gravity_);
642
643 // store potential gradients for further calculations
644 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, potentialDiffW);
645 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, potentialDiffNw);
646
647 // do the upwinding of the mobility depending on the phase potentials
648 Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
649 lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
650 Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
651 lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
652
653 if (compressibility_)
654 {
655 density_[wPhaseIdx] = (potentialDiffW > 0.) ? cellData.density(wPhaseIdx) : densityWBound;
656 density_[wPhaseIdx] =
657 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
658 density_[nPhaseIdx] = (potentialDiffNw > 0.) ? cellData.density(nPhaseIdx) : densityNwBound;
659 density_[nPhaseIdx] =
660 (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + densityNwBound) : density_[nPhaseIdx];
661 }
662
663 Scalar scalarPerm = permeability.two_norm();
664
665 // calculate the gravity term
666 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
667 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
668
669 // calculate unit distVec
670 distVec /= dist;
671 Scalar areaScaling = (unitOuterNormal * distVec);
672 // this treatment of g allows to account for gravity flux through faces where the face normal
673 // has no z component (e.g. parallelepiped grids)
674 Scalar gravityTermW = (gravity_ * distVec) * density_[wPhaseIdx] * areaScaling;
675 Scalar gravityTermNw = (gravity_ * distVec) * density_[nPhaseIdx] * areaScaling;
676
677 // calculate velocity depending on the pressure used -> use pc = pn - pw
678 switch (pressureType_)
679 {
680 case pw:
681 {
682 velocityW *= lambdaW * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermW);
683 velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermNw)
684 + 0.5 * (lambdaNwI + lambdaNwBound) * scalarPerm * (pcI - pcBound) / dist;
685 break;
686 }
687 case pn:
688 {
689 velocityW *= lambdaW * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermW)
690 - 0.5 * (lambdaWI + lambdaWBound) * scalarPerm * (pcI - pcBound) / dist;
691 velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermNw);
692 break;
693 }
694 case pGlobal:
695 {
696 velocityW *= (lambdaW + lambdaNw) * scalarPerm * (cellData.globalPressure() - pressBound) / dist
697 + scalarPerm * (lambdaW * gravityTermW + lambdaNw * gravityTermNw);
698 velocityNw = 0;
699 break;
700 }
701 }
702
703 // store velocities
704 cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
705 cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
706 cellData.fluxData().setVelocityMarker(isIndex);
707
708 } // end Dirichlet boundary
709
710 else if (bcType.isNeumann(eqIdxPress))
711 {
712 problem_.neumann(boundValues, intersection);
713
714 Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
715 Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
716
717 velocityW *= boundValues[wPhaseIdx];
718 velocityNw *= boundValues[nPhaseIdx];
719
720 if (!compressibility_)
721 {
722 velocityW /= density_[wPhaseIdx];
723 velocityNw /= density_[nPhaseIdx];
724 }
725
726 // store potential gradients for further calculations
727 cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, boundValues[wPhaseIdx]);
728 cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, boundValues[nPhaseIdx]);
729
730 cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
731 cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
732 cellData.fluxData().setVelocityMarker(isIndex);
733 } // end Neumann boundary
734 else
735 {
736 DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!");
737 }
738
739// printvector(std::cout, cellData.fluxData().velocity(), "velocity", "row", 4, 1, 3);
740 return;
741}
742} // end namespace Dumux
743#endif
#define GET_PROP_VALUE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:282
#define GET_PROP(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:281
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag EnableCompressibility
Returns whether compressibility is allowed.
Definition: porousmediumflow/2p/sequential/properties.hh:55
Property tag SaturationFormulation
The formulation of the saturation model.
Definition: porousmediumflow/2p/sequential/properties.hh:53
Property tag NumPhases
Number of phases in the system.
Definition: porousmediumflow/sequential/properties.hh:69
Property tag VelocityFormulation
The type of velocity reconstructed for the transport model.
Definition: porousmediumflow/2p/sequential/properties.hh:54
std::string temperature() noexcept
I/O name of temperature for equilibrium models.
Definition: name.hh:51
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Determines the velocity from a finite volume solution of the pressure equation of a sequential model ...
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:61
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:186
void calculateVelocity(const Intersection &, CellData &)
Calculates the velocity at a cell-cell interface.
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:334
bool calculateVelocityInTransport()
Indicates if velocity is reconstructed in the pressure step or in the transport step.
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:172
FVVelocity2P(Problem &problem)
Constructs a FVVelocity2P object.
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:121
void initialize()
For initialization.
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:143
void calculateVelocityOnBoundary(const Intersection &, CellData &)
Calculates the velocity at a boundary.
Definition: 2p/sequential/diffusion/cellcentered/velocity.hh:494
Specifies the properties for immiscible 2p diffusion/pressure models.