3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
lmethod/2dvelocity.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_FVMPFAL2DVELOCITY2P_HH
25#define DUMUX_FVMPFAL2DVELOCITY2P_HH
26
27#include <dune/grid/common/gridenums.hh>
28#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
29#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
30#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh>
31#include "2dtransmissibilitycalculator.hh"
32
33namespace Dumux {
34
57template<class TypeTag> class FvMpfaL2dVelocity2p
58{
59 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
60 enum
61 {
62 dim = GridView::dimension, dimWorld = GridView::dimensionworld
63 };
64
65 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
66 using Problem = GetPropType<TypeTag, Properties::Problem>;
67
68 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
69
70 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
71
72 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
73 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
74
75 using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
76 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
77 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
78 using CellData = GetPropType<TypeTag, Properties::CellData>;
79
80 using Element = typename GridView::Traits::template Codim<0>::Entity;
81 using Grid = typename GridView::Grid;
82 using IndexSet = typename GridView::IndexSet;
83
84 using Geometry = typename Element::Geometry;
85 using JacobianTransposed = typename Geometry::JacobianTransposed ;
86
87 using GridTypeIndices = GetPropType<TypeTag, Properties::GridTypeIndices>;
88
89 using InteractionVolume = FVMPFALInteractionVolume<TypeTag>;
90 using TransmissibilityCalculator = FvMpfaL2dTransmissibilityCalculator<TypeTag>;
91 using InnerBoundaryVolumeFaces = std::vector<Dune::FieldVector<bool, 2*dim> >;
92
93 enum
94 {
95 pw = Indices::pressureW,
96 pn = Indices::pressureNw,
97 pGlobal = Indices::pressureGlobal,
98 sw = Indices::saturationW,
99 sn = Indices::saturationNw,
100 vw = Indices::velocityW,
101 vn = Indices::velocityNw,
102 vt = Indices::velocityTotal
103 };
104 enum
105 {
106 wPhaseIdx = Indices::wPhaseIdx,
107 nPhaseIdx = Indices::nPhaseIdx,
108 pressureIdx = Indices::pressureIdx,
109 saturationIdx = Indices::saturationIdx,
110 pressureEqIdx = Indices::pressureEqIdx,
111 satEqIdx = Indices::satEqIdx,
112 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
113 };
114
115 using LocalPosition = Dune::FieldVector<Scalar, dim>;
116 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
117 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
118 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
119 using DimVector = Dune::FieldVector<Scalar, dim>;
120
121public:
127 FvMpfaL2dVelocity2p(Problem& problem) :
128 problem_(problem), transmissibilityCalculator_(problem), gravity_(problem.gravity())
129 {
130 density_[wPhaseIdx] = 0.;
131 density_[nPhaseIdx] = 0.;
132 viscosity_[wPhaseIdx] = 0.;
133 viscosity_[nPhaseIdx] = 0.;
134
135 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
136 }
137
139 void calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
140 CellData& cellData1, CellData& cellData2, CellData& cellData3, CellData& cellData4,
141 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces);
142 void calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume,
143 CellData& cellData, int elemIdx);
144
146 void initialize()
147 {
148 const auto element = *problem_.gridView().template begin<0>();
149 FluidState fluidState;
150 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
151 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
152 fluidState.setTemperature(problem_.temperature(element));
153 fluidState.setSaturation(wPhaseIdx, 1.);
154 fluidState.setSaturation(nPhaseIdx, 0.);
155 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
156 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
157 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
158 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
159
160 return;
161 }
162
173 template<class MultiWriter>
174 void addOutputVtkFields(MultiWriter &writer)
175 {
176 if (vtkOutputLevel_ > 0)
177 {
178 Dune::BlockVector < DimVector > &velocityWetting
179 = *(writer.template allocateManagedBuffer<Scalar,dim>(problem_.gridView().size(0)));
180 Dune::BlockVector < DimVector > &velocityNonwetting
181 = *(writer.template allocateManagedBuffer<Scalar,dim>(problem_.gridView().size(0)));
182
183 // compute update vector
184 for (const auto& element : elements(problem_.gridView()))
185 {
186 // cell index
187 int eIdxGlobal = problem_.variables().index(element);
188
189 CellData & cellData = problem_.variables().cellData(eIdxGlobal);
190
191 Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
192 Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
193
194 // run through all intersections with neighbors and boundary
195 for (const auto& intersection : intersections(problem_.gridView(), element))
196 {
197 int isIndex = intersection.indexInInside();
198
199 fluxW[isIndex] += intersection.geometry().volume()
200 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
201 fluxNw[isIndex] += intersection.geometry().volume()
202 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
203 }
204
205 DimVector refVelocity(0);
206 refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
207 refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
208
209 const DimVector& localPos = referenceElement(element).position(0, 0);
210
211 // get the transposed Jacobian of the element mapping
212 const JacobianTransposed jacobianT = element.geometry().jacobianTransposed(localPos);
213
214 // calculate the element velocity by the Piola transformation
215 DimVector elementVelocity(0);
216 jacobianT.umtv(refVelocity, elementVelocity);
217 elementVelocity /= element.geometry().integrationElement(localPos);
218
219 velocityWetting[eIdxGlobal] = elementVelocity;
220
221 refVelocity = 0;
222 refVelocity[0] = 0.5 * (fluxNw[1] - fluxNw[0]);
223 refVelocity[1] = 0.5 * (fluxNw[3] - fluxNw[2]);
224
225 // calculate the element velocity by the Piola transformation
226 elementVelocity = 0;
227 jacobianT.umtv(refVelocity, elementVelocity);
228 elementVelocity /= element.geometry().integrationElement(localPos);
229
230 velocityNonwetting[eIdxGlobal] = elementVelocity;
231 }
232
233 writer.attachCellData(velocityWetting, "wetting-velocity", dim);
234 writer.attachCellData(velocityNonwetting, "nonwetting-velocity", dim);
235 }
236
237 return;
238 }
239
240private:
241 Problem& problem_;
242protected:
243 TransmissibilityCalculator transmissibilityCalculator_;
244
245 const GravityVector& gravity_; // vector including the gravity constant
246
247 Scalar density_[numPhases];
248 Scalar viscosity_[numPhases];
249
250 int vtkOutputLevel_;
251
252 static constexpr Scalar threshold_ = 1e-15;
254 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
256 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
258 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
259};
260// end of template
261
275template<class TypeTag>
276void FvMpfaL2dVelocity2p<TypeTag>::calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
277 CellData& cellData1, CellData& cellData2,
278 CellData& cellData3, CellData& cellData4,
279 InnerBoundaryVolumeFaces& innerBoundaryVolumeFaces)
280{
281 auto element1 = interactionVolume.getSubVolumeElement(0);
282 auto element2 = interactionVolume.getSubVolumeElement(1);
283 auto element3 = interactionVolume.getSubVolumeElement(2);
284 auto element4 = interactionVolume.getSubVolumeElement(3);
285
286 int level1 = element1.level();
287 int level2 = element2.level();
288 int level3 = element3.level();
289 int level4 = element4.level();
290
291 // cell index
292 int eIdxGlobal1 = problem_.variables().index(element1);
293 int eIdxGlobal2 = problem_.variables().index(element2);
294 int eIdxGlobal3 = problem_.variables().index(element3);
295 int eIdxGlobal4 = problem_.variables().index(element4);
296
297 // get pressure values
298 Dune::FieldVector < Scalar, 2 * dim > potW(0);
299 Dune::FieldVector < Scalar, 2 * dim > potNw(0);
300
301 potW[0] = cellData1.potential(wPhaseIdx);
302 potW[1] = cellData2.potential(wPhaseIdx);
303 potW[2] = cellData3.potential(wPhaseIdx);
304 potW[3] = cellData4.potential(wPhaseIdx);
305
306 potNw[0] = cellData1.potential(nPhaseIdx);
307 potNw[1] = cellData2.potential(nPhaseIdx);
308 potNw[2] = cellData3.potential(nPhaseIdx);
309 potNw[3] = cellData4.potential(nPhaseIdx);
310
311 //get mobilities of the phases
312 Dune::FieldVector < Scalar, numPhases > lambda1(cellData1.mobility(wPhaseIdx));
313 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
314
315 //compute total mobility of cell 1
316 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
317
318 //get mobilities of the phases
319 Dune::FieldVector < Scalar, numPhases > lambda2(cellData2.mobility(wPhaseIdx));
320 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
321
322 //compute total mobility of cell 1
323 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
324
325 //get mobilities of the phases
326 Dune::FieldVector < Scalar, numPhases > lambda3(cellData3.mobility(wPhaseIdx));
327 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
328
329 //compute total mobility of cell 1
330 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
331
332 //get mobilities of the phases
333 Dune::FieldVector < Scalar, numPhases > lambda4(cellData4.mobility(wPhaseIdx));
334 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
335
336 //compute total mobility of cell 1
337 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
338
339
340 std::vector < DimVector > lambda(2 * dim);
341
342 lambda[0][0] = lambdaTotal1;
343 lambda[0][1] = lambdaTotal1;
344 lambda[1][0] = lambdaTotal2;
345 lambda[1][1] = lambdaTotal2;
346 lambda[2][0] = lambdaTotal3;
347 lambda[2][1] = lambdaTotal3;
348 lambda[3][0] = lambdaTotal4;
349 lambda[3][1] = lambdaTotal4;
350
351 Scalar potentialDiffW12 = 0;
352 Scalar potentialDiffW14 = 0;
353 Scalar potentialDiffW32 = 0;
354 Scalar potentialDiffW34 = 0;
355
356 Scalar potentialDiffNw12 = 0;
357 Scalar potentialDiffNw14 = 0;
358 Scalar potentialDiffNw32 = 0;
359 Scalar potentialDiffNw34 = 0;
360
361 //flux vector
362 Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
363 Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
364
365 Dune::FieldMatrix < Scalar, dim, 2 * dim - dim + 1 > T(0);
366 DimVector Tu(0);
367 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
368
369 int lType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 0, 1, 2, 3);
370
371 if (lType == TransmissibilityCalculator::rightTriangle)
372 {
373 u[0] = potW[1];
374 u[1] = potW[2];
375 u[2] = potW[0];
376
377 T.mv(u, Tu);
378
379 fluxW[0] = Tu[1];
380 potentialDiffW12 = Tu[1];
381
382 u[0] = potNw[1];
383 u[1] = potNw[2];
384 u[2] = potNw[0];
385
386 T.mv(u, Tu);
387
388 fluxNw[0] = Tu[1];
389 potentialDiffNw12 = Tu[1];
390 }
391 else
392 {
393 u[0] = potW[0];
394 u[1] = potW[3];
395 u[2] = potW[1];
396
397 T.mv(u, Tu);
398
399 fluxW[0] = Tu[1];
400 potentialDiffW12 = Tu[1];
401
402 u[0] = potNw[0];
403 u[1] = potNw[3];
404 u[2] = potNw[1];
405
406 T.mv(u, Tu);
407
408 fluxNw[0] = Tu[1];
409 potentialDiffNw12 = Tu[1];
410 }
411
412 lType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 1, 2, 3, 0);
413
414 if (lType == TransmissibilityCalculator::rightTriangle)
415 {
416 u[0] = potW[2];
417 u[1] = potW[3];
418 u[2] = potW[1];
419
420 T.mv(u, Tu);
421
422 fluxW[1] = Tu[1];
423 potentialDiffW32 = -Tu[1];
424
425 u[0] = potNw[2];
426 u[1] = potNw[3];
427 u[2] = potNw[1];
428
429 T.mv(u, Tu);
430
431 fluxNw[1] = Tu[1];
432 potentialDiffNw32 = -Tu[1];
433 }
434 else
435 {
436 u[0] = potW[1];
437 u[1] = potW[0];
438 u[2] = potW[2];
439
440 T.mv(u, Tu);
441
442 fluxW[1] = Tu[1];
443 potentialDiffW32 = -Tu[1];
444
445 u[0] = potNw[1];
446 u[1] = potNw[0];
447 u[2] = potNw[2];
448
449 T.mv(u, Tu);
450
451 fluxNw[1] = Tu[1];
452 potentialDiffNw32 = -Tu[1];
453 }
454
455 lType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 2, 3, 0, 1);
456
457 if (lType == TransmissibilityCalculator::rightTriangle)
458 {
459 u[0] = potW[3];
460 u[1] = potW[0];
461 u[2] = potW[2];
462
463 T.mv(u, Tu);
464
465 fluxW[2] = Tu[1];
466 potentialDiffW34 = Tu[1];
467
468 u[0] = potNw[3];
469 u[1] = potNw[0];
470 u[2] = potNw[2];
471
472 T.mv(u, Tu);
473
474 fluxNw[2] = Tu[1];
475 potentialDiffNw34 = Tu[1];
476 }
477 else
478 {
479 u[0] = potW[2];
480 u[1] = potW[1];
481 u[2] = potW[3];
482
483 T.mv(u, Tu);
484
485 fluxW[2] = Tu[1];
486 potentialDiffW34 = Tu[1];
487
488 u[0] = potNw[2];
489 u[1] = potNw[1];
490 u[2] = potNw[3];
491
492 T.mv(u, Tu);
493
494 fluxNw[2] = Tu[1];
495 potentialDiffNw34 = Tu[1];
496 }
497
498 lType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 3, 0, 1, 2);
499
500 if (lType == TransmissibilityCalculator::rightTriangle)
501 {
502 u[0] = potW[0];
503 u[1] = potW[1];
504 u[2] = potW[3];
505
506 T.mv(u, Tu);
507
508 fluxW[3] = Tu[1];
509 potentialDiffW14 = -Tu[1];
510
511 u[0] = potNw[0];
512 u[1] = potNw[1];
513 u[2] = potNw[3];
514
515 T.mv(u, Tu);
516
517 fluxNw[3] = Tu[1];
518 potentialDiffNw14 = -Tu[1];
519 }
520 else
521 {
522 u[0] = potW[3];
523 u[1] = potW[2];
524 u[2] = potW[0];
525
526 T.mv(u, Tu);
527
528 fluxW[3] = Tu[1];
529 potentialDiffW14 = -Tu[1];
530
531 u[0] = potNw[3];
532 u[1] = potNw[2];
533 u[2] = potNw[0];
534
535 T.mv(u, Tu);
536
537 fluxNw[3] = Tu[1];
538 potentialDiffNw14 = -Tu[1];
539 }
540
541 //store potentials for further calculations (saturation, ...)
542 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffW12);
543 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffNw12);
544 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffW14);
545 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffNw14);
546 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 0), -potentialDiffW32);
547 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), -potentialDiffNw32);
548 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffW12);
549 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffNw12);
550 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffW34);
551 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffNw34);
552 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 1), potentialDiffW32);
553 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 1), potentialDiffNw32);
554 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffW14);
555 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffNw14);
556 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffW34);
557 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffNw34);
558
559 //compute mobilities of face 1
560 Dune::FieldVector < Scalar, numPhases > lambda12Upw(0.0);
561 lambda12Upw[wPhaseIdx] = (potentialDiffW12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
562 lambda12Upw[nPhaseIdx] = (potentialDiffNw12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
563
564 //compute mobilities of face 4
565 Dune::FieldVector < Scalar, numPhases > lambda14Upw(0.0);
566 lambda14Upw[wPhaseIdx] = (potentialDiffW14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
567 lambda14Upw[nPhaseIdx] = (potentialDiffNw14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
568
569 //compute mobilities of face 2
570 Dune::FieldVector < Scalar, numPhases > lambda32Upw(0.0);
571 lambda32Upw[wPhaseIdx] = (potentialDiffW32 >= 0) ? lambda3[wPhaseIdx] : lambda2[wPhaseIdx];
572 lambda32Upw[nPhaseIdx] = (potentialDiffNw32 >= 0) ? lambda3[nPhaseIdx] : lambda2[nPhaseIdx];
573
574 //compute mobilities of face 3
575 Dune::FieldVector < Scalar, numPhases > lambda34Upw(0.0);
576 lambda34Upw[wPhaseIdx] = (potentialDiffW34 >= 0) ? lambda3[wPhaseIdx] : lambda4[wPhaseIdx];
577 lambda34Upw[nPhaseIdx] = (potentialDiffNw34 >= 0) ? lambda3[nPhaseIdx] : lambda4[nPhaseIdx];
578
579 for (int i = 0; i < numPhases; i++)
580 {
581 // evaluate parts of velocity
582 DimVector vel12 = interactionVolume.getNormal(0, 0);
583 DimVector vel14 = interactionVolume.getNormal(3, 0);
584 DimVector vel23 = interactionVolume.getNormal(1, 0);
585 DimVector vel21 = interactionVolume.getNormal(0, 0);
586 DimVector vel34 = interactionVolume.getNormal(2, 0);
587 DimVector vel32 = interactionVolume.getNormal(1, 0);
588 DimVector vel41 = interactionVolume.getNormal(3, 0);
589 DimVector vel43 = interactionVolume.getNormal(2, 0);
590
591 Dune::FieldVector < Scalar, 2 * dim > flux(0);
592 switch (i)
593 {
594 case wPhaseIdx:
595 {
596 flux = fluxW;
597 break;
598 }
599 case nPhaseIdx:
600 {
601 flux = fluxNw;
602 break;
603 }
604 }
605
606 vel12 *= flux[0] / (2 * interactionVolume.getFaceArea(0, 0)); //divide by 2 because the flux is related to the half face!
607 vel14 *= flux[3] / (2 * interactionVolume.getFaceArea(0, 1));
608 vel23 *= flux[1] / (2 * interactionVolume.getFaceArea(1, 0));
609 vel21 *= flux[0] / (2 * interactionVolume.getFaceArea(1, 1));
610 vel34 *= flux[2] / (2 * interactionVolume.getFaceArea(2, 0));
611 vel32 *= flux[1] / (2 * interactionVolume.getFaceArea(2, 1));
612 vel41 *= flux[3] / (2 * interactionVolume.getFaceArea(3, 0));
613 vel43 *= flux[2] / (2 * interactionVolume.getFaceArea(3, 1));
614
615 if (level1 < level2)
616 {
617 vel12 *= 0.5;
618 }
619 else if (level2 < level1)
620 {
621 vel21 *= 0.5;
622 }
623 if (level2 < level3)
624 {
625 vel23 *= 0.5;
626 }
627 else if (level3 < level2)
628 {
629 vel32 *= 0.5;
630 }
631 if (level3 < level4)
632 {
633 vel34 *= 0.5;
634 }
635 else if (level4 < level3)
636 {
637 vel43 *= 0.5;
638 }
639 if (level4 < level1)
640 {
641 vel41 *= 0.5;
642 }
643 else if (level1 < level4)
644 {
645 vel14 *= 0.5;
646 }
647
648 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
649 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
650 Scalar lambdaT32 = lambda32Upw[wPhaseIdx] + lambda32Upw[nPhaseIdx];
651 Scalar lambdaT34 = lambda34Upw[wPhaseIdx] + lambda34Upw[nPhaseIdx];
652 Scalar fracFlow12 = (lambdaT12 > threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
653 Scalar fracFlow14 = (lambdaT14 > threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
654 Scalar fracFlow32 = (lambdaT32 > threshold_) ? lambda32Upw[i] / (lambdaT32) : 0.0;
655 Scalar fracFlow34 = (lambdaT34 > threshold_) ? lambda34Upw[i] / (lambdaT34) : 0.0;
656
657 vel12 *= fracFlow12;
658 vel14 *= fracFlow14;
659 vel23 *= fracFlow32;
660 vel21 *= fracFlow12;
661 vel34 *= fracFlow34;
662 vel32 *= fracFlow32;
663 vel41 *= fracFlow14;
664 vel43 *= fracFlow34;
665
666 if (innerBoundaryVolumeFaces[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
667 {
668 vel12 *= 2;
669 }
670 if (innerBoundaryVolumeFaces[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
671 {
672 vel14 *= 2;
673 }
674 if (innerBoundaryVolumeFaces[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
675 {
676 vel23 *= 2;
677 }
678 if (innerBoundaryVolumeFaces[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
679 {
680 vel21 *= 2;
681 }
682 if (innerBoundaryVolumeFaces[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
683 {
684 vel34 *= 2;
685 }
686 if (innerBoundaryVolumeFaces[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
687 {
688 vel32 *= 2;
689 }
690 if (innerBoundaryVolumeFaces[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
691 {
692 vel41 *= 2;
693 }
694 if (innerBoundaryVolumeFaces[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
695 {
696 vel43 *= 2;
697 }
698
699 //store velocities
700 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 0), vel12);
701 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 1), vel14);
702 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 0), vel23);
703 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 1), vel21);
704 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 0), vel34);
705 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 1), vel32);
706 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 0), vel41);
707 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 1), vel43);
708 }
709 //set velocity marker
710 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 0));
711 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 1));
712 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 0));
713 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 1));
714 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 0));
715 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 1));
716 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 0));
717 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
718}
719
729template<class TypeTag>
730void FvMpfaL2dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume,
731 CellData& cellData, int elemIdx)
732{
733 auto element = interactionVolume.getSubVolumeElement(elemIdx);
734
735 // get global coordinate of cell centers
736 const GlobalPosition& globalPos = element.geometry().center();
737
738 //permeability vector at boundary
739 DimMatrix permeability(problem_.spatialParams().intrinsicPermeability(element));
740
741 //get mobilities of the phases
742 Dune::FieldVector < Scalar, numPhases > lambda(cellData.mobility(wPhaseIdx));
743 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
744
745 for (int fIdx = 0; fIdx < dim; fIdx++)
746 {
747 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
748
749 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
750 {
751 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressureEqIdx))
752 {
753 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
754
755 const auto refElement = referenceElement(element);
756
757 const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
758
759 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
760
761 DimVector distVec(globalPosFace - globalPos);
762 Scalar dist = distVec.two_norm();
763 DimVector unitDistVec(distVec);
764 unitDistVec /= dist;
765
766 // get pc and lambda at the boundary
767 Scalar satWBound = cellData.saturation(wPhaseIdx);
768 //check boundary sat at face 1
769 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
770 {
771 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
772 switch (saturationType_)
773 {
774 case sw:
775 {
776 satWBound = satBound;
777 break;
778 }
779 case sn:
780 {
781 satWBound = 1 - satBound;
782 break;
783 }
784 }
785
786 }
787
788 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
789 Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
790
791 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
792 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
793
794 pcBound += gravityDiffBound;
795
796 Dune::FieldVector <Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
797 lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
798 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
799 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
800
801 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
802 Scalar potentialBoundW = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx] + density_[wPhaseIdx]*gdeltaZ;
803 Scalar potentialBoundNw = potentialBoundW;
804
805 //calculate potential gradients
806 switch (pressureType_)
807 {
808 case pw:
809 {
810 potentialBoundNw += pcBound;
811 break;
812 }
813 case pn:
814 {
815 //calculate potential gradients
816 potentialBoundW -= pcBound;
817 break;
818 }
819 }
820
821 Scalar potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
822 Scalar potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
823
824 //store potentials for further calculations (saturation, ...)
825 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, potentialDiffW);
826 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, potentialDiffNw);
827
828 //calculated phase velocities from advective velocities -> capillary pressure velocity already added in pressure part!
829 DimVector velocityW(0);
830 DimVector velocityNw(0);
831
832 // calculate capillary pressure gradient
833 DimVector pressGradient = unitDistVec;
834 pressGradient *= (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
835 permeability.mv(pressGradient, velocityW);
836
837 pressGradient = unitDistVec;
838 pressGradient *= (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
839 permeability.mv(pressGradient, velocityNw);
840
841 velocityW *= (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
842 velocityNw *= (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
843
844 //velocity is calculated from two vertices of one intersection!
845 velocityW *= 0.5;
846 velocityNw *= 0.5;
847
848 //store velocities
849 velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
850 velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
851 cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
852 cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
853 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
854 }
855 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx))
856 {
857 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
858
859 const auto refElement = referenceElement(element);
860
861 const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
862
863 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
864
865 DimVector distVec(globalPosFace - globalPos);
866 Scalar dist = distVec.two_norm();
867 DimVector unitDistVec(distVec);
868 unitDistVec /= dist;
869
870 // get neumann boundary value
871 PrimaryVariables boundValues(interactionVolume.getNeumannValues(intVolFaceIdx));
872
873 boundValues[wPhaseIdx] /= density_[wPhaseIdx];
874 boundValues[nPhaseIdx] /= density_[nPhaseIdx];
875
876 DimVector velocityW(unitDistVec);
877 DimVector velocityNw(unitDistVec);
878
879 velocityW *= boundValues[wPhaseIdx] / (2 * interactionVolume.getFaceArea(elemIdx, fIdx));
880 velocityNw *= boundValues[nPhaseIdx]
881 / (2 * interactionVolume.getFaceArea(elemIdx, fIdx));
882
883 //store potentials for further calculations (saturation, ...)
884 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, boundValues[wPhaseIdx]);
885 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, boundValues[nPhaseIdx]);
886
887 //store velocities
888 velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
889 velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
890 cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
891 cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
892 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
893 }
894 else
895 {
896 DUNE_THROW(Dune::NotImplemented,
897 "No valid boundary condition type defined for pressure equation!");
898 }
899 }
900 }
901}
902
903} // end namespace Dumux
904#endif
2dtransmissibilitycalculator.hh
Provides methods for transmissibility calculation 2-d.
linteractionvolume.hh
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Dumux
Definition: adapt.hh:29
Dumux::GetProp
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property
Definition: propertysystem.hh:141
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::BoundaryTypes::isNeumann
bool isNeumann(unsigned eqIdx) const
Returns true if an equation is used to specify a Neumann condition.
Definition: common/boundarytypes.hh:273
Dumux::BoundaryTypes::isDirichlet
bool isDirichlet(unsigned eqIdx) const
Returns true if an equation is used to specify a Dirichlet condition.
Definition: common/boundarytypes.hh:236
Dumux::FvMpfaL2dTransmissibilityCalculator
Provides methods for transmissibility calculation in 2-d.
Definition: 2dtransmissibilitycalculator.hh:44
Dumux::FvMpfaL2dVelocity2p
Class for calculating 2-d velocities from cell-wise constant pressure values.
Definition: lmethod/2dvelocity.hh:58
Dumux::FvMpfaL2dVelocity2p::viscosity_
Scalar viscosity_[numPhases]
Definition: lmethod/2dvelocity.hh:248
Dumux::FvMpfaL2dVelocity2p::addOutputVtkFields
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: lmethod/2dvelocity.hh:174
Dumux::FvMpfaL2dVelocity2p::gravity_
const GravityVector & gravity_
Definition: lmethod/2dvelocity.hh:245
Dumux::FvMpfaL2dVelocity2p::vtkOutputLevel_
int vtkOutputLevel_
Definition: lmethod/2dvelocity.hh:250
Dumux::FvMpfaL2dVelocity2p::calculateBoundaryInteractionVolumeVelocity
void calculateBoundaryInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData, int elemIdx)
Calculates the velocity at a boundary flux faces.
Definition: lmethod/2dvelocity.hh:730
Dumux::FvMpfaL2dVelocity2p::initialize
void initialize()
Initializes the velocity model.
Definition: lmethod/2dvelocity.hh:146
Dumux::FvMpfaL2dVelocity2p::velocityType_
static const int velocityType_
gives kind of velocity used ( , , )
Definition: lmethod/2dvelocity.hh:254
Dumux::FvMpfaL2dVelocity2p::saturationType_
static const int saturationType_
gives kind of saturation used ( , )
Definition: lmethod/2dvelocity.hh:258
Dumux::FvMpfaL2dVelocity2p::transmissibilityCalculator_
TransmissibilityCalculator transmissibilityCalculator_
Definition: lmethod/2dvelocity.hh:243
Dumux::FvMpfaL2dVelocity2p::calculateInnerInteractionVolumeVelocity
void calculateInnerInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData1, CellData &cellData2, CellData &cellData3, CellData &cellData4, InnerBoundaryVolumeFaces &innerBoundaryVolumeFaces)
Calculate velocities for flux faces of an interaction volume.
Definition: lmethod/2dvelocity.hh:276
Dumux::FvMpfaL2dVelocity2p::pressureType_
static const int pressureType_
gives kind of pressure used ( , , )
Definition: lmethod/2dvelocity.hh:256
Dumux::FvMpfaL2dVelocity2p::density_
Scalar density_[numPhases]
Definition: lmethod/2dvelocity.hh:247
Dumux::FvMpfaL2dVelocity2p::threshold_
static constexpr Scalar threshold_
Definition: lmethod/2dvelocity.hh:252
Dumux::FvMpfaL2dVelocity2p::FvMpfaL2dVelocity2p
FvMpfaL2dVelocity2p(Problem &problem)
Constructs a FvMpfaL2dVelocity2p object.
Definition: lmethod/2dvelocity.hh:127
Dumux::FVMPFALInteractionVolume
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Definition: linteractionvolume.hh:41
Dumux::FVMPFALInteractionVolume::getIndexOnElement
int getIndexOnElement(int subVolumeIdx, int subVolumeFaceIdx)
Map from local interaction volume numbering to numbering of the Dune reference element.
Definition: linteractionvolume.hh:258
Dumux::FVMPFALInteractionVolume::getNormal
DimVector & getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face normal.
Definition: linteractionvolume.hh:364
Dumux::FVMPFALInteractionVolume::getSubVolumeElement
Element getSubVolumeElement(int subVolumeIdx)
Get an element of the interaction volume.
Definition: linteractionvolume.hh:281
Dumux::FVMPFALInteractionVolume::getBoundaryType
BoundaryTypes & getBoundaryType(int subVolumeFaceIdx)
Get boundary condtion types for a flux face.
Definition: linteractionvolume.hh:292
Dumux::FVMPFALInteractionVolume::getFaceArea
Scalar & getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face area.
Definition: linteractionvolume.hh:388
Dumux::FVMPFALInteractionVolume::getNeumannValues
PrimaryVariables & getNeumannValues(int subVolumeFaceIdx)
Get the Neumann boundary condtions for a flux face.
Definition: linteractionvolume.hh:352
Dumux::FVMPFALInteractionVolume::isBoundaryFace
bool isBoundaryFace(int subVolumeFaceIdx)
Returns true if an interaction volume flux face is a boundary face.
Definition: linteractionvolume.hh:330
Dumux::FVMPFALInteractionVolume::getFaceIndexFromSubVolume
int getFaceIndexFromSubVolume(int subVolumeIdx, int subVolumeFaceIdx)
Map from local interaction volume numbering on element to numbering on interaction volume.
Definition: linteractionvolume.hh:270
Dumux::FVMPFALInteractionVolume::getDirichletValues
PrimaryVariables & getDirichletValues(int subVolumeFaceIdx)
Get the Dirichlet boundary condtions for a flux face.
Definition: linteractionvolume.hh:341
properties.hh
Specifies the properties for immiscible 2p diffusion/pressure models.
properties.hh
Properties for a MPFA method.
Generated by doxygen 1.9.3