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