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