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