3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
2dpressureadaptive.hh
Go to the documentation of this file.
1/*****************************************************************************
2 * See the file COPYING for full copying permissions. *
3 * *
4 * This program is free software: you can redistribute it and/or modify *
5 * it under the terms of the GNU General Public License as published by *
6 * the Free Software Foundation, either version 3 of the License, or *
7 * (at your option) any later version. *
8 * *
9 * This program is distributed in the hope that it will be useful, *
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
12 * GNU General Public License for more details. *
13 * *
14 * You should have received a copy of the GNU General Public License *
15 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
16 *****************************************************************************/
23#ifndef DUMUX_FVMPFAL2DPRESSURE2P_ADAPTIVE_HH
24#define DUMUX_FVMPFAL2DPRESSURE2P_ADAPTIVE_HH
25
26// dumux environment
27#include <dumux/porousmediumflow/sequential/cellcentered/pressure.hh>
28#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh>
29#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
30#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
31#include "2dtransmissibilitycalculator.hh"
32
33#include <dumux/common/deprecated.hh>
34
35namespace Dumux {
36
75template<class TypeTag>
76class FvMpfaL2dPressure2pAdaptive: public FVPressure<TypeTag>
77{
78 using ParentType = FVPressure<TypeTag>;
79 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
80
81 enum
82 {
83 dim = GridView::dimension, dimWorld = GridView::dimensionworld
84 };
85
86 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
87 using Problem = GetPropType<TypeTag, Properties::Problem>;
88
89 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
90
91 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
92
93 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
94 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
95
96 using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
97 using CellData = GetPropType<TypeTag, Properties::CellData>;
98
99 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
100
101 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
102 using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
103
104 using GridTypeIndices = GetPropType<TypeTag, Properties::GridTypeIndices>;
105
106 enum
107 {
108 pw = Indices::pressureW,
109 pn = Indices::pressureNw,
110 sw = Indices::saturationW,
111 sn = Indices::saturationNw
112 };
113 enum
114 {
115 wPhaseIdx = Indices::wPhaseIdx,
116 nPhaseIdx = Indices::nPhaseIdx,
117 pressureIdx = Indices::pressureIdx,
118 saturationIdx = Indices::saturationIdx,
119 pressEqIdx = Indices::pressureEqIdx,
120 satEqIdx = Indices::satEqIdx,
121 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
122 };
123 enum
124 {
125 globalCorner = 2,
126 globalEdge = 3,
127 neumannNeumann = 0,
128 dirichletDirichlet = 1,
129 dirichletNeumann = 2,
130 neumannDirichlet = 3
131 };
132
133
134 using Element = typename GridView::Traits::template Codim<0>::Entity;
135 using IntersectionIterator = typename GridView::IntersectionIterator;
136 using Intersection = typename GridView::Intersection;
137 using Grid = typename GridView::Grid;
138
139 using LocalPosition = Dune::FieldVector<Scalar, dim>;
140 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
141 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
142
143 using DimVector = Dune::FieldVector<Scalar, dim>;
144
145public:
146
147
154 using InteractionVolume = FVMPFALInteractionVolume<TypeTag>;
155 using TransmissibilityCalculator = FvMpfaL2dTransmissibilityCalculator<TypeTag>;
156private:
157
158 using GlobalInteractionVolumeVector = std::vector<InteractionVolume>;
159 using InnerBoundaryVolumeFaces = std::vector<Dune::FieldVector<bool, 2*dim> >;
160
162 Intersection getNextIntersection_(const Element&, const IntersectionIterator&);
163
165 friend class FVPressure<TypeTag>;
166 void initializeMatrix();
167
168 void storeInteractionVolumeInfo();
169
170 void printInteractionVolumes();
171
173 void assemble();
174
175public:
176
178 void updateMaterialLaws();
179
185 void updateInteractionVolumeInfo()
186 {
187 interactionVolumes_.clear();
188 innerBoundaryVolumeFaces_.clear();
189
190 interactionVolumes_.resize(problem_.gridView().size(dim), InteractionVolume(problem_.gridView().grid()));
191 innerBoundaryVolumeFaces_.resize(problem_.gridView().size(0), Dune::FieldVector<bool, 2 * dim>(false));
192
193 storeInteractionVolumeInfo();
194// printInteractionVolumes();
195 }
196
202 void initialize()
203 {
204 ParentType::initialize();
205
206 const auto element = *problem_.gridView().template begin<0>();
207 FluidState fluidState;
208 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
209 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
210 fluidState.setTemperature(problem_.temperature(element));
211 fluidState.setSaturation(wPhaseIdx, 1.);
212 fluidState.setSaturation(nPhaseIdx, 0.);
213 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
214 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
215 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
216 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
217
218 updateMaterialLaws();
219
220 updateInteractionVolumeInfo();
221
222 assemble();
223 this->solve();
224
225 storePressureSolution();
226
227 return;
228 }
229
233 void storePressureSolution()
234 {
235 // iterate through leaf grid an evaluate c0 at cell center
236 for (const auto& element : elements(problem_.gridView()))
237 {
238 storePressureSolution(element);
239 }
240 }
241
247 void storePressureSolution(const Element& element)
248 {
249 int eIdxGlobal = problem_.variables().index(element);
250 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
251
252 switch (pressureType_)
253 {
254 case pw:
255 {
256 Scalar potW = this->pressure()[eIdxGlobal];
257
258 Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
259 Scalar potPc = cellData.capillaryPressure()
260 + gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
261
262 cellData.setPotential(wPhaseIdx, potW);
263 cellData.setPotential(nPhaseIdx, potW + potPc);
264
265 Scalar pressW = potW - gravityDiff * density_[wPhaseIdx];
266
267 cellData.setPressure(wPhaseIdx, pressW);
268 cellData.setPressure(nPhaseIdx, pressW + cellData.capillaryPressure());
269
270 break;
271 }
272 case pn:
273 {
274 Scalar potNw = this->pressure()[eIdxGlobal];
275
276 Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
277 Scalar potPc = cellData.capillaryPressure()
278 + gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
279
280 cellData.setPotential(nPhaseIdx, potNw);
281 cellData.setPotential(wPhaseIdx, potNw - potPc);
282
283 Scalar pressNw = potNw - gravityDiff * density_[nPhaseIdx];
284
285 cellData.setPressure(wPhaseIdx, pressNw - cellData.capillaryPressure());
286 cellData.setPressure(nPhaseIdx, pressNw);
287
288 break;
289 }
290 }
291
292 cellData.fluxData().resetVelocity();
293 }
294
300 void update()
301 {
302 int gridSize = problem_.gridView().size(0);
303
304 //error bounds for error term for incompressible models to correct unphysical saturation
305 //over/undershoots due to saturation transport
306 timeStep_ = problem_.timeManager().timeStepSize();
307 maxError_ = 0.0;
308 for (int i = 0; i < gridSize; i++)
309 {
310 Scalar sat = 0;
311 using std::max;
312 switch (saturationType_)
313 {
314 case sw:
315 sat = problem_.variables().cellData(i).saturation(wPhaseIdx);
316 break;
317 case sn:
318 sat = problem_.variables().cellData(i).saturation(nPhaseIdx);
319 break;
320 }
321 if (sat > 1.0)
322 {
323 maxError_ = max(maxError_, (sat - 1.0) / timeStep_);
324 }
325 if (sat < 0.0)
326 {
327 maxError_ = max(maxError_, (-sat) / timeStep_);
328 }
329 }
330
331 if (problem_.gridAdapt().wasAdapted())
332 {
333 // update RHS vector, matrix
334 this->A_.setSize(gridSize, gridSize); //
335 this->f_.resize(gridSize);
336 this->pressure().resize(gridSize);
337
338 for (int i = 0; i < gridSize; i++)
339 {
340 CellData& cellData = problem_.variables().cellData(i);
341
342 switch (pressureType_)
343 {
344 case pw:
345 this->pressure()[i] = cellData.pressure(wPhaseIdx);
346 break;
347 case pn:
348 this->pressure()[i] = cellData.pressure(nPhaseIdx);
349 break;
350 }
351 }
352
353 initializeMatrix();
354 updateInteractionVolumeInfo();
355 }
356 assemble();
357 this->solve();
358
359
360 storePressureSolution();
361
362 return;
363 }
364
375 template<class MultiWriter>
376 void addOutputVtkFields(MultiWriter &writer)
377 {
378 int size = problem_.gridView().size(0);
379 ScalarSolutionType *potential = writer.allocateManagedBuffer(size);
380
381 (*potential) = this->pressure();
382
383 if (pressureType_ == pw)
384 {
385 writer.attachCellData(*potential, "wetting potential");
386 }
387
388 if (pressureType_ == pn)
389 {
390 writer.attachCellData(*potential, "nonwetting potential");
391 }
392
393 if (vtkOutputLevel_ > 0)
394 {
395 ScalarSolutionType *pressure = writer.allocateManagedBuffer(size);
396 ScalarSolutionType *pressureSecond = writer.allocateManagedBuffer(size);
397 ScalarSolutionType *potentialSecond = writer.allocateManagedBuffer(size);
398 ScalarSolutionType *pc = writer.allocateManagedBuffer(size);
399
400 for (const auto& element : elements(problem_.gridView()))
401 {
402 int idx = problem_.variables().index(element);
403 CellData& cellData = problem_.variables().cellData(idx);
404
405 (*pc)[idx] = cellData.capillaryPressure();
406
407 if (pressureType_ == pw)
408 {
409 (*pressure)[idx] = cellData.pressure(wPhaseIdx);
410 (*potentialSecond)[idx] = cellData.potential(nPhaseIdx);
411 (*pressureSecond)[idx] = cellData.pressure(nPhaseIdx);
412 }
413
414 if (pressureType_ == pn)
415 {
416 (*pressure)[idx] = cellData.pressure(nPhaseIdx);
417 (*potentialSecond)[idx] = cellData.potential(wPhaseIdx);
418 (*pressureSecond)[idx] = cellData.pressure(wPhaseIdx);
419 }
420 }
421
422 if (pressureType_ == pw)
423 {
424 writer.attachCellData(*pressure, "wetting pressure");
425 writer.attachCellData(*pressureSecond, "nonwetting pressure");
426 writer.attachCellData(*potentialSecond, "nonwetting potential");
427 }
428
429 if (pressureType_ == pn)
430 {
431 writer.attachCellData(*pressure, "nonwetting pressure");
432 writer.attachCellData(*pressureSecond, "wetting pressure");
433 writer.attachCellData(*potentialSecond, "wetting potential");
434 }
435
436 writer.attachCellData(*pc, "capillary pressure");
437 }
438
439 return;
440 }
441
446 FvMpfaL2dPressure2pAdaptive(Problem& problem) :
447 ParentType(problem), problem_(problem), transmissibilityCalculator_(problem),
448 gravity_(problem.gravity()),
449 maxError_(0.), timeStep_(1.)
450 {
451 if (pressureType_ != pw && pressureType_ != pn)
452 {
453 DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
454 }
455 if (saturationType_ != sw && saturationType_ != sn)
456 {
457 DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
458 }
459 // The methods to calculate transmissibility are also used by the compressible
460 // compositional models (2p2c), although this implementation does not support it.
461 // Hence warning is not thrown if the file is not used from compositional modules
462 // that require 2p2cproperties.hh.
463 #ifndef DUMUX_2P2CPROPERTIES_HH
464 if (getPropValue<TypeTag, Properties::EnableCompressibility>())
465 {
466 DUNE_THROW(Dune::NotImplemented, "Compressibility not supported!");
467 }
468 #endif
469 if (dim != 2)
470 {
471 DUNE_THROW(Dune::NotImplemented, "Dimension not supported!");
472 }
473
474 ErrorTermFactor_ = getParam<Scalar>("Impet.ErrorTermFactor");
475 ErrorTermLowerBound_ = getParam<Scalar>("Impet.ErrorTermLowerBound");
476 ErrorTermUpperBound_ = getParam<Scalar>("Impet.ErrorTermUpperBound");
477
478 density_[wPhaseIdx] = 0.;
479 density_[nPhaseIdx] = 0.;
480 viscosity_[wPhaseIdx] = 0.;
481 viscosity_[nPhaseIdx] = 0.;
482
483 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
484 }
485
486private:
487 Problem& problem_;
488 TransmissibilityCalculator transmissibilityCalculator_;
489
490protected:
491 GlobalInteractionVolumeVector interactionVolumes_;
492 InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_;
493
494private:
495 const Dune::FieldVector<Scalar, dimWorld>& gravity_;
496
497 Scalar maxError_;
498 Scalar timeStep_;
499 Scalar ErrorTermFactor_;
500 Scalar ErrorTermLowerBound_;
501 Scalar ErrorTermUpperBound_;
502
503 Scalar density_[numPhases];
504 Scalar viscosity_[numPhases];
505
506 int vtkOutputLevel_;
507
508 static constexpr Scalar threshold_ = 1e-15;
510 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
512 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
514 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
515
516 // TODO doc me!
517 Scalar evaluateErrorTerm_(CellData& cellData)
518 {
519 //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
520 // error reduction routine: volumetric error is damped and inserted to right hand side
521 Scalar sat = 0;
522 switch (saturationType_)
523 {
524 case sw:
525 sat = cellData.saturation(wPhaseIdx);
526 break;
527 case sn:
528 sat = cellData.saturation(nPhaseIdx);
529 break;
530 }
531
532 Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
533 if (sat < 0.0)
534 {
535 error = sat;
536 }
537 error /= timeStep_;
538
539 using std::abs;
540 Scalar errorAbs = abs(error);
541
542 if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
543 && (!problem_.timeManager().willBeFinished()))
544 {
545 return ErrorTermFactor_ * error;
546 }
547 return 0.0;
548 }
549
550};
551
552// TODO doc me!
553template<class TypeTag>
554typename FvMpfaL2dPressure2pAdaptive<TypeTag>::Intersection
555 FvMpfaL2dPressure2pAdaptive<TypeTag>::getNextIntersection_(const Element& element,
556 const IntersectionIterator& isIt)
557{
558 auto isItBegin = problem_.gridView().ibegin(element);
559 const auto isEndIt = problem_.gridView().iend(element);
560
561 auto tempIsIt = isIt;
562 auto nextIsIt = ++tempIsIt;
563
564 // get 'nextIsIt'
565 switch (getPropValue<TypeTag, Properties::GridImplementation>())
566 {
567 // for ALUGrid and UGGrid
568 case GridTypeIndices::aluGrid:
569 case GridTypeIndices::ugGrid:
570 {
571 if (nextIsIt == isEndIt)
572 nextIsIt = isItBegin;
573
574 break;
575 }
576 default:
577 {
578 DUNE_THROW(Dune::NotImplemented,
579 "GridType can not be used with adaptive MPFAL!");
580 break;
581 }
582 }
583
584 return *nextIsIt;
585}
586
587// TODO doc me!
588template<class TypeTag>
589void FvMpfaL2dPressure2pAdaptive<TypeTag>::initializeMatrix()
590{
591 // determine matrix row sizes
592 for (const auto& element : elements(problem_.gridView()))
593 {
594 // cell index
595 int eIdxGlobalI = problem_.variables().index(element);
596
597 // initialize row size
598 int rowSize = 1;
599
600 // run through all intersections with neighbors
601 const auto isEndIt = problem_.gridView().iend(element);
602 for (auto isIt = problem_.gridView().ibegin(element); isIt != isEndIt; ++isIt)
603 {
604 const auto& intersection = *isIt;
605
606 if (intersection.neighbor())
607 {
608 rowSize++;
609
610 auto nextIntersection = getNextIntersection_(element, isIt);
611
612 if (nextIntersection.neighbor())
613 {
614 bool isCorner = true;
615
616 for (const auto& innerIntersection
617 : intersections(problem_.gridView(), intersection.outside()))
618 {
619 for (const auto& innerNextIntersection
620 : intersections(problem_.gridView(), nextIntersection.outside()))
621 {
622 if (innerIntersection.neighbor() && innerNextIntersection.neighbor())
623 {
624 if (innerIntersection.outside() == nextIntersection.outside()
625 || innerNextIntersection.outside() == intersection.outside())
626 {
627 isCorner = false;
628 break;
629 }
630 }
631 }
632 if (!isCorner)
633 {
634 break;
635 }
636 }
637 if (isCorner)
638 {
639 rowSize++;
640 }
641 }
642 }
643
644 } // end of intersection loop
645
646 using std::min;
647 rowSize = min(rowSize, 13); //in 2-D
648
649 // set number of indices in row eIdxGlobalI to rowSize
650 this->A_.setrowsize(eIdxGlobalI, rowSize);
651
652 } // end of element loop
653
654 // indicate that size of all rows is defined
655 this->A_.endrowsizes();
656 // determine position of matrix entries
657 for (const auto& element : elements(problem_.gridView()))
658 {
659 // cell index
660 int eIdxGlobalI = problem_.variables().index(element);
661
662 // add diagonal index
663 this->A_.addindex(eIdxGlobalI, eIdxGlobalI);
664
665 // run through all intersections with neighbors
666 const auto isEndIt = problem_.gridView().iend(element);
667 for (auto isIt = problem_.gridView().ibegin(element); isIt != isEndIt; ++isIt)
668 {
669 const auto& intersection = *isIt;
670
671 if (intersection.neighbor())
672 {
673 // access neighbor
674 auto outside = intersection.outside();
675 int eIdxGlobalJ = problem_.variables().index(outside);
676
677 // add off diagonal index
678 // add index (row,col) to the matrix
679 this->A_.addindex(eIdxGlobalI, eIdxGlobalJ);
680
681 if (element.level() < outside.level())
682 {
683 continue;
684 }
685
686 auto nextIntersection = getNextIntersection_(element, isIt);
687
688 if (nextIntersection.neighbor())
689 {
690 // access the common neighbor of intersection's and nextIntersection's outside
691 if (element.level() < nextIntersection.outside().level())
692 {
693 continue;
694 }
695
696 for (const auto& innerIntersection
697 : intersections(problem_.gridView(), intersection.outside()))
698 {
699 for (const auto& innerNextIntersection
700 : intersections(problem_.gridView(), nextIntersection.outside()))
701 {
702 if (innerIntersection.neighbor() && innerNextIntersection.neighbor())
703 {
704 auto innerOutside = innerIntersection.outside();
705 auto nextOutside = nextIntersection.outside();
706
707 if (innerOutside == innerNextIntersection.outside() && innerOutside != element
708 && innerOutside != nextOutside)
709 {
710 int eIdxGlobalCorner = problem_.variables().index(innerOutside);
711
712 this->A_.addindex(eIdxGlobalI, eIdxGlobalCorner);
713
714 if (element.level() > outside.level())
715 {
716 int eIdxGlobalJCorner = problem_.variables().index(nextOutside);
717
718 this->A_.addindex(eIdxGlobalJ, eIdxGlobalJCorner);
719 }
720 if (element.level() > nextOutside.level())
721 {
722 int eIdxGlobalJCorner = problem_.variables().index(nextOutside);
723
724 this->A_.addindex(eIdxGlobalJCorner, eIdxGlobalJ);
725 }
726 if (element.level() > innerOutside.level())
727 {
728 this->A_.addindex(eIdxGlobalCorner, eIdxGlobalI);
729 }
730 }
731 }
732 }
733 }
734 }
735 }
736 } // end of intersection loop
737 } // end of element loop
738
739 // indicate that all indices are defined, check consistency
740 this->A_.endindices();
741
742 return;
743}
744// Indices used in a interaction volume of the MPFA-o method
745// ___________________________________________________
746// | | |
747// | nuxy: cell geometry | nxy: face normal |
748// | vectors (see MPFA) | |
749// | | |
750// | 4-----------3-----------3 |
751// | | --> nu43 | nu34 <-- | |
752// | | |nu41 1|--> n43 ||nu32 |
753// | | v ^ |0 ^ v| |
754// |____________4__0__|n14__|__n23_|_1__2____________|
755// | | 1 0 | 0 |
756// | | ^ |1 nu23 ^ | |
757// | | |nu14 0|--> n12 | | |
758// | | -->nu12 | nu21<-- | |
759// | 1-----------1-----------2 |
760// | elementnumber |inter- |
761// | |face- |
762// | |number |
763// |________________________|________________________|
764
765// only for 2-D general quadrilateral
766// TODO doc me!
767template<class TypeTag>
768void FvMpfaL2dPressure2pAdaptive<TypeTag>::storeInteractionVolumeInfo()
769{
770 BoundaryTypes bcType;
771
772 // run through all elements
773 for (const auto& element : elements(problem_.gridView()))
774 {
775 // get index
776 int eIdxGlobal1 = problem_.variables().index(element);
777
778 const auto refElement = referenceElement(element);
779
780 const auto isEndIt12 = problem_.gridView().iend(element);
781 for (auto isIt12 = problem_.gridView().ibegin(element); isIt12 != isEndIt12; ++isIt12)
782 {
783 const auto& intersection12 = *isIt12;
784
785 int indexInInside12 = intersection12.indexInInside();
786
787 if (intersection12.neighbor())
788 {
789 // access neighbor cell 2 of 'intersection12'
790 auto element2 = intersection12.outside();
791 int eIdxGlobal2 = problem_.variables().index(element2);
792
793 if (element.level() < element2.level())
794 {
795 continue;
796 }
797
798 auto intersection14 = getNextIntersection_(element, isIt12);
799
800 int indexInInside14 = intersection14.indexInInside();
801
802 //get center vertex
803 GlobalPosition corner1234(0);
804 int globalVertIdx1234 = 0;
805 bool finished = false;
806 // get the global coordinate and global vertex index of corner1234
807 for (int i = 0; i < intersection12.geometry().corners(); ++i)
808 {
809 int localVertIdx12corner = refElement.subEntity(indexInInside12, dim - 1, i, dim);
810
811 int globalVertIdx12corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx12corner));
812
813 for (int j = 0; j < intersection14.geometry().corners(); ++j)
814 {
815 int localVertIdx14corner = refElement.subEntity(indexInInside14, dim - 1, j, dim);
816
817 int globalVertIdx14corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx14corner));
818
819 if (globalVertIdx12corner == globalVertIdx14corner)
820 {
821 corner1234 = element.geometry().corner(localVertIdx12corner);
822
823 globalVertIdx1234 = globalVertIdx12corner;
824
825 finished = true;
826 break;
827 }
828 }
829
830 if (finished)
831 {
832 break;
833 }
834 }
835
836 if (interactionVolumes_[globalVertIdx1234].isStored() || !finished)
837 {
838 continue;
839 }
840 else
841 {
842 interactionVolumes_[globalVertIdx1234].setStored();
843 }
844
845 interactionVolumes_[globalVertIdx1234].setCenterPosition(corner1234);
846
847 //store pointer 1
848 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element, 0);
849 interactionVolumes_[globalVertIdx1234].setIndexOnElement(indexInInside12, 0, 0);
850 interactionVolumes_[globalVertIdx1234].setIndexOnElement(indexInInside14, 0, 1);
851
852 // center of face in global coordinates, i.e., the midpoint of edge 'intersection12'
853 const GlobalPosition& globalPosFace12 = intersection12.geometry().center();
854
855 // get face volume
856 Scalar faceVol12 = intersection12.geometry().volume() / 2.0;
857
858 // get outer normal vector scaled with half volume of face 'intersection12'
859 DimVector unitOuterNormal12 = intersection12.centerUnitOuterNormal();
860
861 // center of face in global coordinates, i.e., the midpoint of edge 'intersection14'
862 GlobalPosition globalPosFace41 = intersection14.geometry().center();
863
864 // get face volume
865 Scalar faceVol41 = intersection14.geometry().volume() / 2.0;
866
867 // get outer normal vector scaled with half volume of face 'intersection14': for numbering of n see Aavatsmark, Eigestad
868 DimVector unitOuterNormal14 = intersection14.centerUnitOuterNormal();
869
870 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal12, 0, 0);
871 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 0, 1);
872 //get the normals of from cell 2 and 4
873 unitOuterNormal14 *= -1;
874 unitOuterNormal12 *= -1;
875 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol12, 0, 0);
876 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 0, 1);
877 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace12, 0, 0);
878 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace41, 0, 1);
879
880 //store pointer 2
881 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element2, 1);
882 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection12.indexInOutside(), 1, 1);
883 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal12, 1, 1);
884 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol12, 1, 1);
885 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace12, 1, 1);
886
887 // 'intersection14' is an interior face
888 if (intersection14.neighbor())
889 {
890 // neighbor cell 3
891 // access neighbor cell 3
892 auto element4 = intersection14.outside();
893
894 //store pointer 4
895 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element4, 3);
896 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection14.indexInOutside(), 3, 0);
897
898 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 3, 0);
899 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 3, 0);
900 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace41, 3, 0);
901
902 // cell 3
903 GlobalPosition globalPosFace23(0);
904 GlobalPosition globalPosFace34(0);
905
906 if (element4.level() < element.level())
907 {
908 bool isHangingNode = false;
909
910 for (const auto& intersection2
911 : intersections(problem_.gridView(), element2))
912 {
913 bool breakLoop = false;
914 for (const auto& intersection4
915 : intersections(problem_.gridView(), element4))
916 {
917 if (intersection2.neighbor() && intersection4.neighbor())
918 {
919 auto element32 = intersection2.outside();
920 auto element34 = intersection4.outside();
921
922 //hanging node!
923 if (element32 == element4)
924 {
925 if (element.level() != element2.level())
926 {
927 breakLoop = true;
928 isHangingNode = false;
929 break;
930 }
931
932 isHangingNode = true;
933
934 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInInside(),
935 1, 0);
936 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInOutside(),
937 3, 1);
938
939 globalPosFace23 = intersection2.geometry().center();
940 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace23, 1, 0);
941 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace23, 3, 1);
942
943 Scalar faceVol23 = intersection2.geometry().volume() / 2.0;
944 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 1, 0);
945 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 3, 1);
946
947 // get outer normal vector scaled with half volume of face : for numbering of n see Aavatsmark, Eigestad
948 DimVector unitOuterNormal23 = intersection2.centerUnitOuterNormal();
949
950 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 1, 0);
951 unitOuterNormal23 *= -1;
952 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 3, 1);
953 }
954 else if (element34 == element2)
955 {
956 if (element.level() != element2.level())
957 {
958 breakLoop = true;
959 isHangingNode = false;
960 break;
961 }
962
963 isHangingNode = true;
964
965 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInInside(),
966 3, 1);
967 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInOutside(),
968 1, 0);
969
970 // get outer normal vector scaled with half volume of face : for numbering of n see Aavatsmark, Eigestad
971 DimVector unitOuterNormal43 = intersection4.centerUnitOuterNormal();
972 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 3, 1);
973 }
974 }
975 }
976 if (breakLoop)
977 {
978 break;
979 }
980 }
981 if (!isHangingNode)
982 {
983 bool regularNode = false;
984
985 for (const auto& intersection2
986 : intersections(problem_.gridView(), element2))
987 {
988 for (const auto& intersection4
989 : intersections(problem_.gridView(), element4))
990 {
991 if (intersection4.neighbor())
992 {
993 auto element41 = intersection4.outside();
994
995 if (element41 == element && element41.level() > element.level())
996 {
997 //adjust values of intersection12 in case of hanging nodes
998 globalPosFace41 = intersection4.geometry().center();
999 faceVol41 = intersection4.geometry().volume() / 2.0;
1000
1001 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 0, 1);
1002 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace41, 0,
1003 1);
1004 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 3, 0);
1005 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace41, 3,
1006 0);
1007 }
1008 }
1009
1010 if (intersection2.neighbor() && intersection4.neighbor())
1011 {
1012 auto element32 = intersection2.outside();
1013 auto element34 = intersection4.outside();
1014
1015 //hanging node!
1016 if (element32 == element34 && element32 != element)
1017 {
1018 //store pointer 3
1019 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element32,
1020 2);
1021
1022 interactionVolumes_[globalVertIdx1234].setIndexOnElement(
1023 intersection2.indexInInside(), 1, 0);
1024 interactionVolumes_[globalVertIdx1234].setIndexOnElement(
1025 intersection2.indexInOutside(), 2, 1);
1026 interactionVolumes_[globalVertIdx1234].setIndexOnElement(
1027 intersection4.indexInInside(), 3, 1);
1028 interactionVolumes_[globalVertIdx1234].setIndexOnElement(
1029 intersection4.indexInOutside(), 2, 0);
1030
1031 globalPosFace23 = intersection2.geometry().center();
1032 globalPosFace34 = intersection4.geometry().center();
1033
1034 Scalar faceVol23 = intersection2.geometry().volume() / 2.0;
1035 Scalar faceVol34 = intersection4.geometry().volume() / 2.0;
1036
1037 // get outer normal vector scaled with half volume of face : for numbering of n see Aavatsmark, Eigestad
1038 DimVector unitOuterNormal23 = intersection2.centerUnitOuterNormal();
1039 DimVector unitOuterNormal43 = intersection4.centerUnitOuterNormal();
1040
1041 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 1, 0);
1042 unitOuterNormal23 *= -1;
1043 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 2, 1);
1044 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 3, 1);
1045 unitOuterNormal43 *= -1;
1046 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 2, 0);
1047
1048 if (element32.level() > element2.level())
1049 {
1050 for (const auto& intersection3
1051 : intersections(problem_.gridView(), element32))
1052 {
1053 if (intersection3.neighbor())
1054 {
1055 auto element23 = intersection3.outside();
1056
1057 if (element23 == element2)
1058 {
1059 globalPosFace23 = intersection3.geometry().center();
1060 faceVol23 = intersection3.geometry().volume() / 2.0;
1061
1062 }
1063 }
1064 }
1065 }
1066 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 1, 0);
1067 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 2, 1);
1068 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace23, 1,
1069 0);
1070 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace23, 2,
1071 1);
1072
1073 if (element34.level() > element4.level())
1074 {
1075 for (const auto& intersection3
1076 : intersections(problem_.gridView(), element34))
1077 {
1078 if (intersection3.neighbor())
1079 {
1080 auto element43 = intersection3.outside();
1081
1082 if (element43 == element4)
1083 {
1084 globalPosFace34 = intersection3.geometry().center();
1085 faceVol34 = intersection3.geometry().volume() / 2.0;
1086
1087 }
1088 }
1089 }
1090 }
1091
1092 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 2, 0);
1093 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 3, 1);
1094 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace34, 2,
1095 0);
1096 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace34, 3,
1097 1);
1098
1099 regularNode = true;
1100 break;
1101 }
1102 }
1103 if (regularNode)
1104 {
1105 break;
1106 }
1107 }
1108 }
1109 if (!regularNode)
1110 {
1111 interactionVolumes_[globalVertIdx1234].reset();
1112 continue;
1113 }
1114 }
1115 }
1116 else
1117 {
1118 bool regularNode = false;
1119 for (const auto& intersection2
1120 : intersections(problem_.gridView(), element2))
1121 {
1122 for (const auto& intersection4
1123 : intersections(problem_.gridView(), element4))
1124 {
1125 if (intersection4.neighbor())
1126 {
1127 auto element41 = intersection4.outside();
1128
1129 if (element41 == element && element41.level() > element.level())
1130 {
1131 //adjust values of intersection12 in case of hanging nodes
1132 globalPosFace41 = intersection4.geometry().center();
1133 faceVol41 = intersection4.geometry().volume() / 2.0;
1134
1135 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 0, 1);
1136 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace41, 0, 1);
1137 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 3, 0);
1138 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace41, 3, 0);
1139 }
1140 }
1141
1142 if (intersection2.neighbor() && intersection4.neighbor())
1143 {
1144 auto element32 = intersection2.outside();
1145 auto element34 = intersection4.outside();
1146
1147 // find the common neighbor cell between cell 2 and cell 3, except cell 1
1148 if (element32 == element34 && element32 != element)
1149 {
1150 //store pointer 3
1151 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element32, 2);
1152
1153 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInInside(),
1154 1, 0);
1155 interactionVolumes_[globalVertIdx1234].setIndexOnElement(
1156 intersection2.indexInOutside(), 2, 1);
1157 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInInside(),
1158 3, 1);
1159 interactionVolumes_[globalVertIdx1234].setIndexOnElement(
1160 intersection4.indexInOutside(), 2, 0);
1161
1162 globalPosFace23 = intersection2.geometry().center();
1163 globalPosFace34 = intersection4.geometry().center();
1164
1165 Scalar faceVol23 = intersection2.geometry().volume() / 2.0;
1166 Scalar faceVol34 = intersection4.geometry().volume() / 2.0;
1167
1168 // get outer normal vector scaled with half volume of face : for numbering of n see Aavatsmark, Eigestad
1169 DimVector unitOuterNormal23 = intersection2.centerUnitOuterNormal();
1170
1171 DimVector unitOuterNormal43 = intersection4.centerUnitOuterNormal();
1172
1173 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 1, 0);
1174 unitOuterNormal23 *= -1;
1175 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 2, 1);
1176 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 3, 1);
1177 unitOuterNormal43 *= -1;
1178 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 2, 0);
1179 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 1, 0);
1180 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 2, 1);
1181 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 2, 0);
1182 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 3, 1);
1183 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace23, 1, 0);
1184 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace23, 2, 1);
1185 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace34, 2, 0);
1186 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace34, 3, 1);
1187
1188 regularNode = true;
1189 break;
1190 }
1191 }
1192 }
1193 if (regularNode)
1194 {
1195 break;
1196 }
1197 }
1198 if (!regularNode)
1199 {
1200 interactionVolumes_[globalVertIdx1234].reset();
1201 continue;
1202 }
1203 }
1204
1205 }
1206
1207 // 'intersection14' is on the boundary
1208 else
1209 {
1210 problem_.boundaryTypes(bcType, intersection14);
1211 PrimaryVariables boundValues(0.0);
1212
1213 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 3);
1214 if (bcType.isNeumann(pressEqIdx))
1215 {
1216 problem_.neumann(boundValues, intersection14);
1217 boundValues *= faceVol41;
1218 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 3);
1219 }
1220 if (bcType.hasDirichlet())
1221 {
1222 problem_.dirichlet(boundValues, intersection14);
1223 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 3);
1224 }
1225
1226 // get common geometry information for the following computation
1227 // get the information of the face 'isIt23' between cell2 and cell4 (locally numbered)
1228
1229 // center of face in global coordinates, i.e., the midpoint of edge 'isIt23'
1230 GlobalPosition globalPosFace23(0);
1231
1232 // get face volume
1233 Scalar faceVol23 = 0;
1234
1235 // get outer normal vector scaled with half volume of face 'isIt23'
1236 DimVector unitOuterNormal23(0);
1237
1238 finished = false;
1239
1240 for (const auto& intersection2
1241 : intersections(problem_.gridView(), element2))
1242 {
1243 if (intersection2.boundary())
1244 {
1245 for (int i = 0; i < intersection2.geometry().corners(); ++i)
1246 {
1247 int localVertIdx2corner = refElement.subEntity(intersection2.indexInInside(), dim - 1, i,
1248 dim);
1249
1250 int globalVertIdx2corner = problem_.variables().index(element2.template subEntity<dim>(localVertIdx2corner));
1251
1252 if (globalVertIdx2corner == globalVertIdx1234)
1253 {
1254 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInInside(), 1,
1255 0);
1256
1257 globalPosFace23 = intersection2.geometry().center();
1258
1259 faceVol23 = intersection2.geometry().volume() / 2.0;
1260
1261 unitOuterNormal23 = intersection2.centerUnitOuterNormal();
1262
1263 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 1, 0);
1264 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 1, 0);
1265 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace23, 1, 0);
1266
1267 problem_.boundaryTypes(bcType, intersection2);
1268 boundValues = 0.0;
1269
1270 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 1);
1271 if (bcType.isNeumann(pressEqIdx))
1272 {
1273 problem_.neumann(boundValues, intersection2);
1274 boundValues *= faceVol23;
1275 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 1);
1276 }
1277 if (bcType.hasDirichlet())
1278 {
1279 problem_.dirichlet(boundValues, intersection2);
1280 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 1);
1281 }
1282
1283 interactionVolumes_[globalVertIdx1234].setOutsideFace(2);
1284
1285
1286 if (element.level() == element2.level())
1287 {
1288 innerBoundaryVolumeFaces_[eIdxGlobal1][intersection12.indexInInside()] = true;
1289 innerBoundaryVolumeFaces_[eIdxGlobal2][intersection12.indexInOutside()] = true;
1290 }
1291 else if (element.level() < element2.level())
1292 {
1293 innerBoundaryVolumeFaces_[eIdxGlobal2][intersection12.indexInOutside()] = true;
1294 }
1295 else if (element.level() > element2.level())
1296 {
1297 innerBoundaryVolumeFaces_[eIdxGlobal1][intersection12.indexInInside()] = true;
1298 }
1299
1300 finished = true;
1301
1302 break;
1303 }
1304 }
1305 }
1306 if (finished)
1307 {
1308 break;
1309 }
1310 }
1311 if (!finished)
1312 {
1313 DUNE_THROW(
1314 Dune::NotImplemented,
1315 "fvmpfao2pfaboundpressure2p.hh, l. 997: boundary shape not available as interaction volume shape");
1316 }
1317 }
1318 }
1319
1320 // handle boundary face 'intersection12'
1321 else
1322 {
1323 auto intersection14 = getNextIntersection_(element, isIt12);
1324
1325 int indexInInside14 = intersection14.indexInInside();
1326
1327 //get center vertex
1328 GlobalPosition corner1234(0);
1329 int globalVertIdx1234 = 0;
1330
1331 // get the global coordinate and global vertex index of corner1234
1332 for (int i = 0; i < intersection12.geometry().corners(); ++i)
1333 {
1334 bool finished = false;
1335
1336 int localVertIdx12corner = refElement.subEntity(indexInInside12, dim - 1, i, dim);
1337
1338 int globalVertIdx12corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx12corner));
1339
1340 for (int j = 0; j < intersection14.geometry().corners(); ++j)
1341 {
1342 int localVertIdx14corner = refElement.subEntity(indexInInside14, dim - 1, j, dim);
1343
1344 int globalVertIdx14corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx14corner));
1345
1346 if (globalVertIdx12corner == globalVertIdx14corner)
1347 {
1348 corner1234 = element.geometry().corner(localVertIdx12corner);
1349
1350 globalVertIdx1234 = globalVertIdx12corner;
1351
1352 finished = true;
1353 break;
1354 }
1355 }
1356
1357 if (finished)
1358 {
1359 break;
1360 }
1361 }
1362
1363 if (interactionVolumes_[globalVertIdx1234].isStored())
1364 {
1365 continue;
1366 }
1367 else
1368 {
1369 interactionVolumes_[globalVertIdx1234].setStored();
1370 }
1371
1372 interactionVolumes_[globalVertIdx1234].setCenterPosition(corner1234);
1373
1374 //store pointer 1
1375 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element, 0);
1376 interactionVolumes_[globalVertIdx1234].setIndexOnElement(indexInInside12, 0, 0);
1377 interactionVolumes_[globalVertIdx1234].setIndexOnElement(indexInInside14, 0, 1);
1378
1379 // center of face in global coordinates, i.e., the midpoint of edge 'intersection12'
1380 const GlobalPosition& globalPosFace12 = intersection12.geometry().center();
1381
1382 // get face volume
1383 Scalar faceVol12 = intersection12.geometry().volume() / 2.0;
1384
1385 // get outer normal vector scaled with half volume of face 'intersection12'
1386 DimVector unitOuterNormal12 = intersection12.centerUnitOuterNormal();
1387
1388 // center of face in global coordinates, i.e., the midpoint of edge 'intersection14'
1389 const GlobalPosition& globalPosFace41 = intersection14.geometry().center();
1390
1391 // get face volume
1392 Scalar faceVol41 = intersection14.geometry().volume() / 2.0;
1393
1394 // get outer normal vector scaled with half volume of face 'intersection14': for numbering of n see Aavatsmark, Eigestad
1395 DimVector unitOuterNormal14 = intersection14.centerUnitOuterNormal();
1396
1397 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal12, 0, 0);
1398 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 0, 1);
1399 //get the normals of from cell 2 and 4
1400 unitOuterNormal14 *= -1;
1401 unitOuterNormal12 *= -1;
1402 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol12, 0, 0);
1403 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 0, 1);
1404 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace12, 0, 0);
1405 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace41, 0, 1);
1406
1407 problem_.boundaryTypes(bcType, intersection12);
1408 PrimaryVariables boundValues(0.0);
1409
1410 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 0);
1411 if (bcType.isNeumann(pressEqIdx))
1412 {
1413 problem_.neumann(boundValues, intersection12);
1414 boundValues *= faceVol12;
1415 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 0);
1416 }
1417 if (bcType.hasDirichlet())
1418 {
1419 problem_.dirichlet(boundValues, intersection12);
1420 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 0);
1421 }
1422
1423 // 'intersection14' is on boundary
1424 if (intersection14.boundary())
1425 {
1426 problem_.boundaryTypes(bcType, intersection14);
1427 boundValues = 0.0;
1428
1429 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 3);
1430 if (bcType.isNeumann(pressEqIdx))
1431 {
1432 problem_.neumann(boundValues, intersection14);
1433 boundValues *= faceVol41;
1434 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 3);
1435 }
1436 if (bcType.hasDirichlet())
1437 {
1438 problem_.dirichlet(boundValues, intersection14);
1439 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 3);
1440 }
1441
1442 interactionVolumes_[globalVertIdx1234].setOutsideFace(1);
1443 interactionVolumes_[globalVertIdx1234].setOutsideFace(2);
1444 }
1445
1446 // 'intersection14' is inside
1447 else if (intersection14.neighbor())
1448 {
1449 // neighbor cell 3
1450 // access neighbor cell 3
1451 auto element4 = intersection14.outside();
1452 int eIdxGlobal4 = problem_.variables().index(element4);
1453
1454 if ((element.level() == element4.level() && eIdxGlobal1 > eIdxGlobal4)
1455 || element.level() < element4.level())
1456 {
1457 interactionVolumes_[globalVertIdx1234].reset();
1458 continue;
1459 }
1460
1461 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection14.indexInOutside(), 3, 0);
1462
1463 //store pointer 4
1464 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element4, 3);
1465
1466 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 3, 0);
1467 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 3, 0);
1468 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace41, 3, 0);
1469
1470 bool finished = false;
1471
1472 // get the information of the face 'isIt34' between cell3 and cell4 (locally numbered)
1473 for (const auto& intersection4
1474 : intersections(problem_.gridView(), element4))
1475 {
1476 if (intersection4.boundary())
1477 {
1478 for (int i = 0; i < intersection4.geometry().corners(); ++i)
1479 {
1480 int localVertIdx4corner = refElement.subEntity(intersection4.indexInInside(), dim - 1, i,
1481 dim);
1482
1483 int globalVertIdx4corner = problem_.variables().index(element4.template subEntity<dim>(localVertIdx4corner));
1484
1485 if (globalVertIdx4corner == globalVertIdx1234)
1486 {
1487 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInInside(), 3,
1488 1);
1489
1490 const GlobalPosition& globalPosFace34 = intersection4.geometry().center();
1491 Scalar faceVol34 = intersection4.geometry().volume() / 2.0;
1492
1493 DimVector unitOuterNormal43 = intersection4.centerUnitOuterNormal();
1494
1495 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 3, 1);
1496 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 3, 1);
1497 interactionVolumes_[globalVertIdx1234].setFacePosition(globalPosFace34, 3, 1);
1498
1499 problem_.boundaryTypes(bcType, intersection4);
1500 boundValues = 0.0;
1501
1502 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 2);
1503 if (bcType.isNeumann(pressEqIdx))
1504 {
1505 problem_.neumann(boundValues, intersection4);
1506 boundValues *= faceVol34;
1507 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 2);
1508 }
1509 if (bcType.hasDirichlet())
1510 {
1511 problem_.dirichlet(boundValues, intersection4);
1512 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 2);
1513 }
1514
1515 interactionVolumes_[globalVertIdx1234].setOutsideFace(1);
1516
1517 if (element.level() == element4.level())
1518 {
1519 innerBoundaryVolumeFaces_[eIdxGlobal1][intersection14.indexInInside()] = true;
1520 innerBoundaryVolumeFaces_[eIdxGlobal4][intersection14.indexInOutside()] = true;
1521 }
1522 if (element.level() < element4.level())
1523 {
1524 innerBoundaryVolumeFaces_[eIdxGlobal4][intersection14.indexInOutside()] = true;
1525 }
1526 if (element.level() > element4.level())
1527 {
1528 innerBoundaryVolumeFaces_[eIdxGlobal1][intersection14.indexInInside()] = true;
1529 }
1530
1531
1532 finished = true;
1533
1534 break;
1535 }
1536 }
1537 }
1538 if (finished)
1539 {
1540 break;
1541 }
1542 }
1543 if (!finished)
1544 {
1545 DUNE_THROW(
1546 Dune::NotImplemented,
1547 "fvmpfao2pfaboundpressure2p_adaptive.hh: boundary shape not available as interaction volume shape");
1548 }
1549 }
1550 else
1551 {
1552 DUNE_THROW(Dune::NotImplemented,
1553 "fvmpfao2pfaboundpressure2p_adaptive.hh: interface type not supported!");
1554 }
1555 }
1556
1557 } // end all intersections
1558 } // end grid traversal
1559
1560 return;
1561}
1562
1563// TODO doc me!
1564template<class TypeTag>
1565void FvMpfaL2dPressure2pAdaptive<TypeTag>::printInteractionVolumes()
1566{
1567 for (const auto& vertex : vertices(problem_.gridView()))
1568 {
1569 int vIdxGlobal = problem_.variables().index(vertex);
1570
1571 InteractionVolume& interactionVolume = interactionVolumes_[vIdxGlobal];
1572
1573 if (interactionVolume.getElementNumber() > 2)
1574 {
1575 interactionVolume.printInteractionVolumeInfo();
1576 std::cout << "global vertex index: " << vIdxGlobal << "\n";
1577 if (interactionVolume.getElementNumber() == 3)
1578 {
1579 auto element1 = interactionVolume.getSubVolumeElement(0);
1580 auto element2 = interactionVolume.getSubVolumeElement(1);
1581 auto element4 = interactionVolume.getSubVolumeElement(3);
1582
1583 int eIdxGlobal1 = problem_.variables().index(element1);
1584 int eIdxGlobal2 = problem_.variables().index(element2);
1585 int eIdxGlobal4 = problem_.variables().index(element4);
1586
1587 std::cout << "global element index 1: " << eIdxGlobal1 << "\n";
1588 std::cout << "global element index 2: " << eIdxGlobal2 << "\n";
1589 std::cout << "global element index 4: " << eIdxGlobal4 << "\n";
1590 }
1591 if (interactionVolume.getElementNumber() == 4)
1592 {
1593 auto element1 = interactionVolume.getSubVolumeElement(0);
1594 auto element2 = interactionVolume.getSubVolumeElement(1);
1595 auto element3 = interactionVolume.getSubVolumeElement(2);
1596 auto element4 = interactionVolume.getSubVolumeElement(3);
1597
1598 int eIdxGlobal1 = problem_.variables().index(element1);
1599 int eIdxGlobal2 = problem_.variables().index(element2);
1600 int eIdxGlobal3 = problem_.variables().index(element3);
1601 int eIdxGlobal4 = problem_.variables().index(element4);
1602
1603 std::cout << "global element index 1: " << eIdxGlobal1 << "\n";
1604 std::cout << "global element index 2: " << eIdxGlobal2 << "\n";
1605 std::cout << "global element index 3: " << eIdxGlobal3 << "\n";
1606 std::cout << "global element index 4: " << eIdxGlobal4 << "\n";
1607 }
1608 }
1609 }
1610}
1611
1612// only for 2-D general quadrilateral
1613// TODO doc me!
1614template<class TypeTag>
1615void FvMpfaL2dPressure2pAdaptive<TypeTag>::assemble()
1616{
1617 // initialization: set global matrix this->A_ to zero
1618 this->A_ = 0;
1619 this->f_ = 0;
1620
1621 // run through all vertices
1622 for (const auto& vertex : vertices(problem_.gridView()))
1623 {
1624 int vIdxGlobal = problem_.variables().index(vertex);
1625
1626 InteractionVolume& interactionVolume = interactionVolumes_[vIdxGlobal];
1627
1628 if (interactionVolume.isInnerVolume())
1629 {
1630 if (interactionVolume.getElementNumber() == 4)
1631 {
1632 auto element1 = interactionVolume.getSubVolumeElement(0);
1633 auto element2 = interactionVolume.getSubVolumeElement(1);
1634 auto element3 = interactionVolume.getSubVolumeElement(2);
1635 auto element4 = interactionVolume.getSubVolumeElement(3);
1636
1637 // get global coordinate of cell centers
1638 const GlobalPosition& globalPos1 = element1.geometry().center();
1639 const GlobalPosition& globalPos2 = element2.geometry().center();
1640 const GlobalPosition& globalPos3 = element3.geometry().center();
1641 const GlobalPosition& globalPos4 = element4.geometry().center();
1642
1643 // cell volumes
1644 Scalar volume1 = element1.geometry().volume();
1645 Scalar volume2 = element2.geometry().volume();
1646 Scalar volume3 = element3.geometry().volume();
1647 Scalar volume4 = element4.geometry().volume();
1648
1649 // cell index
1650 int eIdxGlobal1 = problem_.variables().index(element1);
1651 int eIdxGlobal2 = problem_.variables().index(element2);
1652 int eIdxGlobal3 = problem_.variables().index(element3);
1653 int eIdxGlobal4 = problem_.variables().index(element4);
1654
1655 //get the cell Data
1656 CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
1657 CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
1658 CellData& cellData3 = problem_.variables().cellData(eIdxGlobal3);
1659 CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
1660
1661 // evaluate right hand side
1662 PrimaryVariables source(0.0);
1663 problem_.source(source, element1);
1664 this->f_[eIdxGlobal1] += volume1 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1665 problem_.source(source, element2);
1666 this->f_[eIdxGlobal2] += volume2 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1667 problem_.source(source, element3);
1668 this->f_[eIdxGlobal3] += volume3 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1669 problem_.source(source, element4);
1670 this->f_[eIdxGlobal4] += volume4 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1671
1672 this->f_[eIdxGlobal1] += evaluateErrorTerm_(cellData1) * volume1 / (4.0);
1673 this->f_[eIdxGlobal2] += evaluateErrorTerm_(cellData2) * volume2 / (4.0);
1674 this->f_[eIdxGlobal3] += evaluateErrorTerm_(cellData3) * volume3 / (4.0);
1675 this->f_[eIdxGlobal4] += evaluateErrorTerm_(cellData4) * volume4 / (4.0);
1676
1677 //get mobilities of the phases
1678 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
1679 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
1680
1681 //compute total mobility of cell 1
1682 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
1683
1684 //get mobilities of the phases
1685 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
1686 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
1687
1688 //compute total mobility of cell 1
1689 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
1690
1691 //get mobilities of the phases
1692 Dune::FieldVector<Scalar, numPhases> lambda3(cellData3.mobility(wPhaseIdx));
1693 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
1694
1695 //compute total mobility of cell 1
1696 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
1697
1698 //get mobilities of the phases
1699 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
1700 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
1701
1702 //compute total mobility of cell 1
1703 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
1704
1705 std::vector < DimVector > lambda(2 * dim);
1706
1707 lambda[0][0] = lambdaTotal1;
1708 lambda[0][1] = lambdaTotal1;
1709 lambda[1][0] = lambdaTotal2;
1710 lambda[1][1] = lambdaTotal2;
1711 lambda[2][0] = lambdaTotal3;
1712 lambda[2][1] = lambdaTotal3;
1713 lambda[3][0] = lambdaTotal4;
1714 lambda[3][1] = lambdaTotal4;
1715
1716 //add capillary pressure and gravity terms to right-hand-side
1717 //calculate capillary pressure velocity
1718 Dune::FieldVector<Scalar, 2 * dim> pc(0);
1719 pc[0] = cellData1.capillaryPressure();
1720 pc[1] = cellData2.capillaryPressure();
1721 pc[2] = cellData3.capillaryPressure();
1722 pc[3] = cellData4.capillaryPressure();
1723
1724 Dune::FieldVector<Scalar, 2 * dim> gravityDiff(0);
1725
1726 gravityDiff[0] = (problem_.bBoxMax() - globalPos1) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1727 gravityDiff[1] = (problem_.bBoxMax() - globalPos2) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1728 gravityDiff[2] = (problem_.bBoxMax() - globalPos3) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1729 gravityDiff[3] = (problem_.bBoxMax() - globalPos4) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1730
1731 pc += gravityDiff;
1732
1733 Dune::FieldVector<Scalar, 2 * dim> pcFlux(0);
1734
1735 Scalar pcPotential12 = 0;
1736 Scalar pcPotential14 = 0;
1737 Scalar pcPotential32 = 0;
1738 Scalar pcPotential34 = 0;
1739
1740 DimVector Tu(0);
1741 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
1742 Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> T(0);
1743
1744 int transmissibilityType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 0, 1, 2, 3);
1745
1746 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
1747 {
1748 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
1749 {
1750 T *= 2;
1751 }
1752 this->A_[eIdxGlobal1][eIdxGlobal2] += T[1][0];
1753 this->A_[eIdxGlobal1][eIdxGlobal3] += T[1][1];
1754 this->A_[eIdxGlobal1][eIdxGlobal1] += T[1][2];
1755
1756 this->A_[eIdxGlobal2][eIdxGlobal2] -= T[1][0];
1757 this->A_[eIdxGlobal2][eIdxGlobal3] -= T[1][1];
1758 this->A_[eIdxGlobal2][eIdxGlobal1] -= T[1][2];
1759
1760 u[0] = pc[1];
1761 u[1] = pc[2];
1762 u[2] = pc[0];
1763
1764 T.mv(u, Tu);
1765
1766 pcFlux[0] = Tu[1];
1767 pcPotential12 = Tu[1];
1768 }
1769 else if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
1770 {
1771 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
1772 {
1773 T *= 2;
1774 }
1775 this->A_[eIdxGlobal1][eIdxGlobal1] += T[1][0];
1776 this->A_[eIdxGlobal1][eIdxGlobal4] += T[1][1];
1777 this->A_[eIdxGlobal1][eIdxGlobal2] += T[1][2];
1778
1779 this->A_[eIdxGlobal2][eIdxGlobal1] -= T[1][0];
1780 this->A_[eIdxGlobal2][eIdxGlobal4] -= T[1][1];
1781 this->A_[eIdxGlobal2][eIdxGlobal2] -= T[1][2];
1782
1783 u[0] = pc[0];
1784 u[1] = pc[3];
1785 u[2] = pc[1];
1786
1787 T.mv(u, Tu);
1788
1789 pcFlux[0] = Tu[1];
1790 pcPotential12 = Tu[1];
1791 }
1792
1793 transmissibilityType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 1, 2, 3, 0);
1794
1795 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
1796 {
1797 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
1798 {
1799 T *= 2;
1800 }
1801 this->A_[eIdxGlobal2][eIdxGlobal3] += T[1][0];
1802 this->A_[eIdxGlobal2][eIdxGlobal4] += T[1][1];
1803 this->A_[eIdxGlobal2][eIdxGlobal2] += T[1][2];
1804
1805 this->A_[eIdxGlobal3][eIdxGlobal3] -= T[1][0];
1806 this->A_[eIdxGlobal3][eIdxGlobal4] -= T[1][1];
1807 this->A_[eIdxGlobal3][eIdxGlobal2] -= T[1][2];
1808
1809 u[0] = pc[2];
1810 u[1] = pc[3];
1811 u[2] = pc[1];
1812
1813 T.mv(u, Tu);
1814
1815 pcFlux[1] = Tu[1];
1816 pcPotential32 = -Tu[1];
1817 }
1818 else if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
1819 {
1820 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
1821 {
1822 T *= 2;
1823 }
1824 this->A_[eIdxGlobal2][eIdxGlobal2] += T[1][0];
1825 this->A_[eIdxGlobal2][eIdxGlobal1] += T[1][1];
1826 this->A_[eIdxGlobal2][eIdxGlobal3] += T[1][2];
1827
1828 this->A_[eIdxGlobal3][eIdxGlobal2] -= T[1][0];
1829 this->A_[eIdxGlobal3][eIdxGlobal1] -= T[1][1];
1830 this->A_[eIdxGlobal3][eIdxGlobal3] -= T[1][2];
1831
1832 u[0] = pc[1];
1833 u[1] = pc[0];
1834 u[2] = pc[2];
1835
1836 T.mv(u, Tu);
1837
1838 pcFlux[1] = Tu[1];
1839 pcPotential32 = -Tu[1];
1840 }
1841
1842 transmissibilityType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 2, 3, 0, 1);
1843
1844 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
1845 {
1846 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
1847 {
1848 T *= 2;
1849 }
1850 this->A_[eIdxGlobal3][eIdxGlobal4] += T[1][0];
1851 this->A_[eIdxGlobal3][eIdxGlobal1] += T[1][1];
1852 this->A_[eIdxGlobal3][eIdxGlobal3] += T[1][2];
1853
1854 this->A_[eIdxGlobal4][eIdxGlobal4] -= T[1][0];
1855 this->A_[eIdxGlobal4][eIdxGlobal1] -= T[1][1];
1856 this->A_[eIdxGlobal4][eIdxGlobal3] -= T[1][2];
1857
1858 u[0] = pc[3];
1859 u[1] = pc[0];
1860 u[2] = pc[2];
1861
1862 T.mv(u, Tu);
1863
1864 pcFlux[2] = Tu[1];
1865 pcPotential34 = Tu[1];
1866 }
1867 else if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
1868 {
1869 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
1870 {
1871 T *= 2;
1872 }
1873 this->A_[eIdxGlobal3][eIdxGlobal3] += T[1][0];
1874 this->A_[eIdxGlobal3][eIdxGlobal2] += T[1][1];
1875 this->A_[eIdxGlobal3][eIdxGlobal4] += T[1][2];
1876
1877 this->A_[eIdxGlobal4][eIdxGlobal3] -= T[1][0];
1878 this->A_[eIdxGlobal4][eIdxGlobal2] -= T[1][1];
1879 this->A_[eIdxGlobal4][eIdxGlobal4] -= T[1][2];
1880
1881 u[0] = pc[2];
1882 u[1] = pc[1];
1883 u[2] = pc[3];
1884
1885 T.mv(u, Tu);
1886
1887 pcFlux[2] = Tu[1];
1888 pcPotential34 = Tu[1];
1889 }
1890
1891 transmissibilityType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 3, 0, 1, 2);
1892
1893 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
1894 {
1895 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
1896 {
1897 T *= 2;
1898 }
1899 this->A_[eIdxGlobal4][eIdxGlobal1] += T[1][0];
1900 this->A_[eIdxGlobal4][eIdxGlobal2] += T[1][1];
1901 this->A_[eIdxGlobal4][eIdxGlobal4] += T[1][2];
1902
1903 this->A_[eIdxGlobal1][eIdxGlobal1] -= T[1][0];
1904 this->A_[eIdxGlobal1][eIdxGlobal2] -= T[1][1];
1905 this->A_[eIdxGlobal1][eIdxGlobal4] -= T[1][2];
1906
1907 u[0] = pc[0];
1908 u[1] = pc[1];
1909 u[2] = pc[3];
1910
1911 T.mv(u, Tu);
1912
1913 pcFlux[3] = Tu[1];
1914 pcPotential14 = -Tu[1];
1915 }
1916 else if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
1917 {
1918 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
1919 {
1920 T *= 2;
1921 }
1922 this->A_[eIdxGlobal4][eIdxGlobal4] += T[1][0];
1923 this->A_[eIdxGlobal4][eIdxGlobal3] += T[1][1];
1924 this->A_[eIdxGlobal4][eIdxGlobal1] += T[1][2];
1925
1926 this->A_[eIdxGlobal1][eIdxGlobal4] -= T[1][0];
1927 this->A_[eIdxGlobal1][eIdxGlobal3] -= T[1][1];
1928 this->A_[eIdxGlobal1][eIdxGlobal1] -= T[1][2];
1929
1930 u[0] = pc[3];
1931 u[1] = pc[2];
1932 u[2] = pc[0];
1933
1934 T.mv(u, Tu);
1935
1936 pcFlux[3] = Tu[1];
1937 pcPotential14 = -Tu[1];
1938 }
1939
1940 if (pc[0] == 0 && pc[1] == 0 && pc[2] == 0 && pc[3] == 0)
1941 {
1942 continue;
1943 }
1944
1945 //compute mobilities of face 1
1946 Dune::FieldVector<Scalar, numPhases> lambda12Upw(0.0);
1947 lambda12Upw[wPhaseIdx] = (pcPotential12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
1948 lambda12Upw[nPhaseIdx] = (pcPotential12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
1949
1950 //compute mobilities of face 4
1951 Dune::FieldVector<Scalar, numPhases> lambda14Upw(0.0);
1952 lambda14Upw[wPhaseIdx] = (pcPotential14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
1953 lambda14Upw[nPhaseIdx] = (pcPotential14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
1954
1955 //compute mobilities of face 2
1956 Dune::FieldVector<Scalar, numPhases> lambda32Upw(0.0);
1957 lambda32Upw[wPhaseIdx] = (pcPotential32 >= 0) ? lambda3[wPhaseIdx] : lambda2[wPhaseIdx];
1958 lambda32Upw[nPhaseIdx] = (pcPotential32 >= 0) ? lambda3[nPhaseIdx] : lambda2[nPhaseIdx];
1959
1960 //compute mobilities of face 3
1961 Dune::FieldVector<Scalar, numPhases> lambda34Upw(0.0);
1962 lambda34Upw[wPhaseIdx] = (pcPotential34 >= 0) ? lambda3[wPhaseIdx] : lambda4[wPhaseIdx];
1963 lambda34Upw[nPhaseIdx] = (pcPotential34 >= 0) ? lambda3[nPhaseIdx] : lambda4[nPhaseIdx];
1964
1965 for (int i = 0; i < numPhases; i++)
1966 {
1967 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
1968 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
1969 Scalar lambdaT32 = lambda32Upw[wPhaseIdx] + lambda32Upw[nPhaseIdx];
1970 Scalar lambdaT34 = lambda34Upw[wPhaseIdx] + lambda34Upw[nPhaseIdx];
1971 Scalar fracFlow12 = (lambdaT12 > threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
1972 Scalar fracFlow14 = (lambdaT14 > threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
1973 Scalar fracFlow32 = (lambdaT32 > threshold_) ? lambda32Upw[i] / (lambdaT32) : 0.0;
1974 Scalar fracFlow34 = (lambdaT34 > threshold_) ? lambda34Upw[i] / (lambdaT34) : 0.0;
1975
1976 Dune::FieldVector<Scalar, 2 * dim> pcFluxReal(pcFlux);
1977
1978 pcFluxReal[0] *= fracFlow12;
1979 pcFluxReal[1] *= fracFlow32;
1980 pcFluxReal[2] *= fracFlow34;
1981 pcFluxReal[3] *= fracFlow14;
1982
1983 switch (pressureType_)
1984 {
1985 case pw:
1986 {
1987 if (i == nPhaseIdx)
1988 {
1989 //add capillary pressure term to right hand side
1990 this->f_[eIdxGlobal1] -= (pcFluxReal[0] - pcFluxReal[3]);
1991 this->f_[eIdxGlobal2] -= (pcFluxReal[1] - pcFluxReal[0]);
1992 this->f_[eIdxGlobal3] -= (pcFluxReal[2] - pcFluxReal[1]);
1993 this->f_[eIdxGlobal4] -= (pcFluxReal[3] - pcFluxReal[2]);
1994
1995 }
1996 break;
1997 }
1998 case pn:
1999 {
2000 if (i == wPhaseIdx)
2001 {
2002 //add capillary pressure term to right hand side
2003 this->f_[eIdxGlobal1] += (pcFluxReal[0] - pcFluxReal[3]);
2004 this->f_[eIdxGlobal2] += (pcFluxReal[1] - pcFluxReal[0]);
2005 this->f_[eIdxGlobal3] += (pcFluxReal[2] - pcFluxReal[1]);
2006 this->f_[eIdxGlobal4] += (pcFluxReal[3] - pcFluxReal[2]);
2007 }
2008 break;
2009 }
2010 }
2011 }
2012 }
2013 else if (interactionVolume.getElementNumber() == 3)
2014 {
2015 auto element1 = interactionVolume.getSubVolumeElement(0);
2016 auto element2 = interactionVolume.getSubVolumeElement(1);
2017 auto element4 = interactionVolume.getSubVolumeElement(3);
2018
2019 // get global coordinate of cell centers
2020 const GlobalPosition& globalPos1 = element1.geometry().center();
2021 const GlobalPosition& globalPos2 = element2.geometry().center();
2022 const GlobalPosition& globalPos4 = element4.geometry().center();
2023
2024 // cell volumes
2025 Scalar volume1 = element1.geometry().volume();
2026 Scalar volume2 = element2.geometry().volume();
2027
2028 // cell index
2029 int eIdxGlobal1 = problem_.variables().index(element1);
2030 int eIdxGlobal2 = problem_.variables().index(element2);
2031 int eIdxGlobal4 = problem_.variables().index(element4);
2032
2033 //get the cell Data
2034 CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
2035 CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
2036 CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
2037
2038 // evaluate right hand side -> only add source for the cells without hanging node!
2039 // In doing so every cell gets the source from 4 vertices and the division by 4 is correct!
2040 PrimaryVariables source(0.0);
2041 problem_.source(source, element1);
2042 this->f_[eIdxGlobal1] += volume1 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
2043 problem_.source(source, element2);
2044 this->f_[eIdxGlobal2] += volume2 / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
2045
2046 this->f_[eIdxGlobal1] += evaluateErrorTerm_(cellData1) * volume1 / (4.0);
2047 this->f_[eIdxGlobal2] += evaluateErrorTerm_(cellData2) * volume2 / (4.0);
2048
2049 //get mobilities of the phases
2050 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
2051 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
2052
2053 //compute total mobility of cell 1
2054 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
2055
2056 //get mobilities of the phases
2057 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
2058 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
2059
2060 //compute total mobility of cell 1
2061 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
2062
2063 //get mobilities of the phases
2064 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
2065 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
2066
2067 //compute total mobility of cell 1
2068 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
2069
2070 std::vector < DimVector > lambda(4);
2071
2072 lambda[0][0] = lambdaTotal1;
2073 lambda[0][1] = lambdaTotal1;
2074 lambda[1][0] = lambdaTotal2;
2075 lambda[1][1] = lambdaTotal2;
2076 lambda[3][0] = lambdaTotal4;
2077 lambda[3][1] = lambdaTotal4;
2078
2079 //add capillary pressure and gravity terms to right-hand-side
2080 //calculate capillary pressure velocity
2081 Dune::FieldVector<Scalar, 3> pc(0);
2082 pc[0] = cellData1.capillaryPressure();
2083 pc[1] = cellData2.capillaryPressure();
2084 pc[2] = cellData4.capillaryPressure();
2085
2086 Dune::FieldVector<Scalar, 3> gravityDiff(0);
2087
2088 gravityDiff[0] = (problem_.bBoxMax() - globalPos1) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
2089 gravityDiff[1] = (problem_.bBoxMax() - globalPos2) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
2090 gravityDiff[2] = (problem_.bBoxMax() - globalPos4) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
2091
2092 pc += gravityDiff;
2093
2094 Dune::FieldVector<Scalar, 3> pcFlux(0);
2095
2096 Scalar pcPotential12 = 0;
2097 Scalar pcPotential14 = 0;
2098 Scalar pcPotential24 = 0;
2099
2100 DimVector Tu(0);
2101 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
2102 Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> T(0);
2103
2104 int transmissibilityType = transmissibilityCalculator_.calculateTransmissibility(T, interactionVolume, lambda, 0, 1, 3, 3);
2105
2106 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
2107 {
2108 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
2109 {
2110 T *= 2;
2111 }
2112 this->A_[eIdxGlobal1][eIdxGlobal2] += T[1][0];
2113 this->A_[eIdxGlobal1][eIdxGlobal4] += T[1][1];
2114 this->A_[eIdxGlobal1][eIdxGlobal1] += T[1][2];
2115
2116 this->A_[eIdxGlobal2][eIdxGlobal2] -= T[1][0];
2117 this->A_[eIdxGlobal2][eIdxGlobal4] -= T[1][1];
2118 this->A_[eIdxGlobal2][eIdxGlobal1] -= T[1][2];
2119
2120 u[0] = pc[1];
2121 u[1] = pc[2];
2122 u[2] = pc[0];
2123
2124 T.mv(u, Tu);
2125
2126 pcFlux[0] = Tu[1];
2127 pcPotential12 = Tu[1];
2128 }
2129 else if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
2130 {
2131 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
2132 {
2133 T *= 2;
2134 }
2135 this->A_[eIdxGlobal1][eIdxGlobal1] += T[1][0];
2136 this->A_[eIdxGlobal1][eIdxGlobal4] += T[1][1];
2137 this->A_[eIdxGlobal1][eIdxGlobal2] += T[1][2];
2138
2139 this->A_[eIdxGlobal2][eIdxGlobal1] -= T[1][0];
2140 this->A_[eIdxGlobal2][eIdxGlobal4] -= T[1][1];
2141 this->A_[eIdxGlobal2][eIdxGlobal2] -= T[1][2];
2142
2143 u[0] = pc[0];
2144 u[1] = pc[2];
2145 u[2] = pc[1];
2146
2147 T.mv(u, Tu);
2148
2149 pcFlux[0] = Tu[1];
2150 pcPotential12 = Tu[1];
2151 }
2152 else
2153 {
2154 DUNE_THROW(Dune::RangeError, "Could not calculate Tranmissibility!");
2155 }
2156
2157 transmissibilityType = transmissibilityCalculator_.calculateLeftHNTransmissibility(T, interactionVolume, lambda, 1, 3, 0);
2158
2159 if (transmissibilityType == TransmissibilityCalculator::leftTriangle)
2160 {
2161 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
2162 {
2163 T *= 2;
2164 }
2165 this->A_[eIdxGlobal2][eIdxGlobal2] += T[1][0];
2166 this->A_[eIdxGlobal2][eIdxGlobal1] += T[1][1];
2167 this->A_[eIdxGlobal2][eIdxGlobal4] += T[1][2];
2168
2169 this->A_[eIdxGlobal4][eIdxGlobal2] -= T[1][0];
2170 this->A_[eIdxGlobal4][eIdxGlobal1] -= T[1][1];
2171 this->A_[eIdxGlobal4][eIdxGlobal4] -= T[1][2];
2172
2173 u[0] = pc[1];
2174 u[1] = pc[0];
2175 u[2] = pc[2];
2176
2177 T.mv(u, Tu);
2178
2179 pcFlux[1] = Tu[1];
2180 pcPotential24 = Tu[1];
2181 }
2182 else
2183 {
2184 DUNE_THROW(Dune::RangeError, "Could not calculate Tranmissibility!");
2185 }
2186
2187 transmissibilityType = transmissibilityCalculator_.calculateRightHNTransmissibility(T, interactionVolume, lambda, 3, 0, 1);
2188
2189 if (transmissibilityType == TransmissibilityCalculator::rightTriangle)
2190 {
2191 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
2192 {
2193 T *= 2;
2194 }
2195
2196 this->A_[eIdxGlobal4][eIdxGlobal1] += T[1][0];
2197 this->A_[eIdxGlobal4][eIdxGlobal2] += T[1][1];
2198 this->A_[eIdxGlobal4][eIdxGlobal4] += T[1][2];
2199
2200 this->A_[eIdxGlobal1][eIdxGlobal1] -= T[1][0];
2201 this->A_[eIdxGlobal1][eIdxGlobal2] -= T[1][1];
2202 this->A_[eIdxGlobal1][eIdxGlobal4] -= T[1][2];
2203
2204 u[0] = pc[0];
2205 u[1] = pc[1];
2206 u[2] = pc[2];
2207
2208 T.mv(u, Tu);
2209
2210 pcFlux[2] = Tu[1];
2211 pcPotential14 = -Tu[1];
2212 }
2213 else
2214 {
2215 DUNE_THROW(Dune::RangeError, "Could not calculate Tranmissibility!");
2216 }
2217
2218 if (pc[0] == 0 && pc[1] == 0 && pc[2] == 0)
2219 {
2220 continue;
2221 }
2222
2223 //compute mobilities of face 1
2224 Dune::FieldVector<Scalar, numPhases> lambda12Upw(0.0);
2225 lambda12Upw[wPhaseIdx] = (pcPotential12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
2226 lambda12Upw[nPhaseIdx] = (pcPotential12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
2227
2228 //compute mobilities of face 4
2229 Dune::FieldVector<Scalar, numPhases> lambda14Upw(0.0);
2230 lambda14Upw[wPhaseIdx] = (pcPotential14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
2231 lambda14Upw[nPhaseIdx] = (pcPotential14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
2232
2233 //compute mobilities of face 2
2234 Dune::FieldVector<Scalar, numPhases> lambda24Upw(0.0);
2235 lambda24Upw[wPhaseIdx] = (pcPotential24 >= 0) ? lambda2[wPhaseIdx] : lambda4[wPhaseIdx];
2236 lambda24Upw[nPhaseIdx] = (pcPotential24 >= 0) ? lambda2[nPhaseIdx] : lambda4[nPhaseIdx];
2237
2238 for (int i = 0; i < numPhases; i++)
2239 {
2240 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
2241 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
2242 Scalar lambdaT24 = lambda24Upw[wPhaseIdx] + lambda24Upw[nPhaseIdx];
2243 Scalar fracFlow12 = (lambdaT12 > threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
2244 Scalar fracFlow14 = (lambdaT14 > threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
2245 Scalar fracFlow24 = (lambdaT24 > threshold_) ? lambda24Upw[i] / (lambdaT24) : 0.0;
2246
2247 Dune::FieldVector<Scalar, 3> pcFluxReal(pcFlux);
2248
2249 pcFluxReal[0] *= fracFlow12;
2250 pcFluxReal[1] *= fracFlow24;
2251 pcFluxReal[2] *= fracFlow14;
2252
2253 switch (pressureType_)
2254 {
2255 case pw:
2256 {
2257 if (i == nPhaseIdx)
2258 {
2259 //add capillary pressure term to right hand side
2260 this->f_[eIdxGlobal1] -= (pcFluxReal[0] - pcFluxReal[2]);
2261 this->f_[eIdxGlobal2] -= (pcFluxReal[1] - pcFluxReal[0]);
2262 this->f_[eIdxGlobal4] -= (pcFluxReal[2] - pcFluxReal[1]);
2263
2264 }
2265 break;
2266 }
2267 case pn:
2268 {
2269 if (i == wPhaseIdx)
2270 {
2271 //add capillary pressure term to right hand side
2272 this->f_[eIdxGlobal1] += (pcFluxReal[0] - pcFluxReal[2]);
2273 this->f_[eIdxGlobal2] += (pcFluxReal[1] - pcFluxReal[0]);
2274 this->f_[eIdxGlobal4] += (pcFluxReal[2] - pcFluxReal[1]);
2275 }
2276 break;
2277 }
2278 }
2279 }
2280 }
2281 else
2282 {
2283 DUNE_THROW(Dune::NotImplemented, "Unknown interactionvolume type!");
2284 }
2285 }
2286
2287 // at least one face on boundary!
2288 else
2289 {
2290 for (int elemIdx = 0; elemIdx < 2 * dim; elemIdx++)
2291 {
2292 bool isOutside = false;
2293 for (int fIdx = 0; fIdx < dim; fIdx++)
2294 {
2295 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
2296 if (interactionVolume.isOutsideFace(intVolFaceIdx))
2297 {
2298 isOutside = true;
2299 break;
2300 }
2301 }
2302 if (isOutside)
2303 {
2304 continue;
2305 }
2306
2307 auto element = interactionVolume.getSubVolumeElement(elemIdx);
2308
2309 // get global coordinate of cell centers
2310 const GlobalPosition& globalPos = element.geometry().center();
2311
2312 // cell volumes
2313 Scalar volume = element.geometry().volume();
2314
2315 // cell index
2316 int eIdxGlobal = problem_.variables().index(element);
2317
2318 //get the cell Data
2319 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
2320
2321 //permeability vector at boundary
2322 DimMatrix permeability(problem_.spatialParams().intrinsicPermeability(element));
2323
2324 // evaluate right hand side
2325 PrimaryVariables source(0);
2326 problem_.source(source, element);
2327 this->f_[eIdxGlobal] += volume / (4.0) * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
2328
2329 this->f_[eIdxGlobal] += evaluateErrorTerm_(cellData) * volume / (4.0);
2330
2331 //get mobilities of the phases
2332 Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
2333 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
2334
2335 Scalar pc = cellData.capillaryPressure();
2336
2337 Scalar gravityDiff = (problem_.bBoxMax() - globalPos) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
2338
2339 pc += gravityDiff; //minus because of gravity definition!
2340
2341 for (int fIdx = 0; fIdx < dim; fIdx++)
2342 {
2343 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
2344
2345 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
2346 {
2347
2348 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressEqIdx))
2349 {
2350 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
2351
2352 const auto refElement = referenceElement(element);
2353
2354 const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
2355
2356 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
2357
2358 DimVector distVec(globalPosFace - globalPos);
2359 Scalar dist = distVec.two_norm();
2360 DimVector unitDistVec(distVec);
2361 unitDistVec /= dist;
2362
2363 Scalar faceArea = interactionVolume.getFaceArea(elemIdx, fIdx);
2364
2365 // get pc and lambda at the boundary
2366 Scalar satWBound = cellData.saturation(wPhaseIdx);
2367 //check boundary sat at face 1
2368 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
2369 {
2370 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
2371 switch (saturationType_)
2372 {
2373 case sw:
2374 {
2375 satWBound = satBound;
2376 break;
2377 }
2378 case sn:
2379 {
2380 satWBound = 1 - satBound;
2381 break;
2382 }
2383 }
2384
2385 }
2386
2387 // old material law interface is deprecated: Replace this by
2388 // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
2389 // after the release of 3.3, when the deprecated interface is no longer supported
2390 const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
2391
2392 Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
2393
2394 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
2395 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
2396
2397 pcBound += gravityDiffBound;
2398
2399 Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
2400 lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
2401 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
2402 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
2403
2404 Scalar potentialBound = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx];
2405 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
2406
2407 //calculate potential gradients
2408 Scalar potentialDiffW = 0;
2409 Scalar potentialDiffNw = 0;
2410 switch (pressureType_)
2411 {
2412 case pw:
2413 {
2414 potentialBound += density_[wPhaseIdx]*gdeltaZ;
2415 potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound) / dist;
2416 potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound - pcBound) / dist;
2417 break;
2418 }
2419 case pn:
2420 {
2421 potentialBound += density_[nPhaseIdx]*gdeltaZ;
2422 potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound + pcBound) / dist;
2423 potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound) / dist;
2424 break;
2425 }
2426 }
2427
2428 Scalar lambdaTotal = (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
2429 lambdaTotal += (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
2430
2431 DimVector permTimesNormal(0);
2432 permeability.mv(unitDistVec, permTimesNormal);
2433
2434 //calculate current matrix entry
2435 Scalar entry = lambdaTotal * (unitDistVec * permTimesNormal) / dist * faceArea;
2436
2437 //calculate right hand side
2438 Scalar pcFlux = 0;
2439
2440 switch (pressureType_)
2441 {
2442 case pw:
2443 {
2444 // calculate capillary pressure gradient
2445 DimVector pcGradient = unitDistVec;
2446 pcGradient *= (pc - pcBound) / dist;
2447
2448 //add capillary pressure term to right hand side
2449 pcFlux = 0.5 * (lambda[nPhaseIdx] + lambdaBound[nPhaseIdx])
2450 * (permTimesNormal * pcGradient) * faceArea;
2451
2452 break;
2453 }
2454 case pn:
2455 {
2456 // calculate capillary pressure gradient
2457 DimVector pcGradient = unitDistVec;
2458 pcGradient *= (pc - pcBound) / dist;
2459
2460 //add capillary pressure term to right hand side
2461 pcFlux = 0.5 * (lambda[wPhaseIdx] + lambdaBound[wPhaseIdx])
2462 * (permTimesNormal * pcGradient) * faceArea;
2463
2464 break;
2465
2466 }
2467 }
2468
2469 // set diagonal entry and right hand side entry
2470 this->A_[eIdxGlobal][eIdxGlobal] += entry;
2471 this->f_[eIdxGlobal] += entry * potentialBound;
2472
2473 if (pc == 0 && pcBound == 0)
2474 {
2475 continue;
2476 }
2477
2478 for (int i = 0; i < numPhases; i++)
2479 {
2480 switch (pressureType_)
2481 {
2482 case pw:
2483 {
2484 if (i == nPhaseIdx)
2485 {
2486 //add capillary pressure term to right hand side
2487 this->f_[eIdxGlobal] -= pcFlux;
2488 }
2489 break;
2490 }
2491 case pn:
2492 {
2493 if (i == wPhaseIdx)
2494 {
2495 //add capillary pressure term to right hand side
2496 this->f_[eIdxGlobal] += pcFlux;
2497 }
2498 break;
2499 }
2500 }
2501 }
2502
2503 }
2504 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx))
2505 {
2506 Scalar J = interactionVolume.getNeumannValues(intVolFaceIdx)[wPhaseIdx] / density_[wPhaseIdx];
2507 J += interactionVolume.getNeumannValues(intVolFaceIdx)[nPhaseIdx] / density_[nPhaseIdx];
2508
2509 this->f_[eIdxGlobal] -= J;
2510 }
2511 else
2512 {
2513 std::cout << "interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx)"
2514 << interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx) << "\n";
2515 DUNE_THROW(Dune::NotImplemented,
2516 "No valid boundary condition type defined for pressure equation!");
2517 }
2518 }
2519 }
2520 }
2521 } // end boundaries
2522
2523 } // end vertex iterator
2524
2525 // only do more if we have more than one process
2526 if (problem_.gridView().comm().size() > 1)
2527 {
2528 // set ghost and overlap element entries
2529 for (const auto& element : elements(problem_.gridView()))
2530 {
2531 if (element.partitionType() == Dune::InteriorEntity)
2532 continue;
2533
2534 // get the global index of the cell
2535 int eIdxGlobalI = problem_.variables().index(element);
2536
2537 this->A_[eIdxGlobalI] = 0.0;
2538 this->A_[eIdxGlobalI][eIdxGlobalI] = 1.0;
2539 this->f_[eIdxGlobalI] = this->pressure()[eIdxGlobalI];
2540 }
2541 }
2542
2543 return;
2544}
2545
2551template<class TypeTag>
2552void FvMpfaL2dPressure2pAdaptive<TypeTag>::updateMaterialLaws()
2553{
2554 // iterate through leaf grid an evaluate c0 at cell center
2555 for (const auto& element : elements(problem_.gridView()))
2556 {
2557 int eIdxGlobal = problem_.variables().index(element);
2558
2559 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
2560
2561 const Scalar satW = cellData.saturation(wPhaseIdx);
2562
2563 // old material law interface is deprecated: Replace this by
2564 // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
2565 // after the release of 3.3, when the deprecated interface is no longer supported
2566 const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
2567
2568 const Scalar pc = fluidMatrixInteraction.pc(satW);
2569
2570 cellData.setCapillaryPressure(pc);
2571
2572 // initialize mobilities
2573 const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
2574 const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
2575
2576 // initialize mobilities
2577 cellData.setMobility(wPhaseIdx, mobilityW);
2578 cellData.setMobility(nPhaseIdx, mobilityNw);
2579
2580 //initialize fractional flow functions
2581 cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
2582 cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
2583 }
2584 return;
2585}
2586
2587} // end namespace Dumux
2588#endif
deprecated.hh
Helpers for deprecation.
2dtransmissibilitycalculator.hh
Provides methods for transmissibility calculation 2-d.
linteractionvolume.hh
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Dumux
Definition: adapt.hh:29
Dumux::GetProp
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
Dumux::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Dumux::IOName::viscosity
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
Dumux::IOName::permeability
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
Dumux::IOName::pressure
std::string pressure(int phaseIdx) noexcept
I/O name of pressure for multiphase systems.
Definition: name.hh:34
Dumux::IOName::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Dumux::FvMpfaL2dPressure2pAdaptive
Grid adaptive finite volume MPFA L-method discretization of a two-phase flow pressure equation of the...
Definition: 2dpressureadaptive.hh:77
Dumux::FvMpfaL2dPressure2pAdaptive::storePressureSolution
void storePressureSolution(const Element &element)
Stores the pressure solution of a cell.
Definition: 2dpressureadaptive.hh:247
Dumux::FvMpfaL2dPressure2pAdaptive::addOutputVtkFields
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: 2dpressureadaptive.hh:376
Dumux::FvMpfaL2dPressure2pAdaptive::updateInteractionVolumeInfo
void updateInteractionVolumeInfo()
Updates interaction volumes.
Definition: 2dpressureadaptive.hh:185
Dumux::FvMpfaL2dPressure2pAdaptive::innerBoundaryVolumeFaces_
InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_
Vector marking faces which intersect the boundary.
Definition: 2dpressureadaptive.hh:492
Dumux::FvMpfaL2dPressure2pAdaptive::storePressureSolution
void storePressureSolution()
Globally stores the pressure solution.
Definition: 2dpressureadaptive.hh:233
Dumux::FvMpfaL2dPressure2pAdaptive::FvMpfaL2dPressure2pAdaptive
FvMpfaL2dPressure2pAdaptive(Problem &problem)
Constructs a FvMpfaL2dPressure2pAdaptive object.
Definition: 2dpressureadaptive.hh:446
Dumux::FvMpfaL2dPressure2pAdaptive::initialize
void initialize()
Initializes the pressure model.
Definition: 2dpressureadaptive.hh:202
Dumux::FvMpfaL2dPressure2pAdaptive::TransmissibilityCalculator
FvMpfaL2dTransmissibilityCalculator< TypeTag > TransmissibilityCalculator
Definition: 2dpressureadaptive.hh:155
Dumux::FvMpfaL2dPressure2pAdaptive::updateMaterialLaws
void updateMaterialLaws()
constitutive functions are initialized and stored in the variables object
Definition: 2dpressureadaptive.hh:2552
Dumux::FvMpfaL2dPressure2pAdaptive::update
void update()
Pressure update.
Definition: 2dpressureadaptive.hh:300
Dumux::FvMpfaL2dPressure2pAdaptive::interactionVolumes_
GlobalInteractionVolumeVector interactionVolumes_
Global Vector of interaction volumes.
Definition: 2dpressureadaptive.hh:491
Dumux::FvMpfaL2dTransmissibilityCalculator
Provides methods for transmissibility calculation in 2-d.
Definition: 2dtransmissibilitycalculator.hh:44
Dumux::FVMPFALInteractionVolume
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Definition: linteractionvolume.hh:41
Dumux::FVPressure
The finite volume base class for the solution of a pressure equation.
Definition: sequential/cellcentered/pressure.hh:49
Dumux::FVPressure::initialize
void initialize()
Initialize pressure model.
Definition: sequential/cellcentered/pressure.hh:213
Dumux::FVPressure::pressure
PressureSolution & pressure()
Returns the vector containing the pressure solution.
Definition: sequential/cellcentered/pressure.hh:120
Dumux::FVPressure::A_
Matrix A_
Global stiffnes matrix (sparse matrix which is build by the initializeMatrix() function)
Definition: sequential/cellcentered/pressure.hh:324
Dumux::FVPressure::solve
void solve()
Solves the global system of equations to get the spatial distribution of the pressure.
Definition: sequential/cellcentered/pressure.hh:527
Dumux::FVPressure::f_
RHSVector f_
Right hand side vector.
Definition: sequential/cellcentered/pressure.hh:325
properties.hh
Specifies the properties for immiscible 2p diffusion/pressure models.
properties.hh
Properties for a MPFA method.
pressure.hh
Finite Volume Diffusion Model.
Generated by doxygen 1.9.3