3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
3dvelocity.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_FVMPFAL2PVELOCITY2P_HH
25#define DUMUX_FVMPFAL2PVELOCITY2P_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 "3dtransmissibilitycalculator.hh"
31
32namespace Dumux {
33
54template<class TypeTag> class FvMpfaL3dVelocity2p
55{
56 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
57 enum
58 {
59 dim = GridView::dimension, dimWorld = GridView::dimensionworld
60 };
61
62 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
63 using Problem = GetPropType<TypeTag, Properties::Problem>;
64
65 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
66
67 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
68
69 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
70 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
71
72 using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
73 using CellData = GetPropType<TypeTag, Properties::CellData>;
74
75 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
76 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
77
78 using Element = typename GridView::Traits::template Codim<0>::Entity;
79 using Grid = typename GridView::Grid;
80 using IndexSet = typename GridView::IndexSet;
81
82 using Geometry = typename Element::Geometry;
83 using JacobianTransposed = typename Geometry::JacobianTransposed ;
84
85 using GridTypeIndices = GetPropType<TypeTag, Properties::GridTypeIndices>;
86
87 using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
88 using InteractionVolumeContainer = GetPropType<TypeTag, Properties::MPFAInteractionVolumeContainer>;
89 using TransmissibilityCalculator = FvMpfaL3dTransmissibilityCalculator<TypeTag>;
90 using TransmissibilityType = typename TransmissibilityCalculator::TransmissibilityType;
91
92 enum
93 {
94 pw = Indices::pressureW,
95 pn = Indices::pressureNw,
96 pGlobal = Indices::pressureGlobal,
97 sw = Indices::saturationW,
98 sn = Indices::saturationNw,
99 vw = Indices::velocityW,
100 vn = Indices::velocityNw,
101 vt = Indices::velocityTotal
102 };
103 enum
104 {
105 wPhaseIdx = Indices::wPhaseIdx,
106 nPhaseIdx = Indices::nPhaseIdx,
107 pressureIdx = Indices::pressureIdx,
108 saturationIdx = Indices::saturationIdx,
109 pressureEqIdx = Indices::pressureEqIdx,
110 satEqIdx = Indices::satEqIdx,
111 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
112 };
113 enum
114 {
115 globalCorner = 2,
116 globalEdge = 3,
117 neumannNeumann = 0,
118 dirichletDirichlet = 1,
119 dirichletNeumann = 2,
120 neumannDirichlet = 3
121 };
122 enum
123 {
124 innerEdgeFace = 2, innerSideFace = 1
125 };
126
127 using GlobalPosition = typename Geometry::GlobalCoordinate;
128 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
129 using DimVector = Dune::FieldVector<Scalar, dim>;
130
131public:
137 FvMpfaL3dVelocity2p(Problem& problem) :
138 problem_(problem), gravity_(problem.gravity())
139 {
140 density_[wPhaseIdx] = 0.;
141 density_[nPhaseIdx] = 0.;
142 viscosity_[wPhaseIdx] = 0.;
143 viscosity_[nPhaseIdx] = 0.;
144
145 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
146 }
147
149 void calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
150 CellData & cellData1, CellData & cellData2, CellData & cellData3, CellData & cellData4,
151 CellData & cellData5, CellData & cellData6, CellData & cellData7, CellData & cellData8,
152 InteractionVolumeContainer& interactionVolumes,
153 TransmissibilityCalculator& transmissibilityCalculator, int fIdx = -1);
154 void calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume,
155 CellData& cellData, int elemIdx);
156
158 void initialize(bool solveTwice = true)
159 {
160 const auto element = *problem_.gridView().template begin<0>();
161 FluidState fluidState;
162 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
163 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
164 fluidState.setTemperature(problem_.temperature(element));
165 fluidState.setSaturation(wPhaseIdx, 1.);
166 fluidState.setSaturation(nPhaseIdx, 0.);
167 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
168 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
169 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
170 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
171
172 return;
173 }
174
185 template<class MultiWriter>
186 void addOutputVtkFields(MultiWriter &writer)
187 {
188 if (vtkOutputLevel_ > 0)
189 {
190 Dune::BlockVector < DimVector > &velocityWetting = *(writer.template allocateManagedBuffer<Scalar,
191 dim>(problem_.gridView().size(0)));
192 Dune::BlockVector < DimVector > &velocityNonwetting = *(writer.template allocateManagedBuffer<Scalar,
193 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 Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
204 Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
205
206 // run through all intersections with neighbors and boundary
207 for (const auto& intersection : intersections(problem_.gridView(), element))
208 {
209 int isIndex = intersection.indexInInside();
210
211 fluxW[isIndex] += intersection.geometry().volume()
212 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
213 fluxNw[isIndex] += intersection.geometry().volume()
214 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
215 }
216
217 DimVector refVelocity(0);
218 refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
219 refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
220 refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
221
222 const DimVector& localPos = referenceElement(element).position(0, 0);
223
224 // get the transposed Jacobian of the element mapping
225 const JacobianTransposed jacobianT = element.geometry().jacobianTransposed(localPos);
226
227 // calculate the element velocity by the Piola transformation
228 DimVector elementVelocity(0);
229 jacobianT.umtv(refVelocity, elementVelocity);
230 elementVelocity /= element.geometry().integrationElement(localPos);
231
232 velocityWetting[eIdxGlobal] = elementVelocity;
233
234 refVelocity = 0;
235 refVelocity[0] = 0.5 * (fluxNw[1] - fluxNw[0]);
236 refVelocity[1] = 0.5 * (fluxNw[3] - fluxNw[2]);
237 refVelocity[2] = 0.5 * (fluxNw[5] - fluxNw[4]);
238
239 // calculate the element velocity by the Piola transformation
240 elementVelocity = 0;
241 jacobianT.umtv(refVelocity, elementVelocity);
242 elementVelocity /= element.geometry().integrationElement(localPos);
243
244 velocityNonwetting[eIdxGlobal] = elementVelocity;
245 }
246
247 writer.attachCellData(velocityWetting, "wetting-velocity", dim);
248 writer.attachCellData(velocityNonwetting, "nonwetting-velocity", dim);
249 }
250
251 return;
252 }
253
254private:
255 Problem& problem_;
256 const Dune::FieldVector<Scalar, dimWorld>& gravity_;
257
258 Scalar density_[numPhases];
259 Scalar viscosity_[numPhases];
260
261 int vtkOutputLevel_;
262
263 static constexpr Scalar threshold_ = 1e-15;
265 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
267 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
269 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
270};
271// end of template
272
293template<class TypeTag>
294void FvMpfaL3dVelocity2p<TypeTag>::calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
295 CellData & cellData1, CellData & cellData2, CellData & cellData3, CellData & cellData4,
296 CellData & cellData5, CellData & cellData6, CellData & cellData7, CellData & cellData8,
297 InteractionVolumeContainer& interactionVolumes, TransmissibilityCalculator& transmissibilityCalculator, int fIdx)
298 {
299 auto element1 = interactionVolume.getSubVolumeElement(0);
300 auto element2 = interactionVolume.getSubVolumeElement(1);
301 auto element3 = interactionVolume.getSubVolumeElement(2);
302 auto element4 = interactionVolume.getSubVolumeElement(3);
303 auto element5 = interactionVolume.getSubVolumeElement(4);
304 auto element6 = interactionVolume.getSubVolumeElement(5);
305 auto element7 = interactionVolume.getSubVolumeElement(6);
306 auto element8 = interactionVolume.getSubVolumeElement(7);
307
308 // cell index
309 int eIdxGlobal1 = problem_.variables().index(element1);
310 int eIdxGlobal2 = problem_.variables().index(element2);
311 int eIdxGlobal3 = problem_.variables().index(element3);
312 int eIdxGlobal4 = problem_.variables().index(element4);
313 int eIdxGlobal5 = problem_.variables().index(element5);
314 int eIdxGlobal6 = problem_.variables().index(element6);
315 int eIdxGlobal7 = problem_.variables().index(element7);
316 int eIdxGlobal8 = problem_.variables().index(element8);
317
318 // pressures flux calculation
319 Dune::FieldVector<Scalar, 8> potW(0);
320 Dune::FieldVector<Scalar, 8> potNw(0);
321
322 potW[0] = cellData1.potential(wPhaseIdx);
323 potW[1] = cellData2.potential(wPhaseIdx);
324 potW[2] = cellData3.potential(wPhaseIdx);
325 potW[3] = cellData4.potential(wPhaseIdx);
326 potW[4] = cellData5.potential(wPhaseIdx);
327 potW[5] = cellData6.potential(wPhaseIdx);
328 potW[6] = cellData7.potential(wPhaseIdx);
329 potW[7] = cellData8.potential(wPhaseIdx);
330
331 potNw[0] = cellData1.potential(nPhaseIdx);
332 potNw[1] = cellData2.potential(nPhaseIdx);
333 potNw[2] = cellData3.potential(nPhaseIdx);
334 potNw[3] = cellData4.potential(nPhaseIdx);
335 potNw[4] = cellData5.potential(nPhaseIdx);
336 potNw[5] = cellData6.potential(nPhaseIdx);
337 potNw[6] = cellData7.potential(nPhaseIdx);
338 potNw[7] = cellData8.potential(nPhaseIdx);
339
340 // get mobilities of the phases
341 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
342 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
343
344 // compute total mobility of cell 1
345 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
346
347 // get mobilities of the phases
348 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
349 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
350
351 // compute total mobility of cell 2
352 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
353
354 // get mobilities of the phases
355 Dune::FieldVector<Scalar, numPhases> lambda3(cellData3.mobility(wPhaseIdx));
356 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
357
358 // compute total mobility of cell 3
359 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
360
361 // get mobilities of the phases
362 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
363 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
364
365 // compute total mobility of cell 4
366 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
367
368 // get mobilities of the phases
369 Dune::FieldVector<Scalar, numPhases> lambda5(cellData5.mobility(wPhaseIdx));
370 lambda5[nPhaseIdx] = cellData5.mobility(nPhaseIdx);
371
372 // compute total mobility of cell 5
373 Scalar lambdaTotal5 = lambda5[wPhaseIdx] + lambda5[nPhaseIdx];
374
375 // get mobilities of the phases
376 Dune::FieldVector<Scalar, numPhases> lambda6(cellData6.mobility(wPhaseIdx));
377 lambda6[nPhaseIdx] = cellData6.mobility(nPhaseIdx);
378
379 // compute total mobility of cell 6
380 Scalar lambdaTotal6 = lambda6[wPhaseIdx] + lambda6[nPhaseIdx];
381
382 // get mobilities of the phases
383 Dune::FieldVector<Scalar, numPhases> lambda7(cellData7.mobility(wPhaseIdx));
384 lambda7[nPhaseIdx] = cellData7.mobility(nPhaseIdx);
385
386 // compute total mobility of cell 7
387 Scalar lambdaTotal7 = lambda7[wPhaseIdx] + lambda7[nPhaseIdx];
388
389 // get mobilities of the phases
390 Dune::FieldVector<Scalar, numPhases> lambda8(cellData8.mobility(wPhaseIdx));
391 lambda8[nPhaseIdx] = cellData8.mobility(nPhaseIdx);
392
393 // compute total mobility of cell 8
394 Scalar lambdaTotal8 = lambda8[wPhaseIdx] + lambda8[nPhaseIdx];
395
396 std::vector<DimVector> lambda(8);
397 lambda[0][0] = lambdaTotal1;
398 lambda[0][1] = lambdaTotal1;
399 lambda[0][2] = lambdaTotal1;
400 lambda[1][0] = lambdaTotal2;
401 lambda[1][1] = lambdaTotal2;
402 lambda[1][2] = lambdaTotal2;
403 lambda[2][0] = lambdaTotal3;
404 lambda[2][1] = lambdaTotal3;
405 lambda[2][2] = lambdaTotal3;
406 lambda[3][0] = lambdaTotal4;
407 lambda[3][1] = lambdaTotal4;
408 lambda[3][2] = lambdaTotal4;
409 lambda[4][0] = lambdaTotal5;
410 lambda[4][1] = lambdaTotal5;
411 lambda[4][2] = lambdaTotal5;
412 lambda[5][0] = lambdaTotal6;
413 lambda[5][1] = lambdaTotal6;
414 lambda[5][2] = lambdaTotal6;
415 lambda[6][0] = lambdaTotal7;
416 lambda[6][1] = lambdaTotal7;
417 lambda[6][2] = lambdaTotal7;
418 lambda[7][0] = lambdaTotal8;
419 lambda[7][1] = lambdaTotal8;
420 lambda[7][2] = lambdaTotal8;
421
422 Scalar potentialDiffW0 = 0;
423 Scalar potentialDiffW1 = 0;
424 Scalar potentialDiffW2 = 0;
425 Scalar potentialDiffW3 = 0;
426 Scalar potentialDiffW4 = 0;
427 Scalar potentialDiffW5 = 0;
428 Scalar potentialDiffW6 = 0;
429 Scalar potentialDiffW7 = 0;
430 Scalar potentialDiffW8 = 0;
431 Scalar potentialDiffW9 = 0;
432 Scalar potentialDiffW10 = 0;
433 Scalar potentialDiffW11 = 0;
434
435 Scalar potentialDiffNw0 = 0;
436 Scalar potentialDiffNw1 = 0;
437 Scalar potentialDiffNw2 = 0;
438 Scalar potentialDiffNw3 = 0;
439 Scalar potentialDiffNw4 = 0;
440 Scalar potentialDiffNw5 = 0;
441 Scalar potentialDiffNw6 = 0;
442 Scalar potentialDiffNw7 = 0;
443 Scalar potentialDiffNw8 = 0;
444 Scalar potentialDiffNw9 = 0;
445 Scalar potentialDiffNw10 = 0;
446 Scalar potentialDiffNw11 = 0;
447
448 //flux vector
449 Dune::FieldVector<Scalar, 12> fluxW(0);
450 Dune::FieldVector<Scalar, 12> fluxNw(0);
451
452 DimVector Tu(0);
453 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
454 TransmissibilityType T(0);
455
456 if (fIdx < 0 || fIdx == 0)
457 {
458 // calculate the flux through the subvolumeface 1 (subVolumeFaceIdx = 0)
459 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 0, 1, 2, 3,
460 4, 5);
461
462 if (caseL == 1)
463 {
464 u[0] = potW[0];
465 u[1] = potW[1];
466 u[2] = potW[2];
467 u[3] = potW[4];
468
469 T.mv(u, Tu);
470
471 fluxW[0] = Tu[0];
472 potentialDiffW0 = Tu[0];
473
474 u[0] = potNw[0];
475 u[1] = potNw[1];
476 u[2] = potNw[2];
477 u[3] = potNw[4];
478
479 T.mv(u, Tu);
480
481 fluxNw[0] = Tu[0];
482 potentialDiffNw0 = Tu[0];
483 }
484 else if (caseL == 2)
485 {
486 u[0] = potW[0];
487 u[1] = potW[1];
488 u[2] = potW[3];
489 u[3] = potW[5];
490
491 T.mv(u, Tu);
492
493 fluxW[0] = Tu[0];
494 potentialDiffW0 = Tu[0];
495
496 u[0] = potNw[0];
497 u[1] = potNw[1];
498 u[2] = potNw[3];
499 u[3] = potNw[5];
500
501 T.mv(u, Tu);
502
503 fluxNw[0] = Tu[0];
504 potentialDiffNw0 = Tu[0];
505 }
506 else if (caseL == 3)
507 {
508 u[0] = potW[0];
509 u[1] = potW[1];
510 u[2] = potW[3];
511 u[3] = potW[4];
512
513 T.mv(u, Tu);
514
515 fluxW[0] = Tu[0];
516 potentialDiffW0 = Tu[0];
517
518 u[0] = potNw[0];
519 u[1] = potNw[1];
520 u[2] = potNw[3];
521 u[3] = potNw[4];
522
523 T.mv(u, Tu);
524
525 fluxNw[0] = Tu[0];
526 potentialDiffNw0 = Tu[0];
527 }
528 else
529 {
530 u[0] = potW[0];
531 u[1] = potW[1];
532 u[2] = potW[2];
533 u[3] = potW[5];
534
535 T.mv(u, Tu);
536
537 fluxW[0] = Tu[0];
538 potentialDiffW0 = Tu[0];
539
540 u[0] = potNw[0];
541 u[1] = potNw[1];
542 u[2] = potNw[2];
543 u[3] = potNw[5];
544
545 T.mv(u, Tu);
546
547 fluxNw[0] = Tu[0];
548 potentialDiffNw0 = Tu[0];
549 }
550 }
551
552 if (fIdx < 0 || fIdx == 1)
553 {
554 // calculate the flux through the subvolumeface 2 (subVolumeFaceIdx = 1)
555 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 1, 3, 0, 2, 5,
556 7);
557
558 if (caseL == 1)
559 {
560 u[0] = potW[1];
561 u[1] = potW[3];
562 u[2] = potW[0];
563 u[3] = potW[5];
564
565 T.mv(u, Tu);
566
567 fluxW[1] = Tu[0];
568 potentialDiffW1 = Tu[0];
569
570 u[0] = potNw[1];
571 u[1] = potNw[3];
572 u[2] = potNw[0];
573 u[3] = potNw[5];
574
575 T.mv(u, Tu);
576
577 fluxNw[1] = Tu[0];
578 potentialDiffNw1 = Tu[0];
579 }
580 else if (caseL == 2)
581 {
582 u[0] = potW[1];
583 u[1] = potW[3];
584 u[2] = potW[2];
585 u[3] = potW[7];
586
587 T.mv(u, Tu);
588
589 fluxW[1] = Tu[0];
590 potentialDiffW1 = Tu[0];
591
592 u[0] = potNw[1];
593 u[1] = potNw[3];
594 u[2] = potNw[2];
595 u[3] = potNw[7];
596
597 T.mv(u, Tu);
598
599 fluxNw[1] = Tu[0];
600 potentialDiffNw1 = Tu[0];
601 }
602 else if (caseL == 3)
603 {
604 u[0] = potW[1];
605 u[1] = potW[3];
606 u[2] = potW[2];
607 u[3] = potW[5];
608
609 T.mv(u, Tu);
610
611 fluxW[1] = Tu[0];
612 potentialDiffW1 = Tu[0];
613
614 u[0] = potNw[1];
615 u[1] = potNw[3];
616 u[2] = potNw[2];
617 u[3] = potNw[5];
618
619 T.mv(u, Tu);
620
621 fluxNw[1] = Tu[0];
622 potentialDiffNw1 = Tu[0];
623 }
624 else
625 {
626 u[0] = potW[1];
627 u[1] = potW[3];
628 u[2] = potW[0];
629 u[3] = potW[7];
630
631 T.mv(u, Tu);
632
633 fluxW[1] = Tu[0];
634 potentialDiffW1 = Tu[0];
635
636 u[0] = potNw[1];
637 u[1] = potNw[3];
638 u[2] = potNw[0];
639 u[3] = potNw[7];
640
641 T.mv(u, Tu);
642
643 fluxNw[1] = Tu[0];
644 potentialDiffNw1 = Tu[0];
645 }
646 }
647
648 if (fIdx < 0 || fIdx == 2)
649 {
650 // calculate the flux through the subvolumeface 3 (subVolumeFaceIdx = 2)
651 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 3, 2, 1, 0, 7,
652 6);
653
654 if (caseL == 1)
655 {
656 u[0] = potW[3];
657 u[1] = potW[2];
658 u[2] = potW[1];
659 u[3] = potW[7];
660
661 T.mv(u, Tu);
662
663 fluxW[2] = Tu[0];
664 potentialDiffW2 = Tu[0];
665
666 u[0] = potNw[3];
667 u[1] = potNw[2];
668 u[2] = potNw[1];
669 u[3] = potNw[7];
670
671 T.mv(u, Tu);
672
673 fluxNw[2] = Tu[0];
674 potentialDiffNw2 = Tu[0];
675 }
676 else if (caseL == 2)
677 {
678 u[0] = potW[3];
679 u[1] = potW[2];
680 u[2] = potW[0];
681 u[3] = potW[6];
682
683 T.mv(u, Tu);
684
685 fluxW[2] = Tu[0];
686 potentialDiffW2 = Tu[0];
687
688 u[0] = potNw[3];
689 u[1] = potNw[2];
690 u[2] = potNw[0];
691 u[3] = potNw[6];
692
693 T.mv(u, Tu);
694
695 fluxNw[2] = Tu[0];
696 potentialDiffNw2 = Tu[0];
697 }
698 else if (caseL == 3)
699 {
700 u[0] = potW[3];
701 u[1] = potW[2];
702 u[2] = potW[0];
703 u[3] = potW[7];
704
705 T.mv(u, Tu);
706
707 fluxW[2] = Tu[0];
708 potentialDiffW2 = Tu[0];
709
710 u[0] = potNw[3];
711 u[1] = potNw[2];
712 u[2] = potNw[0];
713 u[3] = potNw[7];
714
715 T.mv(u, Tu);
716
717 fluxNw[2] = Tu[0];
718 potentialDiffNw2 = Tu[0];
719 }
720 else
721 {
722 u[0] = potW[3];
723 u[1] = potW[2];
724 u[2] = potW[1];
725 u[3] = potW[6];
726
727 T.mv(u, Tu);
728
729 fluxW[2] = Tu[0];
730 potentialDiffW2 = Tu[0];
731
732 u[0] = potNw[3];
733 u[1] = potNw[2];
734 u[2] = potNw[1];
735 u[3] = potNw[6];
736
737 T.mv(u, Tu);
738
739 fluxNw[2] = Tu[0];
740 potentialDiffNw2 = Tu[0];
741 }
742 }
743
744 if (fIdx < 0 || fIdx == 3)
745 {
746 // calculate the flux through the subvolumeface 4 (subVolumeFaceIdx = 3)
747 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 2, 0, 3, 1, 6,
748 4);
749
750 if (caseL == 1)
751 {
752 u[0] = potW[2];
753 u[1] = potW[0];
754 u[2] = potW[3];
755 u[3] = potW[6];
756
757 T.mv(u, Tu);
758
759 fluxW[3] = Tu[0];
760 potentialDiffW3 = Tu[0];
761
762 u[0] = potNw[2];
763 u[1] = potNw[0];
764 u[2] = potNw[3];
765 u[3] = potNw[6];
766
767 T.mv(u, Tu);
768
769 fluxNw[3] = Tu[0];
770 potentialDiffNw3 = Tu[0];
771 }
772 else if (caseL == 2)
773 {
774 u[0] = potW[2];
775 u[1] = potW[0];
776 u[2] = potW[1];
777 u[3] = potW[4];
778
779 T.mv(u, Tu);
780
781 fluxW[3] = Tu[0];
782 potentialDiffW3 = Tu[0];
783
784 u[0] = potNw[2];
785 u[1] = potNw[0];
786 u[2] = potNw[1];
787 u[3] = potNw[4];
788
789 T.mv(u, Tu);
790
791 fluxNw[3] = Tu[0];
792 potentialDiffNw3 = Tu[0];
793 }
794 else if (caseL == 3)
795 {
796 u[0] = potW[2];
797 u[1] = potW[0];
798 u[2] = potW[1];
799 u[3] = potW[6];
800
801 T.mv(u, Tu);
802
803 fluxW[3] = Tu[0];
804 potentialDiffW3 = Tu[0];
805
806 u[0] = potNw[2];
807 u[1] = potNw[0];
808 u[2] = potNw[1];
809 u[3] = potNw[6];
810
811 T.mv(u, Tu);
812
813 fluxNw[3] = Tu[0];
814 potentialDiffNw3 = Tu[0];
815 }
816 else
817 {
818 u[0] = potW[2];
819 u[1] = potW[0];
820 u[2] = potW[3];
821 u[3] = potW[4];
822
823 T.mv(u, Tu);
824
825 fluxW[3] = Tu[0];
826 potentialDiffW3 = Tu[0];
827
828 u[0] = potNw[2];
829 u[1] = potNw[0];
830 u[2] = potNw[3];
831 u[3] = potNw[4];
832
833 T.mv(u, Tu);
834
835 fluxNw[3] = Tu[0];
836 potentialDiffNw3 = Tu[0];
837 }
838 }
839
840 if (fIdx < 0 || fIdx == 4)
841 {
842 // calculate the flux through the subvolumeface 5 (subVolumeFaceIdx = 4)
843 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 5, 4, 7, 6, 1,
844 0);
845
846 if (caseL == 1)
847 {
848 u[0] = potW[5];
849 u[1] = potW[4];
850 u[2] = potW[7];
851 u[3] = potW[1];
852
853 T.mv(u, Tu);
854
855 fluxW[4] = Tu[0];
856 potentialDiffW4 = Tu[0];
857
858 u[0] = potNw[5];
859 u[1] = potNw[4];
860 u[2] = potNw[7];
861 u[3] = potNw[1];
862
863 T.mv(u, Tu);
864
865 fluxNw[4] = Tu[0];
866 potentialDiffNw4 = Tu[0];
867 }
868 else if (caseL == 2)
869 {
870 u[0] = potW[5];
871 u[1] = potW[4];
872 u[2] = potW[6];
873 u[3] = potW[0];
874
875 T.mv(u, Tu);
876
877 fluxW[4] = Tu[0];
878 potentialDiffW4 = Tu[0];
879
880 u[0] = potNw[5];
881 u[1] = potNw[4];
882 u[2] = potNw[6];
883 u[3] = potNw[0];
884
885 T.mv(u, Tu);
886
887 fluxNw[4] = Tu[0];
888 potentialDiffNw4 = Tu[0];
889 }
890 else if (caseL == 3)
891 {
892 u[0] = potW[5];
893 u[1] = potW[4];
894 u[2] = potW[6];
895 u[3] = potW[1];
896
897 T.mv(u, Tu);
898
899 fluxW[4] = Tu[0];
900 potentialDiffW4 = Tu[0];
901
902 u[0] = potNw[5];
903 u[1] = potNw[4];
904 u[2] = potNw[6];
905 u[3] = potNw[1];
906
907 T.mv(u, Tu);
908
909 fluxNw[4] = Tu[0];
910 potentialDiffNw4 = Tu[0];
911 }
912 else
913 {
914 u[0] = potW[5];
915 u[1] = potW[4];
916 u[2] = potW[7];
917 u[3] = potW[0];
918
919 T.mv(u, Tu);
920
921 fluxW[4] = Tu[0];
922 potentialDiffW4 = Tu[0];
923
924 u[0] = potNw[5];
925 u[1] = potNw[4];
926 u[2] = potNw[7];
927 u[3] = potNw[0];
928
929 T.mv(u, Tu);
930
931 fluxNw[4] = Tu[0];
932 potentialDiffNw4 = Tu[0];
933 }
934 }
935
936 if (fIdx < 0 || fIdx == 5)
937 {
938 // calculate the flux through the subvolumeface 6 (subVolumeFaceIdx = 5)
939 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 7, 5, 6, 4, 3,
940 1);
941
942 if (caseL == 1)
943 {
944 u[0] = potW[7];
945 u[1] = potW[5];
946 u[2] = potW[6];
947 u[3] = potW[3];
948
949 T.mv(u, Tu);
950
951 fluxW[5] = Tu[0];
952 potentialDiffW5 = Tu[0];
953
954 u[0] = potNw[7];
955 u[1] = potNw[5];
956 u[2] = potNw[6];
957 u[3] = potNw[3];
958
959 T.mv(u, Tu);
960
961 fluxNw[5] = Tu[0];
962 potentialDiffNw5 = Tu[0];
963 }
964 else if (caseL == 2)
965 {
966 u[0] = potW[7];
967 u[1] = potW[5];
968 u[2] = potW[4];
969 u[3] = potW[1];
970
971 T.mv(u, Tu);
972
973 fluxW[5] = Tu[0];
974 potentialDiffW5 = Tu[0];
975
976 u[0] = potNw[7];
977 u[1] = potNw[5];
978 u[2] = potNw[4];
979 u[3] = potNw[1];
980
981 T.mv(u, Tu);
982
983 fluxNw[5] = Tu[0];
984 potentialDiffNw5 = Tu[0];
985 }
986 else if (caseL == 3)
987 {
988 u[0] = potW[7];
989 u[1] = potW[5];
990 u[2] = potW[4];
991 u[3] = potW[3];
992
993 T.mv(u, Tu);
994
995 fluxW[5] = Tu[0];
996 potentialDiffW5 = Tu[0];
997
998 u[0] = potNw[7];
999 u[1] = potNw[5];
1000 u[2] = potNw[4];
1001 u[3] = potNw[3];
1002
1003 T.mv(u, Tu);
1004
1005 fluxNw[5] = Tu[0];
1006 potentialDiffNw5 = Tu[0];
1007 }
1008 else
1009 {
1010 u[0] = potW[7];
1011 u[1] = potW[5];
1012 u[2] = potW[6];
1013 u[3] = potW[1];
1014
1015 T.mv(u, Tu);
1016
1017 fluxW[5] = Tu[0];
1018 potentialDiffW5 = Tu[0];
1019
1020 u[0] = potNw[7];
1021 u[1] = potNw[5];
1022 u[2] = potNw[6];
1023 u[3] = potNw[1];
1024
1025 T.mv(u, Tu);
1026
1027 fluxNw[5] = Tu[0];
1028 potentialDiffNw5 = Tu[0];
1029 }
1030 }
1031
1032 if (fIdx < 0 || fIdx == 6)
1033 {
1034 // calculate the flux through the subvolumeface 7 (subVolumeFaceIdx = 6)
1035 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 6, 7, 4, 5, 2,
1036 3);
1037
1038 if (caseL == 1)
1039 {
1040 u[0] = potW[6];
1041 u[1] = potW[7];
1042 u[2] = potW[4];
1043 u[3] = potW[2];
1044
1045 T.mv(u, Tu);
1046
1047 fluxW[6] = Tu[0];
1048 potentialDiffW6 = Tu[0];
1049
1050 u[0] = potNw[6];
1051 u[1] = potNw[7];
1052 u[2] = potNw[4];
1053 u[3] = potNw[2];
1054
1055 T.mv(u, Tu);
1056
1057 fluxNw[6] = Tu[0];
1058 potentialDiffNw6 = Tu[0];
1059 }
1060 else if (caseL == 2)
1061 {
1062 u[0] = potW[6];
1063 u[1] = potW[7];
1064 u[2] = potW[5];
1065 u[3] = potW[3];
1066
1067 T.mv(u, Tu);
1068
1069 fluxW[6] = Tu[0];
1070 potentialDiffW6 = Tu[0];
1071
1072 u[0] = potNw[6];
1073 u[1] = potNw[7];
1074 u[2] = potNw[5];
1075 u[3] = potNw[3];
1076
1077 T.mv(u, Tu);
1078
1079 fluxNw[6] = Tu[0];
1080 potentialDiffNw6 = Tu[0];
1081 }
1082 else if (caseL == 3)
1083 {
1084 u[0] = potW[6];
1085 u[1] = potW[7];
1086 u[2] = potW[5];
1087 u[3] = potW[2];
1088
1089 T.mv(u, Tu);
1090
1091 fluxW[6] = Tu[0];
1092 potentialDiffW6 = Tu[0];
1093
1094 u[0] = potNw[6];
1095 u[1] = potNw[7];
1096 u[2] = potNw[5];
1097 u[3] = potNw[2];
1098
1099 T.mv(u, Tu);
1100
1101 fluxNw[6] = Tu[0];
1102 potentialDiffNw6 = Tu[0];
1103 }
1104 else
1105 {
1106 u[0] = potW[6];
1107 u[1] = potW[7];
1108 u[2] = potW[4];
1109 u[3] = potW[3];
1110
1111 T.mv(u, Tu);
1112
1113 fluxW[6] = Tu[0];
1114 potentialDiffW6 = Tu[0];
1115
1116 u[0] = potNw[6];
1117 u[1] = potNw[7];
1118 u[2] = potNw[4];
1119 u[3] = potNw[3];
1120
1121 T.mv(u, Tu);
1122
1123 fluxNw[6] = Tu[0];
1124 potentialDiffNw6 = Tu[0];
1125 }
1126 }
1127
1128 if (fIdx < 0 || fIdx == 7)
1129 {
1130 // calculate the flux through the subvolumeface 8 (subVolumeFaceIdx = 7)
1131 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 4, 6, 5, 7, 0,
1132 2);
1133
1134 if (caseL == 1)
1135 {
1136 u[0] = potW[4];
1137 u[1] = potW[6];
1138 u[2] = potW[5];
1139 u[3] = potW[0];
1140
1141 T.mv(u, Tu);
1142
1143 fluxW[7] = Tu[0];
1144 potentialDiffW7 = Tu[0];
1145
1146 u[0] = potNw[4];
1147 u[1] = potNw[6];
1148 u[2] = potNw[5];
1149 u[3] = potNw[0];
1150
1151 T.mv(u, Tu);
1152
1153 fluxNw[7] = Tu[0];
1154 potentialDiffNw7 = Tu[0];
1155 }
1156 else if (caseL == 2)
1157 {
1158 u[0] = potW[4];
1159 u[1] = potW[6];
1160 u[2] = potW[7];
1161 u[3] = potW[2];
1162
1163 T.mv(u, Tu);
1164
1165 fluxW[7] = Tu[0];
1166 potentialDiffW7 = Tu[0];
1167
1168 u[0] = potNw[4];
1169 u[1] = potNw[6];
1170 u[2] = potNw[7];
1171 u[3] = potNw[2];
1172
1173 T.mv(u, Tu);
1174
1175 fluxNw[7] = Tu[0];
1176 potentialDiffNw7 = Tu[0];
1177 }
1178 else if (caseL == 3)
1179 {
1180 u[0] = potW[4];
1181 u[1] = potW[6];
1182 u[2] = potW[7];
1183 u[3] = potW[0];
1184
1185 T.mv(u, Tu);
1186
1187 fluxW[7] = Tu[0];
1188 potentialDiffW7 = Tu[0];
1189
1190 u[0] = potNw[4];
1191 u[1] = potNw[6];
1192 u[2] = potNw[7];
1193 u[3] = potNw[0];
1194
1195 T.mv(u, Tu);
1196
1197 fluxNw[7] = Tu[0];
1198 potentialDiffNw7 = Tu[0];
1199 }
1200 else
1201 {
1202 u[0] = potW[4];
1203 u[1] = potW[6];
1204 u[2] = potW[5];
1205 u[3] = potW[2];
1206
1207 T.mv(u, Tu);
1208
1209 fluxW[7] = Tu[0];
1210 potentialDiffW7 = Tu[0];
1211
1212 u[0] = potNw[4];
1213 u[1] = potNw[6];
1214 u[2] = potNw[5];
1215 u[3] = potNw[2];
1216
1217 T.mv(u, Tu);
1218
1219 fluxNw[7] = Tu[0];
1220 potentialDiffNw7 = Tu[0];
1221 }
1222 }
1223
1224 if (fIdx < 0 || fIdx == 8)
1225 {
1226 // calculate the flux through the subvolumeface 9 (subVolumeFaceIdx = 8)
1227 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 4, 0, 6, 2, 5,
1228 1);
1229
1230 if (caseL == 1)
1231 {
1232 u[0] = potW[4];
1233 u[1] = potW[0];
1234 u[2] = potW[6];
1235 u[3] = potW[5];
1236
1237 T.mv(u, Tu);
1238
1239 fluxW[8] = Tu[0];
1240 potentialDiffW8 = Tu[0];
1241
1242 u[0] = potNw[4];
1243 u[1] = potNw[0];
1244 u[2] = potNw[6];
1245 u[3] = potNw[5];
1246
1247 T.mv(u, Tu);
1248
1249 fluxNw[8] = Tu[0];
1250 potentialDiffNw8 = Tu[0];
1251 }
1252 else if (caseL == 2)
1253 {
1254 u[0] = potW[4];
1255 u[1] = potW[0];
1256 u[2] = potW[2];
1257 u[3] = potW[1];
1258
1259 T.mv(u, Tu);
1260
1261 fluxW[8] = Tu[0];
1262 potentialDiffW8 = Tu[0];
1263
1264 u[0] = potNw[4];
1265 u[1] = potNw[0];
1266 u[2] = potNw[2];
1267 u[3] = potNw[1];
1268
1269 T.mv(u, Tu);
1270
1271 fluxNw[8] = Tu[0];
1272 potentialDiffNw8 = Tu[0];
1273 }
1274 else if (caseL == 3)
1275 {
1276 u[0] = potW[4];
1277 u[1] = potW[0];
1278 u[2] = potW[2];
1279 u[3] = potW[5];
1280
1281 T.mv(u, Tu);
1282
1283 fluxW[8] = Tu[0];
1284 potentialDiffW8 = Tu[0];
1285
1286 u[0] = potNw[4];
1287 u[1] = potNw[0];
1288 u[2] = potNw[2];
1289 u[3] = potNw[5];
1290
1291 T.mv(u, Tu);
1292
1293 fluxNw[8] = Tu[0];
1294 potentialDiffNw8 = Tu[0];
1295 }
1296 else
1297 {
1298 u[0] = potW[4];
1299 u[1] = potW[0];
1300 u[2] = potW[6];
1301 u[3] = potW[1];
1302
1303 T.mv(u, Tu);
1304
1305 fluxW[8] = Tu[0];
1306 potentialDiffW8 = Tu[0];
1307
1308 u[0] = potNw[4];
1309 u[1] = potNw[0];
1310 u[2] = potNw[6];
1311 u[3] = potNw[1];
1312
1313 T.mv(u, Tu);
1314
1315 fluxNw[8] = Tu[0];
1316 potentialDiffNw8 = Tu[0];
1317 }
1318 }
1319
1320 if (fIdx < 0 || fIdx == 9)
1321 {
1322 // calculate the flux through the subvolumeface 10 (subVolumeFaceIdx = 9)
1323 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 1, 5, 3, 7, 0,
1324 4);
1325
1326 if (caseL == 1)
1327 {
1328 u[0] = potW[1];
1329 u[1] = potW[5];
1330 u[2] = potW[3];
1331 u[3] = potW[0];
1332
1333 T.mv(u, Tu);
1334
1335 fluxW[9] = Tu[0];
1336 potentialDiffW9 = Tu[0];
1337
1338 u[0] = potNw[1];
1339 u[1] = potNw[5];
1340 u[2] = potNw[3];
1341 u[3] = potNw[0];
1342
1343 T.mv(u, Tu);
1344
1345 fluxNw[9] = Tu[0];
1346 potentialDiffNw9 = Tu[0];
1347 }
1348 else if (caseL == 2)
1349 {
1350 u[0] = potW[1];
1351 u[1] = potW[5];
1352 u[2] = potW[7];
1353 u[3] = potW[4];
1354
1355 T.mv(u, Tu);
1356
1357 fluxW[9] = Tu[0];
1358 potentialDiffW9 = Tu[0];
1359
1360 u[0] = potNw[1];
1361 u[1] = potNw[5];
1362 u[2] = potNw[7];
1363 u[3] = potNw[4];
1364
1365 T.mv(u, Tu);
1366
1367 fluxNw[9] = Tu[0];
1368 potentialDiffNw9 = Tu[0];
1369 }
1370 else if (caseL == 3)
1371 {
1372 u[0] = potW[1];
1373 u[1] = potW[5];
1374 u[2] = potW[7];
1375 u[3] = potW[0];
1376
1377 T.mv(u, Tu);
1378
1379 fluxW[9] = Tu[0];
1380 potentialDiffW9 = Tu[0];
1381
1382 u[0] = potNw[1];
1383 u[1] = potNw[5];
1384 u[2] = potNw[7];
1385 u[3] = potNw[0];
1386
1387 T.mv(u, Tu);
1388
1389 fluxNw[9] = Tu[0];
1390 potentialDiffNw9 = Tu[0];
1391 }
1392 else
1393 {
1394 u[0] = potW[1];
1395 u[1] = potW[5];
1396 u[2] = potW[3];
1397 u[3] = potW[4];
1398
1399 T.mv(u, Tu);
1400
1401 fluxW[9] = Tu[0];
1402 potentialDiffW9 = Tu[0];
1403
1404 u[0] = potNw[1];
1405 u[1] = potNw[5];
1406 u[2] = potNw[3];
1407 u[3] = potNw[4];
1408
1409 T.mv(u, Tu);
1410
1411 fluxNw[9] = Tu[0];
1412 potentialDiffNw9 = Tu[0];
1413 }
1414 }
1415
1416 if (fIdx < 0 || fIdx == 10)
1417 {
1418 // calculate the flux through the subvolumeface 11 (subVolumeFaceIdx = 10)
1419 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 7, 3, 5, 1, 6,
1420 2);
1421
1422 if (caseL == 1)
1423 {
1424 u[0] = potW[7];
1425 u[1] = potW[3];
1426 u[2] = potW[5];
1427 u[3] = potW[6];
1428
1429 T.mv(u, Tu);
1430
1431 fluxW[10] = Tu[0];
1432 potentialDiffW10 = Tu[0];
1433
1434 u[0] = potNw[7];
1435 u[1] = potNw[3];
1436 u[2] = potNw[5];
1437 u[3] = potNw[6];
1438
1439 T.mv(u, Tu);
1440
1441 fluxNw[10] = Tu[0];
1442 potentialDiffNw10 = Tu[0];
1443 }
1444 else if (caseL == 2)
1445 {
1446 u[0] = potW[7];
1447 u[1] = potW[3];
1448 u[2] = potW[1];
1449 u[3] = potW[2];
1450
1451 T.mv(u, Tu);
1452
1453 fluxW[10] = Tu[0];
1454 potentialDiffW10 = Tu[0];
1455
1456 u[0] = potNw[7];
1457 u[1] = potNw[3];
1458 u[2] = potNw[1];
1459 u[3] = potNw[2];
1460
1461 T.mv(u, Tu);
1462
1463 fluxNw[10] = Tu[0];
1464 potentialDiffNw10 = Tu[0];
1465 }
1466 else if (caseL == 3)
1467 {
1468 u[0] = potW[7];
1469 u[1] = potW[3];
1470 u[2] = potW[1];
1471 u[3] = potW[6];
1472
1473 T.mv(u, Tu);
1474
1475 fluxW[10] = Tu[0];
1476 potentialDiffW10 = Tu[0];
1477
1478 u[0] = potNw[7];
1479 u[1] = potNw[3];
1480 u[2] = potNw[1];
1481 u[3] = potNw[6];
1482
1483 T.mv(u, Tu);
1484
1485 fluxNw[10] = Tu[0];
1486 potentialDiffNw10 = Tu[0];
1487 }
1488 else
1489 {
1490 u[0] = potW[7];
1491 u[1] = potW[3];
1492 u[2] = potW[5];
1493 u[3] = potW[2];
1494
1495 T.mv(u, Tu);
1496
1497 fluxW[10] = Tu[0];
1498 potentialDiffW10 = Tu[0];
1499
1500 u[0] = potNw[7];
1501 u[1] = potNw[3];
1502 u[2] = potNw[5];
1503 u[3] = potNw[2];
1504
1505 T.mv(u, Tu);
1506
1507 fluxNw[10] = Tu[0];
1508 potentialDiffNw10 = Tu[0];
1509 }
1510 }
1511
1512 if (fIdx < 0 || fIdx == 11)
1513 {
1514 // calculate the flux through the subvolumeface 12 (subVolumeFaceIdx = 11)
1515 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 2, 6, 0, 4, 3,
1516 7);
1517
1518 if (caseL == 1)
1519 {
1520 u[0] = potW[2];
1521 u[1] = potW[6];
1522 u[2] = potW[0];
1523 u[3] = potW[3];
1524
1525 T.mv(u, Tu);
1526
1527 fluxW[11] = Tu[0];
1528 potentialDiffW11 = Tu[0];
1529
1530 u[0] = potNw[2];
1531 u[1] = potNw[6];
1532 u[2] = potNw[0];
1533 u[3] = potNw[3];
1534
1535 T.mv(u, Tu);
1536
1537 fluxNw[11] = Tu[0];
1538 potentialDiffNw11 = Tu[0];
1539 }
1540 else if (caseL == 2)
1541 {
1542 u[0] = potW[2];
1543 u[1] = potW[6];
1544 u[2] = potW[4];
1545 u[3] = potW[7];
1546
1547 T.mv(u, Tu);
1548
1549 fluxW[11] = Tu[0];
1550 potentialDiffW11 = Tu[0];
1551
1552 u[0] = potNw[2];
1553 u[1] = potNw[6];
1554 u[2] = potNw[4];
1555 u[3] = potNw[7];
1556
1557 T.mv(u, Tu);
1558
1559 fluxNw[11] = Tu[0];
1560 potentialDiffNw11 = Tu[0];
1561 }
1562 else if (caseL == 3)
1563 {
1564 u[0] = potW[2];
1565 u[1] = potW[6];
1566 u[2] = potW[4];
1567 u[3] = potW[3];
1568
1569 T.mv(u, Tu);
1570
1571 fluxW[11] = Tu[0];
1572 potentialDiffW11 = Tu[0];
1573
1574 u[0] = potNw[2];
1575 u[1] = potNw[6];
1576 u[2] = potNw[4];
1577 u[3] = potNw[3];
1578
1579 T.mv(u, Tu);
1580
1581 fluxNw[11] = Tu[0];
1582 potentialDiffNw11 = Tu[0];
1583 }
1584 else
1585 {
1586 u[0] = potW[2];
1587 u[1] = potW[6];
1588 u[2] = potW[0];
1589 u[3] = potW[7];
1590
1591 T.mv(u, Tu);
1592
1593 fluxW[11] = Tu[0];
1594 potentialDiffW11 = Tu[0];
1595
1596 u[0] = potNw[2];
1597 u[1] = potNw[6];
1598 u[2] = potNw[0];
1599 u[3] = potNw[7];
1600
1601 T.mv(u, Tu);
1602
1603 fluxNw[11] = Tu[0];
1604 potentialDiffNw11 = Tu[0];
1605 }
1606 }
1607
1608 //store potentials for further calculations (saturation, ...)
1609 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffW0);
1610 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffNw0);
1611 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 1), -potentialDiffW3);
1612 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), -potentialDiffNw3);
1613 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 2), -potentialDiffW8);
1614 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 2), -potentialDiffNw8);
1615
1616 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 0), potentialDiffW1);
1617 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), potentialDiffNw1);
1618 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffW0);
1619 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffNw0);
1620 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 2), potentialDiffW9);
1621 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 2), potentialDiffNw9);
1622
1623 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffW3);
1624 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffNw3);
1625 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 1), -potentialDiffW2);
1626 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 1), -potentialDiffNw2);
1627 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 2), potentialDiffW11);
1628 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 2), potentialDiffNw11);
1629
1630 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 0), potentialDiffW2);
1631 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), potentialDiffNw2);
1632 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffW1);
1633 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffNw1);
1634 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 2), -potentialDiffW10);
1635 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 2), -potentialDiffNw10);
1636
1637 cellData5.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(4, 0), potentialDiffW8);
1638 cellData5.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(4, 0), potentialDiffNw8);
1639 cellData5.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(4, 1), -potentialDiffW4);
1640 cellData5.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(4, 1), -potentialDiffNw4);
1641 cellData5.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(4, 2), potentialDiffW7);
1642 cellData5.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(4, 2), potentialDiffNw7);
1643
1644 cellData6.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(5, 0), -potentialDiffW9);
1645 cellData6.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(5, 0), -potentialDiffNw9);
1646 cellData6.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(5, 1), -potentialDiffW5);
1647 cellData6.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(5, 1), -potentialDiffNw5);
1648 cellData6.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(5, 2), potentialDiffW4);
1649 cellData6.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(5, 2), potentialDiffNw4);
1650
1651 cellData7.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(6, 0), -potentialDiffW11);
1652 cellData7.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(6, 0), -potentialDiffNw11);
1653 cellData7.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(6, 1), -potentialDiffW7);
1654 cellData7.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(6, 1), -potentialDiffNw7);
1655 cellData7.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(6, 2), potentialDiffW6);
1656 cellData7.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(6, 2), potentialDiffNw6);
1657
1658 cellData8.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(7, 0), potentialDiffW10);
1659 cellData8.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(7, 0), potentialDiffNw10);
1660 cellData8.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(7, 1), -potentialDiffW6);
1661 cellData8.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(7, 1), -potentialDiffNw6);
1662 cellData8.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(7, 2), potentialDiffW5);
1663 cellData8.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(7, 2), potentialDiffNw5);
1664
1665 // compute mobilities of subvolumeface 1 (subVolumeFaceIdx = 0)
1666 Dune::FieldVector<Scalar, numPhases> lambda0Upw(0.0);
1667 lambda0Upw[wPhaseIdx] = (potentialDiffW0 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
1668 lambda0Upw[nPhaseIdx] = (potentialDiffNw0 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
1669
1670 // compute mobilities of subvolumeface 2 (subVolumeFaceIdx = 1)
1671 Dune::FieldVector<Scalar, numPhases> lambda1Upw(0.0);
1672 lambda1Upw[wPhaseIdx] = (potentialDiffW1 >= 0) ? lambda2[wPhaseIdx] : lambda4[wPhaseIdx];
1673 lambda1Upw[nPhaseIdx] = (potentialDiffNw1 >= 0) ? lambda2[nPhaseIdx] : lambda4[nPhaseIdx];
1674
1675 // compute mobilities of subvolumeface 3 (subVolumeFaceIdx = 2)
1676 Dune::FieldVector<Scalar, numPhases> lambda2Upw(0.0);
1677 lambda2Upw[wPhaseIdx] = (potentialDiffW2 >= 0) ? lambda4[wPhaseIdx] : lambda3[wPhaseIdx];
1678 lambda2Upw[nPhaseIdx] = (potentialDiffNw2 >= 0) ? lambda4[nPhaseIdx] : lambda3[nPhaseIdx];
1679
1680 // compute mobilities of subvolumeface 4 (subVolumeFaceIdx = 3)
1681 Dune::FieldVector<Scalar, numPhases> lambda3Upw(0.0);
1682 lambda3Upw[wPhaseIdx] = (potentialDiffW3 >= 0) ? lambda3[wPhaseIdx] : lambda1[wPhaseIdx];
1683 lambda3Upw[nPhaseIdx] = (potentialDiffNw3 >= 0) ? lambda3[nPhaseIdx] : lambda1[nPhaseIdx];
1684
1685 // compute mobilities of subvolumeface 5 (subVolumeFaceIdx = 4)
1686 Dune::FieldVector<Scalar, numPhases> lambda4Upw(0.0);
1687 lambda4Upw[wPhaseIdx] = (potentialDiffW4 >= 0) ? lambda6[wPhaseIdx] : lambda5[wPhaseIdx];
1688 lambda4Upw[nPhaseIdx] = (potentialDiffNw4 >= 0) ? lambda6[nPhaseIdx] : lambda5[nPhaseIdx];
1689
1690 // compute mobilities of subvolumeface 6 (subVolumeFaceIdx = 5)
1691 Dune::FieldVector<Scalar, numPhases> lambda5Upw(0.0);
1692 lambda5Upw[wPhaseIdx] = (potentialDiffW5 >= 0) ? lambda8[wPhaseIdx] : lambda6[wPhaseIdx];
1693 lambda5Upw[nPhaseIdx] = (potentialDiffNw5 >= 0) ? lambda8[nPhaseIdx] : lambda6[nPhaseIdx];
1694
1695 // compute mobilities of subvolumeface 7 (subVolumeFaceIdx = 6)
1696 Dune::FieldVector<Scalar, numPhases> lambda6Upw(0.0);
1697 lambda6Upw[wPhaseIdx] = (potentialDiffW6 >= 0) ? lambda7[wPhaseIdx] : lambda8[wPhaseIdx];
1698 lambda6Upw[nPhaseIdx] = (potentialDiffNw6 >= 0) ? lambda7[nPhaseIdx] : lambda8[nPhaseIdx];
1699
1700 // compute mobilities of subvolumeface 8 (subVolumeFaceIdx = 7)
1701 Dune::FieldVector<Scalar, numPhases> lambda7Upw(0.0);
1702 lambda7Upw[wPhaseIdx] = (potentialDiffW7 >= 0) ? lambda5[wPhaseIdx] : lambda7[wPhaseIdx];
1703 lambda7Upw[nPhaseIdx] = (potentialDiffNw7 >= 0) ? lambda5[nPhaseIdx] : lambda7[nPhaseIdx];
1704
1705 // compute mobilities of subvolumeface 9 (subVolumeFaceIdx = 8)
1706 Dune::FieldVector<Scalar, numPhases> lambda8Upw(0.0);
1707 lambda8Upw[wPhaseIdx] = (potentialDiffW8 >= 0) ? lambda5[wPhaseIdx] : lambda1[wPhaseIdx];
1708 lambda8Upw[nPhaseIdx] = (potentialDiffNw8 >= 0) ? lambda5[nPhaseIdx] : lambda1[nPhaseIdx];
1709
1710 // compute mobilities of subvolumeface 10 (subVolumeFaceIdx = 9)
1711 Dune::FieldVector<Scalar, numPhases> lambda9Upw(0.0);
1712 lambda9Upw[wPhaseIdx] = (potentialDiffW9 >= 0) ? lambda2[wPhaseIdx] : lambda6[wPhaseIdx];
1713 lambda9Upw[nPhaseIdx] = (potentialDiffNw9 >= 0) ? lambda2[nPhaseIdx] : lambda6[nPhaseIdx];
1714
1715 // compute mobilities of subvolumeface 11 (subVolumeFaceIdx = 10)
1716 Dune::FieldVector<Scalar, numPhases> lambda10Upw(0.0);
1717 lambda10Upw[wPhaseIdx] = (potentialDiffW10 >= 0) ? lambda8[wPhaseIdx] : lambda4[wPhaseIdx];
1718 lambda10Upw[nPhaseIdx] = (potentialDiffNw10 >= 0) ? lambda8[nPhaseIdx] : lambda4[nPhaseIdx];
1719
1720 // compute mobilities of subvolumeface 12 (subVolumeFaceIdx = 11)
1721 Dune::FieldVector<Scalar, numPhases> lambda11Upw(0.0);
1722 lambda11Upw[wPhaseIdx] = (potentialDiffW11 >= 0) ? lambda3[wPhaseIdx] : lambda7[wPhaseIdx];
1723 lambda11Upw[nPhaseIdx] = (potentialDiffNw11 >= 0) ? lambda3[nPhaseIdx] : lambda7[nPhaseIdx];
1724
1725 for (int i = 0; i < numPhases; i++)
1726 {
1727 // evaluate parts of velocity --> always take the normal for which the flux is calculated!
1728 DimVector vel12 = interactionVolume.getNormal(0, 0);
1729 DimVector vel21 = interactionVolume.getNormal(1, 1);
1730 DimVector vel24 = interactionVolume.getNormal(1, 0);
1731 DimVector vel42 = interactionVolume.getNormal(3, 1);
1732 DimVector vel43 = interactionVolume.getNormal(3, 0);
1733 DimVector vel34 = interactionVolume.getNormal(2, 1);
1734 DimVector vel31 = interactionVolume.getNormal(2, 0);
1735 DimVector vel13 = interactionVolume.getNormal(0, 1);
1736 DimVector vel65 = interactionVolume.getNormal(5, 2);
1737 DimVector vel56 = interactionVolume.getNormal(4, 1);
1738 DimVector vel86 = interactionVolume.getNormal(7, 2);
1739 DimVector vel68 = interactionVolume.getNormal(5, 1);
1740 DimVector vel78 = interactionVolume.getNormal(6, 2);
1741 DimVector vel87 = interactionVolume.getNormal(7, 1);
1742 DimVector vel57 = interactionVolume.getNormal(4, 2);
1743 DimVector vel75 = interactionVolume.getNormal(6, 1);
1744 DimVector vel51 = interactionVolume.getNormal(4, 0);
1745 DimVector vel15 = interactionVolume.getNormal(0, 2);
1746 DimVector vel26 = interactionVolume.getNormal(1, 2);
1747 DimVector vel62 = interactionVolume.getNormal(5, 0);
1748 DimVector vel84 = interactionVolume.getNormal(7, 0);
1749 DimVector vel48 = interactionVolume.getNormal(3, 2);
1750 DimVector vel37 = interactionVolume.getNormal(2, 2);
1751 DimVector vel73 = interactionVolume.getNormal(6, 0);
1752 vel21 *= -1;
1753 vel42 *= -1;
1754 vel34 *= -1;
1755 vel13 *= -1;
1756 vel56 *= -1;
1757 vel68 *= -1;
1758 vel87 *= -1;
1759 vel75 *= -1;
1760 vel15 *= -1;
1761 vel62 *= -1;
1762 vel48 *= -1;
1763 vel73 *= -1;
1764
1765 Dune::FieldVector<Scalar, 12> flux(0);
1766 switch (i)
1767 {
1768 case wPhaseIdx:
1769 {
1770 flux = fluxW;
1771 break;
1772 }
1773 case nPhaseIdx:
1774 {
1775 flux = fluxNw;
1776 break;
1777 }
1778 }
1779
1780 vel12 *= flux[0] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal1, 0, 0));
1781 vel21 *= flux[0] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal2, 1, 1));
1782 vel24 *= flux[1] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal2, 1, 0));
1783 vel42 *= flux[1] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal4, 3, 1));
1784 vel43 *= flux[2] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal4, 3, 0));
1785 vel34 *= flux[2] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal3, 2, 1));
1786 vel31 *= flux[3] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal3, 2, 0));
1787 vel13 *= flux[3] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal1, 0, 1));
1788 vel65 *= flux[4] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal6, 5, 2));
1789 vel56 *= flux[4] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal5, 4, 1));
1790 vel86 *= flux[5] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal8, 7, 2));
1791 vel68 *= flux[5] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal6, 5, 1));
1792 vel78 *= flux[6] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal7, 6, 2));
1793 vel87 *= flux[6] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal8, 7, 1));
1794 vel57 *= flux[7] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal5, 4, 2));
1795 vel75 *= flux[7] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal7, 6, 1));
1796 vel51 *= flux[8] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal5, 4, 0));
1797 vel15 *= flux[8] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal1, 0, 2));
1798 vel26 *= flux[9] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal2, 1, 2));
1799 vel62 *= flux[9] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal6, 5, 0));
1800 vel84 *= flux[10] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal8, 7, 0));
1801 vel48 *= flux[10] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal4, 3, 2));
1802 vel37 *= flux[11] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal3, 2, 2));
1803 vel73 *= flux[11] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal7, 6, 0));
1804
1805
1806 Scalar lambdaT0 = lambda0Upw[wPhaseIdx] + lambda0Upw[nPhaseIdx];
1807 Scalar lambdaT1 = lambda1Upw[wPhaseIdx] + lambda1Upw[nPhaseIdx];
1808 Scalar lambdaT2 = lambda2Upw[wPhaseIdx] + lambda2Upw[nPhaseIdx];
1809 Scalar lambdaT3 = lambda3Upw[wPhaseIdx] + lambda3Upw[nPhaseIdx];
1810 Scalar lambdaT4 = lambda4Upw[wPhaseIdx] + lambda4Upw[nPhaseIdx];
1811 Scalar lambdaT5 = lambda5Upw[wPhaseIdx] + lambda5Upw[nPhaseIdx];
1812 Scalar lambdaT6 = lambda6Upw[wPhaseIdx] + lambda6Upw[nPhaseIdx];
1813 Scalar lambdaT7 = lambda7Upw[wPhaseIdx] + lambda7Upw[nPhaseIdx];
1814 Scalar lambdaT8 = lambda8Upw[wPhaseIdx] + lambda8Upw[nPhaseIdx];
1815 Scalar lambdaT9 = lambda9Upw[wPhaseIdx] + lambda9Upw[nPhaseIdx];
1816 Scalar lambdaT10 = lambda10Upw[wPhaseIdx] + lambda10Upw[nPhaseIdx];
1817 Scalar lambdaT11 = lambda11Upw[wPhaseIdx] + lambda11Upw[nPhaseIdx];
1818 Scalar fracFlow0 = (lambdaT0 > threshold_) ? lambda0Upw[i] / (lambdaT0) : 0.0;
1819 Scalar fracFlow1 = (lambdaT1 > threshold_) ? lambda1Upw[i] / (lambdaT1) : 0.0;
1820 Scalar fracFlow2 = (lambdaT2 > threshold_) ? lambda2Upw[i] / (lambdaT2) : 0.0;
1821 Scalar fracFlow3 = (lambdaT3 > threshold_) ? lambda3Upw[i] / (lambdaT3) : 0.0;
1822 Scalar fracFlow4 = (lambdaT4 > threshold_) ? lambda4Upw[i] / (lambdaT4) : 0.0;
1823 Scalar fracFlow5 = (lambdaT5 > threshold_) ? lambda5Upw[i] / (lambdaT5) : 0.0;
1824 Scalar fracFlow6 = (lambdaT6 > threshold_) ? lambda6Upw[i] / (lambdaT6) : 0.0;
1825 Scalar fracFlow7 = (lambdaT7 > threshold_) ? lambda7Upw[i] / (lambdaT7) : 0.0;
1826 Scalar fracFlow8 = (lambdaT8 > threshold_) ? lambda8Upw[i] / (lambdaT8) : 0.0;
1827 Scalar fracFlow9 = (lambdaT9 > threshold_) ? lambda9Upw[i] / (lambdaT9) : 0.0;
1828 Scalar fracFlow10 = (lambdaT10 > threshold_) ? lambda10Upw[i] / (lambdaT10) : 0.0;
1829 Scalar fracFlow11 = (lambdaT11 > threshold_) ? lambda11Upw[i] / (lambdaT11) : 0.0;
1830
1831 vel12 *= fracFlow0;
1832 vel21 *= fracFlow0;
1833 vel24 *= fracFlow1;
1834 vel42 *= fracFlow1;
1835 vel43 *= fracFlow2;
1836 vel34 *= fracFlow2;
1837 vel31 *= fracFlow3;
1838 vel13 *= fracFlow3;
1839 vel65 *= fracFlow4;
1840 vel56 *= fracFlow4;
1841 vel86 *= fracFlow5;
1842 vel68 *= fracFlow5;
1843 vel78 *= fracFlow6;
1844 vel87 *= fracFlow6;
1845 vel57 *= fracFlow7;
1846 vel75 *= fracFlow7;
1847 vel51 *= fracFlow8;
1848 vel15 *= fracFlow8;
1849 vel26 *= fracFlow9;
1850 vel62 *= fracFlow9;
1851 vel84 *= fracFlow10;
1852 vel48 *= fracFlow10;
1853 vel37 *= fracFlow11;
1854 vel73 *= fracFlow11;
1855
1856 //store velocities
1857 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 0), vel12);
1858 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 1), vel13);
1859 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 2), vel15);
1860 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 0), vel24);
1861 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 1), vel21);
1862 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 2), vel26);
1863 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 0), vel31);
1864 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 1), vel34);
1865 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 2), vel37);
1866 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 0), vel43);
1867 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 1), vel42);
1868 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 2), vel48);
1869 cellData5.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(4, 0), vel51);
1870 cellData5.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(4, 1), vel56);
1871 cellData5.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(4, 2), vel57);
1872 cellData6.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(5, 0), vel62);
1873 cellData6.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(5, 1), vel68);
1874 cellData6.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(5, 2), vel65);
1875 cellData7.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(6, 0), vel73);
1876 cellData7.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(6, 1), vel75);
1877 cellData7.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(6, 2), vel78);
1878 cellData8.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(7, 0), vel84);
1879 cellData8.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(7, 1), vel87);
1880 cellData8.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(7, 2), vel86);
1881 }
1882 //set velocity marker
1883 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 0));
1884 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 1));
1885 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 2));
1886 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 0));
1887 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 1));
1888 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 2));
1889 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 0));
1890 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 1));
1891 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 2));
1892 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 0));
1893 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
1894 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 2));
1895 cellData5.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(4, 0));
1896 cellData5.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(4, 1));
1897 cellData5.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(4, 2));
1898 cellData6.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(5, 0));
1899 cellData6.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(5, 1));
1900 cellData6.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(5, 2));
1901 cellData7.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(6, 0));
1902 cellData7.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(6, 1));
1903 cellData7.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(6, 2));
1904 cellData8.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(7, 0));
1905 cellData8.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(7, 1));
1906 cellData8.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(7, 2));
1907 }
1908
1918template<class TypeTag>
1919void FvMpfaL3dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData, int elemIdx)
1920{
1921 auto element = interactionVolume.getSubVolumeElement(elemIdx);
1922
1923 // get global coordinate of cell centers
1924 const GlobalPosition& globalPos = element.geometry().center();
1925
1926 // permeability vector at boundary
1927 DimMatrix permeability(problem_.spatialParams().intrinsicPermeability(element));
1928
1929 //get mobilities of the phases
1930 Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
1931 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
1932
1933 Scalar potW = cellData.potential(wPhaseIdx);
1934 Scalar potNw = cellData.potential(nPhaseIdx);
1935
1936 for (int fIdx = 0; fIdx < dim; fIdx++)
1937 {
1938 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
1939
1940 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
1941 {
1942 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressureEqIdx))
1943 {
1944 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
1945
1946 const GlobalPosition& globalPosFace = interactionVolume.getFacePosition(elemIdx, fIdx);
1947
1948 DimVector distVec(globalPosFace - globalPos);
1949 Scalar dist = distVec.two_norm();
1950 DimVector& normal = interactionVolume.getNormal(elemIdx, fIdx);
1951
1952 // get pc and lambda at the boundary
1953 Scalar satWBound = cellData.saturation(wPhaseIdx);
1954 //check boundary sat at face 1
1955 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
1956 {
1957 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
1958 switch (saturationType_)
1959 {
1960 case sw:
1961 {
1962 satWBound = satBound;
1963 break;
1964 }
1965 case sn:
1966 {
1967 satWBound = 1 - satBound;
1968 break;
1969 }
1970 }
1971
1972 }
1973
1974 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
1975 Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
1976
1977 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
1978 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1979
1980 pcBound += gravityDiffBound;
1981
1982 Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
1983 lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
1984 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
1985 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
1986
1987 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
1988 Scalar potentialBoundW = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx] + density_[wPhaseIdx]*gdeltaZ;
1989 Scalar potentialBoundNw = potentialBoundW;
1990
1991 //calculate potential gradients
1992 switch (pressureType_)
1993 {
1994 case pw:
1995 {
1996 potentialBoundNw += pcBound;
1997 break;
1998 }
1999 case pn:
2000 {
2001 //calculate potential gradients
2002 potentialBoundW -= pcBound;
2003 break;
2004 }
2005 }
2006
2007 Scalar potentialDiffW = (potW - potentialBoundW) / dist;
2008 Scalar potentialDiffNw = (potNw - potentialBoundNw) / dist;
2009
2010 //store potentials for further calculations (saturation, ...)
2011 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, potentialDiffW);
2012 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, potentialDiffNw);
2013
2014 //calculated phase velocities from advective velocities -> capillary pressure velocity already added in pressure part!
2015 DimVector velocityW(0);
2016 DimVector velocityNw(0);
2017
2018 // calculate capillary pressure gradient
2019 DimVector potentialGradient = normal;
2020 potentialGradient *= (potW - potentialBoundW) / dist;
2021 permeability.mv(potentialGradient, velocityW);
2022
2023 potentialGradient = normal;
2024 potentialGradient *= (potNw - potentialBoundNw) / dist;
2025 permeability.mv(potentialGradient, velocityNw);
2026
2027 velocityW *= (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
2028 velocityNw *= (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
2029
2030 //velocity is calculated from two vertices of one intersection!
2031 velocityW *= 0.25;
2032 velocityNw *= 0.25;
2033
2034 //store velocities
2035 cellData.fluxData().addVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
2036 cellData.fluxData().addVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
2037 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
2038 }
2039 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx))
2040 {
2041 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
2042
2043 DimVector& normal = interactionVolume.getNormal(elemIdx, fIdx);
2044
2045 // get neumann boundary value
2046 PrimaryVariables boundValues(interactionVolume.getNeumannValues(intVolFaceIdx));
2047
2048 boundValues[wPhaseIdx] /= density_[wPhaseIdx];
2049 boundValues[nPhaseIdx] /= density_[nPhaseIdx];
2050
2051 DimVector velocityW(normal);
2052 DimVector velocityNw(normal);
2053
2054 velocityW *= boundValues[wPhaseIdx] / (4.0*interactionVolume.getFaceArea(elemIdx, fIdx));
2055 velocityNw *= boundValues[nPhaseIdx]
2056 / (4.0*interactionVolume.getFaceArea(elemIdx, fIdx));
2057
2058 //store potentials for further calculations (saturation, ...)
2059 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, boundValues[wPhaseIdx]);
2060 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, boundValues[nPhaseIdx]);
2061
2062 //store velocities
2063 cellData.fluxData().addVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
2064 cellData.fluxData().addVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
2065 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
2066 }
2067 else
2068 {
2069 std::cout << "interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx)"
2070 << interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx) << "\n";
2071 DUNE_THROW(Dune::NotImplemented,
2072 "No valid boundary condition type defined for pressure equation!");
2073 }
2074 }
2075 }
2076}
2077
2078} // end namespace Dumux
2079#endif
3dtransmissibilitycalculator.hh
Provides methods for transmissibility calculation in 3-d.
Dumux::normal
Vector normal(const Vector &v)
Create a vector normal to the given one (v is expected to be non-zero)
Definition: normal.hh:36
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::FvMpfaL3dTransmissibilityCalculator
Provides methods for transmissibility calculation in 3-d.
Definition: 3dtransmissibilitycalculator.hh:51
Dumux::FvMpfaL3dTransmissibilityCalculator::TransmissibilityType
Dune::FieldMatrix< Scalar, dim, 2 *dim - dim+1 > TransmissibilityType
Type of the transmissibility matrix.
Definition: 3dtransmissibilitycalculator.hh:78
Dumux::FvMpfaL3dTransmissibilityCalculator::transmissibility
int transmissibility(Dune::FieldMatrix< Scalar, dim, 2 *dim-dim+1 > &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx3, int idx4, int idx5, int idx6)
Dumux::FvMpfaL3dVelocity2p
Class for calculating 3-d velocities from cell-wise constant pressure values.
Definition: 3dvelocity.hh:55
Dumux::FvMpfaL3dVelocity2p::FvMpfaL3dVelocity2p
FvMpfaL3dVelocity2p(Problem &problem)
Constructs a FvMpfaL3dVelocity2p object.
Definition: 3dvelocity.hh:137
Dumux::FvMpfaL3dVelocity2p::initialize
void initialize(bool solveTwice=true)
Initializes the velocity model.
Definition: 3dvelocity.hh:158
Dumux::FvMpfaL3dVelocity2p::addOutputVtkFields
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: 3dvelocity.hh:186
Dumux::FvMpfaL3dVelocity2p::calculateInnerInteractionVolumeVelocity
void calculateInnerInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData1, CellData &cellData2, CellData &cellData3, CellData &cellData4, CellData &cellData5, CellData &cellData6, CellData &cellData7, CellData &cellData8, InteractionVolumeContainer &interactionVolumes, TransmissibilityCalculator &transmissibilityCalculator, int fIdx=-1)
Calculate velocities for flux faces of an interaction volume.
Definition: 3dvelocity.hh:294
Dumux::FvMpfaL3dVelocity2p::calculateBoundaryInteractionVolumeVelocity
void calculateBoundaryInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData, int elemIdx)
Calculates the velocity at a boundary flux faces.
Definition: 3dvelocity.hh:1919
properties.hh
Specifies the properties for immiscible 2p diffusion/pressure models.
properties.hh
Properties for a MPFA method.
Generated by doxygen 1.9.3