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