3.2-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
omethod/2dpressure.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 *****************************************************************************/
22#ifndef DUMUX_FVMPFAO2DPRESSURE2P_HH
23#define DUMUX_FVMPFAO2DPRESSURE2P_HH
24
25// dumux environment
26#include <dumux/porousmediumflow/sequential/cellcentered/pressure.hh>
27#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh>
28#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
29#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
30
31namespace Dumux {
65template<class TypeTag>
66class FvMpfaO2dPressure2p: public FVPressure<TypeTag>
67{
68 using ParentType = FVPressure<TypeTag>;
69 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
70
71 enum
72 {
73 dim = GridView::dimension, dimWorld = GridView::dimensionworld
74 };
75
76 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
77 using Problem = GetPropType<TypeTag, Properties::Problem>;
78
79 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
80
81 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
82 using MaterialLaw = typename SpatialParams::MaterialLaw;
83
84 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
85
86 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
87 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
88
89 using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
90 using CellData = GetPropType<TypeTag, Properties::CellData>;
91 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
92 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
93 using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
94
95 using GridTypeIndices = GetPropType<TypeTag, Properties::GridTypeIndices>;
96
97 enum
98 {
99 pw = Indices::pressureW,
100 pn = Indices::pressureNw,
101 sw = Indices::saturationW,
102 sn = Indices::saturationNw
103 };
104 enum
105 {
106 wPhaseIdx = Indices::wPhaseIdx,
107 nPhaseIdx = Indices::nPhaseIdx,
108 pressureIdx = Indices::pressureIdx,
109 saturationIdx = Indices::saturationIdx,
110 pressEqIdx = Indices::pressureEqIdx,
111 satEqIdx = Indices::satEqIdx,
112 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
113 };
114 enum
115 {
116 globalCorner = 2,
117 globalEdge = 3,
118 neumannNeumann = 0,
119 dirichletDirichlet = 1,
120 dirichletNeumann = 2,
121 neumannDirichlet = 3
122 };
123
124 using Element = typename GridView::Traits::template Codim<0>::Entity;
125 using IntersectionIterator = typename GridView::IntersectionIterator;
126 using Intersection = typename GridView::Intersection;
127
128 using LocalPosition = Dune::FieldVector<Scalar, dim>;
129 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
130 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
131 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
132
133 using DimVector = Dune::FieldVector<Scalar, dim>;
134
135 using InteractionVolume = FVMPFAOInteractionVolume<TypeTag>;
136
137 using GlobalInteractionVolumeVector = std::vector<InteractionVolume>;
138 using InnerBoundaryVolumeFaces = std::vector<Dune::FieldVector<bool, 2*dim> >;
139
141 Intersection getNextIntersection_(const Element&, const IntersectionIterator&);
142
144 friend class FVPressure<TypeTag>;
145 void initializeMatrix();
146
147 void storeInteractionVolumeInfo();
148
150 void assemble();
151public:
152
154 void updateMaterialLaws();
155
161 void updateInteractionVolumeInfo()
162 {
163 interactionVolumes_.clear();
164 innerBoundaryVolumeFaces_.clear();
165
166 interactionVolumes_.resize(problem_.gridView().size(dim), InteractionVolume(problem_.gridView().grid()));
167 innerBoundaryVolumeFaces_.resize(problem_.gridView().size(0), Dune::FieldVector<bool, 2 * dim>(false));
168
169 storeInteractionVolumeInfo();
170 }
171
177 void initialize()
178 {
179 ParentType::initialize();
180
181 const auto element = *problem_.gridView().template begin<0>();
182 FluidState fluidState;
183 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
184 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
185 fluidState.setTemperature(problem_.temperature(element));
186 fluidState.setSaturation(wPhaseIdx, 1.);
187 fluidState.setSaturation(nPhaseIdx, 0.);
188 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
189 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
190 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
191 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
192
193 updateMaterialLaws();
194
195 interactionVolumes_.resize(problem_.gridView().size(dim), InteractionVolume(problem_.gridView().grid()));
196 innerBoundaryVolumeFaces_.resize(problem_.gridView().size(0), Dune::FieldVector<bool, 2 * dim>(false));
197
198 storeInteractionVolumeInfo();
199
200 assemble();
201 this->solve();
202
203 storePressureSolution();
204
205 return;
206 }
207
213 void update()
214 {
215 //error bounds for error term for incompressible models to correct unphysical saturation
216 //over/undershoots due to saturation transport
217 timeStep_ = problem_.timeManager().timeStepSize();
218 maxError_ = 0.0;
219 int size = problem_.gridView().size(0);
220 for (int i = 0; i < size; i++)
221 {
222 Scalar sat = 0;
223 using std::max;
224 switch (saturationType_)
225 {
226 case sw:
227 sat = problem_.variables().cellData(i).saturation(wPhaseIdx);
228 break;
229 case sn:
230 sat = problem_.variables().cellData(i).saturation(nPhaseIdx);
231 break;
232 }
233 if (sat > 1.0)
234 {
235 maxError_ = max(maxError_, (sat - 1.0) / timeStep_);
236 }
237 if (sat < 0.0)
238 {
239 maxError_ = max(maxError_, (-sat) / timeStep_);
240 }
241 }
242
243 assemble();
244 this->solve();
245
246 storePressureSolution();
247
248 return;
249 }
250
254 void storePressureSolution()
255 {
256 // iterate through leaf grid an evaluate c0 at cell center
257 for (const auto& element : elements(problem_.gridView()))
258 {
259 storePressureSolution(element);
260 }
261 }
262
268 void storePressureSolution(const Element& element)
269 {
270 int eIdxGlobal = problem_.variables().index(element);
271 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
272
273 switch (pressureType_)
274 {
275 case pw:
276 {
277 Scalar potW = this->pressure()[eIdxGlobal];
278
279 Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
280 Scalar potPc = cellData.capillaryPressure()
281 + gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
282
283 cellData.setPotential(wPhaseIdx, potW);
284 cellData.setPotential(nPhaseIdx, potW + potPc);
285
286 Scalar pressW = potW - gravityDiff * density_[wPhaseIdx];
287
288 cellData.setPressure(wPhaseIdx, pressW);
289 cellData.setPressure(nPhaseIdx, pressW + cellData.capillaryPressure());
290
291 break;
292 }
293 case pn:
294 {
295 Scalar potNw = this->pressure()[eIdxGlobal];
296
297 Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
298 Scalar potPc = cellData.capillaryPressure()
299 + gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
300
301 cellData.setPotential(nPhaseIdx, potNw);
302 cellData.setPotential(wPhaseIdx, potNw - potPc);
303
304 Scalar pressNw = potNw - gravityDiff * density_[nPhaseIdx];
305
306 cellData.setPressure(wPhaseIdx, pressNw - cellData.capillaryPressure());
307 cellData.setPressure(nPhaseIdx, pressNw);
308
309 break;
310 }
311 }
312
313 cellData.fluxData().resetVelocity();
314 }
315
326 template<class MultiWriter>
327 void addOutputVtkFields(MultiWriter &writer)
328 {
329 int size = problem_.gridView().size(0);
330 ScalarSolutionType *potential = writer.allocateManagedBuffer(size);
331
332 (*potential) = this->pressure();
333
334 if (pressureType_ == pw)
335 {
336 writer.attachCellData(*potential, "wetting potential");
337 }
338
339 if (pressureType_ == pn)
340 {
341 writer.attachCellData(*potential, "nonwetting potential");
342 }
343
344 if (vtkOutputLevel_ > 0)
345 {
346 ScalarSolutionType *pressure = writer.allocateManagedBuffer(size);
347 ScalarSolutionType *pressureSecond = writer.allocateManagedBuffer(size);
348 ScalarSolutionType *potentialSecond = writer.allocateManagedBuffer(size);
349 ScalarSolutionType *pc = writer.allocateManagedBuffer(size);
350
351 for (const auto& element : elements(problem_.gridView()))
352 {
353 int idx = problem_.variables().index(element);
354 CellData& cellData = problem_.variables().cellData(idx);
355
356 (*pc)[idx] = cellData.capillaryPressure();
357
358 if (pressureType_ == pw)
359 {
360 (*pressure)[idx] = cellData.pressure(wPhaseIdx);
361 (*potentialSecond)[idx] = cellData.potential(nPhaseIdx);
362 (*pressureSecond)[idx] = cellData.pressure(nPhaseIdx);
363 }
364
365 if (pressureType_ == pn)
366 {
367 (*pressure)[idx] = cellData.pressure(nPhaseIdx);
368 (*potentialSecond)[idx] = cellData.potential(wPhaseIdx);
369 (*pressureSecond)[idx] = cellData.pressure(wPhaseIdx);
370 }
371 }
372
373 if (pressureType_ == pw)
374 {
375 writer.attachCellData(*pressure, "wetting pressure");
376 writer.attachCellData(*pressureSecond, "nonwetting pressure");
377 writer.attachCellData(*potentialSecond, "nonwetting potential");
378 }
379
380 if (pressureType_ == pn)
381 {
382 writer.attachCellData(*pressure, "nonwetting pressure");
383 writer.attachCellData(*pressureSecond, "wetting pressure");
384 writer.attachCellData(*potentialSecond, "wetting potential");
385 }
386
387 writer.attachCellData(*pc, "capillary pressure");
388 }
389
390 return;
391 }
392
397 FvMpfaO2dPressure2p(Problem& problem) :
398 ParentType(problem), problem_(problem), gravity_(problem.gravity()), maxError_(
399 0.), timeStep_(1.)
400 {
401 if (pressureType_ != pw && pressureType_ != pn)
402 {
403 DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
404 }
405 if (saturationType_ != sw && saturationType_ != sn)
406 {
407 DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
408 }
409 if (getPropValue<TypeTag, Properties::EnableCompressibility>())
410 {
411 DUNE_THROW(Dune::NotImplemented, "Compressibility not supported!");
412 }
413 if (dim != 2)
414 {
415 DUNE_THROW(Dune::NotImplemented, "Dimension not supported!");
416 }
417
418 ErrorTermFactor_ = getParam<Scalar>("Impet.ErrorTermFactor");
419 ErrorTermLowerBound_ = getParam<Scalar>("Impet.ErrorTermLowerBound");
420 ErrorTermUpperBound_ = getParam<Scalar>("Impet.ErrorTermUpperBound");
421
422 density_[wPhaseIdx] = 0.;
423 density_[nPhaseIdx] = 0.;
424 viscosity_[wPhaseIdx] = 0.;
425 viscosity_[nPhaseIdx] = 0.;
426
427 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
428 }
429
430private:
431 Problem& problem_;
432
449 Scalar evaluateErrorTerm_(CellData& cellData)
450 {
451 //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
452 // error reduction routine: volumetric error is damped and inserted to right hand side
453 Scalar sat = 0;
454 switch (saturationType_)
455 {
456 case sw:
457 sat = cellData.saturation(wPhaseIdx);
458 break;
459 case sn:
460 sat = cellData.saturation(nPhaseIdx);
461 break;
462 }
463
464 Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
465 if (sat < 0.0)
466 {
467 error = sat;
468 }
469 error /= timeStep_;
470
471 using std::abs;
472 Scalar errorAbs = abs(error);
473
474 if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
475 && (!problem_.timeManager().willBeFinished()))
476 {
477 return ErrorTermFactor_ * error;
478 }
479 return 0.0;
480 }
481
482protected:
483 GlobalInteractionVolumeVector interactionVolumes_;
484 InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_;
485
486private:
487 const GravityVector& gravity_;
488
489 Scalar maxError_;
490 Scalar timeStep_;
491 Scalar ErrorTermFactor_;
492 Scalar ErrorTermLowerBound_;
493 Scalar ErrorTermUpperBound_;
494
495 Scalar density_[numPhases];
496 Scalar viscosity_[numPhases];
497
498 int vtkOutputLevel_;
499
500 static constexpr Scalar threshold_ = 1e-15;
502 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
504 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
506 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
507};
508
509// TODO doc me!
510template<class TypeTag>
511typename FvMpfaO2dPressure2p<TypeTag>::Intersection
512 FvMpfaO2dPressure2p<TypeTag>::getNextIntersection_(const Element& element,
513 const IntersectionIterator& isIt)
514{
515 auto isItBegin = problem_.gridView().ibegin(element);
516 const auto isEndIt = problem_.gridView().iend(element);
517
518 auto tempIsIt = isIt;
519 auto nextIsIt = ++tempIsIt;
520
521 // get 'nextIsIt'
522 switch (getPropValue<TypeTag, Properties::GridImplementation>())
523 {
524 // for YaspGrid
525 case GridTypeIndices::yaspGrid:
526 {
527 if (nextIsIt == isEndIt)
528 {
529 nextIsIt = isItBegin;
530 }
531 else
532 {
533 nextIsIt = ++tempIsIt;
534
535 if (nextIsIt == isEndIt)
536 {
537 auto tempIsItBegin = isItBegin;
538 nextIsIt = ++tempIsItBegin;
539 }
540 }
541
542 break;
543 }
544 // for ALUGrid and UGGrid
545 case GridTypeIndices::aluGrid:
546 case GridTypeIndices::ugGrid:
547 {
548 if (nextIsIt == isEndIt)
549 nextIsIt = isItBegin;
550
551 break;
552 }
553 default:
554 {
555 DUNE_THROW(Dune::NotImplemented, "GridType can not be used with MPFAO implementation!");
556 break;
557 }
558 }
559
560 return *nextIsIt;
561}
562
563// TODO doc me!
564template<class TypeTag>
565void FvMpfaO2dPressure2p<TypeTag>::initializeMatrix()
566{
567 // determine matrix row sizes
568 for (const auto& element : elements(problem_.gridView()))
569 {
570 // cell index
571 int eIdxGlobalI = problem_.variables().index(element);
572
573 // initialize row size
574 int rowSize = 1;
575
576 // run through all intersections with neighbors
577 const auto isEndIt = problem_.gridView().iend(element);
578 for (auto isIt = problem_.gridView().ibegin(element); isIt != isEndIt; ++isIt)
579 {
580 const auto& intersection = *isIt;
581 auto nextIntersection = getNextIntersection_(element, isIt);
582
583 if (intersection.neighbor())
584 rowSize++;
585
586 // try to find the common neighbor of intersection's and nextIntersection's outside
587 if (intersection.neighbor() && nextIntersection.neighbor())
588 {
589 for (const auto& innerIntersection
590 : intersections(problem_.gridView(), intersection.outside()))
591 for (const auto& innerNextIntersection
592 : intersections(problem_.gridView(), nextIntersection.outside()))
593 {
594 if (innerIntersection.neighbor() && innerNextIntersection.neighbor())
595 {
596 if (innerIntersection.outside() == innerNextIntersection.outside()
597 && innerIntersection.outside() != intersection.inside())
598 {
599 rowSize++;
600 }
601 }
602 }
603 }
604 } // end of intersection loop
605
606 // set number of indices in row eIdxGlobalI to rowSize
607 this->A_.setrowsize(eIdxGlobalI, rowSize);
608
609 } // end of element loop
610
611 // indicate that size of all rows is defined
612 this->A_.endrowsizes();
613
614 // determine position of matrix entries
615 for (const auto& element : elements(problem_.gridView()))
616 {
617 // cell index
618 int eIdxGlobalI = problem_.variables().index(element);
619
620 // add diagonal index
621 this->A_.addindex(eIdxGlobalI, eIdxGlobalI);
622
623 // run through all intersections with neighbors
624 const auto isEndIt = problem_.gridView().iend(element);
625 for (auto isIt = problem_.gridView().ibegin(element); isIt != isEndIt; ++isIt)
626 {
627 const auto& intersection = *isIt;
628 auto nextIntersection = getNextIntersection_(element, isIt);
629
630 if (intersection.neighbor())
631 {
632 // access neighbor
633 int eIdxGlobalJ = problem_.variables().index(intersection.outside());
634
635 // add off diagonal index
636 // add index (row,col) to the matrix
637 this->A_.addindex(eIdxGlobalI, eIdxGlobalJ);
638 }
639
640 if (intersection.neighbor() && nextIntersection.neighbor())
641 {
642 for (const auto& innerIntersection
643 : intersections(problem_.gridView(), intersection.outside()))
644 for (const auto& innerNextIntersection
645 : intersections(problem_.gridView(), nextIntersection.outside()))
646 {
647 if (innerIntersection.neighbor() && innerNextIntersection.neighbor())
648 {
649 auto innerOutside = innerIntersection.outside();
650
651 if (innerOutside == innerNextIntersection.outside()
652 && innerOutside != intersection.inside())
653 {
654 int eIdxGlobalJ = problem_.variables().index(innerOutside);
655
656 this->A_.addindex(eIdxGlobalI, eIdxGlobalJ);
657 }
658 }
659 }
660 }
661 } // end of intersection loop
662 } // end of element loop
663
664 // indicate that all indices are defined, check consistency
665 this->A_.endindices();
666
667 return;
668}
669
670// Indices used in a interaction volume of the MPFA-o method
671// ___________________________________________________
672// | | |
673// | nuxy: cell geometry | nxy: face normal |
674// | vectors (see MPFA) | |
675// | | |
676// | 4-----------3-----------3 |
677// | | --> nu43 | nu34 <-- | |
678// | | |nu41 1|--> n43 ||nu32 |
679// | | v ^ |0 ^ v| |
680// |____________4__0__|n14__|__n23_|_1__2____________|
681// | | 1 0 | 0 |
682// | | ^ |1 nu23 ^ | |
683// | | |nu14 0|--> n12 | | |
684// | | -->nu12 | nu21<-- | |
685// | 1-----------1-----------2 |
686// | elementnumber |inter- |
687// | |face- |
688// | |number |
689// |________________________|________________________|
690
691// only for 2-D general quadrilateral
692// TODO doc me!
693template<class TypeTag>
694void FvMpfaO2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
695{
696 // introduce matrix R for vector rotation and R is initialized as zero matrix
697 DimMatrix R(0);
698
699 // evaluate matrix R
700 if (dim == 2)
701 {
702 R[0][1] = 1;
703 R[1][0] = -1;
704 }
705
706 BoundaryTypes bcType;
707
708 // run through all elements
709 for (const auto& element : elements(problem_.gridView()))
710 {
711 // get common geometry information for the following computation
712
713 int eIdxGlobal1 = problem_.variables().index(element);
714 // get global coordinate of cell 1 center
715 const GlobalPosition& globalPos1 = element.geometry().center();
716
717 // get absolute permeability of neighbor cell 1
718 DimMatrix K1(problem_.spatialParams().intrinsicPermeability(element));
719
720 const auto isIt12End = problem_.gridView().iend(element);
721 for (auto isIt12 = problem_.gridView().ibegin(element); isIt12 != isIt12End; ++isIt12)
722 {
723 const auto& intersection12 = *isIt12;
724 auto intersection14 = getNextIntersection_(element, isIt12);
725
726 int indexInInside12 = intersection12.indexInInside();
727 int indexInInside14 = intersection14.indexInInside();
728
729 // get the intersection node /bar^{x_3} between 'isIt12' and 'isIt14', denoted as 'corner1234'
730 // initialization of corner1234
731
732 const auto referenceElement = ReferenceElements::general(element.type());
733
734 GlobalPosition corner1234(0);
735
736 int globalVertIdx1234 = 0;
737
738 // get the global coordinate and global vertex index of corner1234
739 for (int i = 0; i < intersection12.geometry().corners(); ++i)
740 {
741 bool finished = false;
742
743 const GlobalPosition& isIt12corner = intersection12.geometry().corner(i);
744
745 int localVertIdx12corner = referenceElement.subEntity(indexInInside12, dim - 1, i, dim);
746
747 int globalVertIdx12corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx12corner));
748
749 for (int j = 0; j < intersection14.geometry().corners(); ++j)
750 {
751 int localVertIdx14corner = referenceElement.subEntity(indexInInside14, dim - 1, j, dim);
752
753 int globalVertIdx14corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx14corner));
754
755 if (globalVertIdx12corner == globalVertIdx14corner)
756 {
757 corner1234 = isIt12corner;
758
759 globalVertIdx1234 = globalVertIdx12corner;
760
761 finished = true;
762 break;
763 }
764 }
765
766 if (finished)
767 {
768 break;
769 }
770 }
771
772 if (interactionVolumes_[globalVertIdx1234].isStored())
773 {
774 continue;
775 }
776 else
777 {
778 interactionVolumes_[globalVertIdx1234].setStored();
779 // std::cout << "vertIdx = " << globalVertIdx1234 << "\n";
780 }
781
782 //store pointer 1
783 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element, 0);
784 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection12.indexInInside(), 0, 0);
785 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection14.indexInInside(), 0, 1);
786
787 // center of face in global coordinates, i.e., the midpoint of edge 'isIt12'
788 const GlobalPosition& globalPosFace12 = intersection12.geometry().center();
789
790 // get face volume
791 Scalar faceVol12 = intersection12.geometry().volume() / 2.0;
792
793 // get outer normal vector scaled with half volume of face 'isIt12'
794 DimVector unitOuterNormal12 = intersection12.centerUnitOuterNormal();
795
796 // center of face in global coordinates, i.e., the midpoint of edge 'isIt14'
797 const GlobalPosition& globalPosFace41 = intersection14.geometry().center();
798
799 // get face volume
800 Scalar faceVol41 = intersection14.geometry().volume() / 2.0;
801
802 // get outer normal vector scaled with half volume of face 'isIt14': for numbering of n see Aavatsmark, Eigestad
803 DimVector unitOuterNormal14 = intersection14.centerUnitOuterNormal();
804
805 // compute normal vectors nu14,nu12
806 DimVector nu14(0);
807 R.mv(globalPos1 - globalPosFace12, nu14);
808
809 DimVector nu12(0);
810 R.mv(globalPosFace41 - globalPos1, nu12);
811
812 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu12, K1, 0, 0);
813 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu14, K1, 0, 1);
814 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal12, 0, 0);
815 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 0, 1);
816 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol12, 0, 0);
817 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 0, 1);
818
819 // compute dF1, the area of quadrilateral made by normal vectors 'nu'
820 DimVector Rnu12(0);
821 R.umv(nu12, Rnu12);
822 using std::abs;
823 interactionVolumes_[globalVertIdx1234].setDF(abs(nu14 * Rnu12), 0);
824
825 // handle interior face
826 if (intersection12.neighbor())
827 {
828 // access neighbor cell 2 of 'isIt12'
829 auto element2 = intersection12.outside();
830
831 int eIdxGlobal2 = problem_.variables().index(element2);
832
833 //store pointer 2
834 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element2, 1);
835 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection12.indexInOutside(), 1, 1);
836
837 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal12, 1, 1);
838 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol12, 1, 1);
839
840 // get global coordinate of neighbor cell 2 center
841 const GlobalPosition& globalPos2 = element2.geometry().center();
842
843 // get absolute permeability of neighbor cell 2
844 DimMatrix K2(problem_.spatialParams().intrinsicPermeability(element2));
845
846 // 'isIt14' is an interior face
847 if (intersection14.neighbor())
848 {
849 // neighbor cell 3
850 // access neighbor cell 3
851 auto element4 = intersection14.outside();
852
853 //store pointer 4
854 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element4, 3);
855 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection14.indexInOutside(), 3, 0);
856
857 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 3, 0);
858 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 3, 0);
859
860 // get basic information of cell 1,2's neighbor cell 3,4
861
862 // get global coordinate of neighbor cell 4 center
863 const GlobalPosition& globalPos4 = element4.geometry().center();
864
865 // get absolute permeability of neighbor cell 2
866 DimMatrix K4(problem_.spatialParams().intrinsicPermeability(element4));
867
868 // cell 3
869 GlobalPosition globalPos3(0);
870
871 GlobalPosition globalPosFace23(0);
872 GlobalPosition globalPosFace34(0);
873
874 for (const auto& intersection2
875 : intersections(problem_.gridView(), element2))
876 {
877 bool finished = false;
878
879 for (const auto& intersection4
880 : intersections(problem_.gridView(), element4))
881 {
882 if (intersection2.neighbor() && intersection4.neighbor())
883 {
884 auto element32 = intersection2.outside();
885 auto element34 = intersection4.outside();
886
887 // find the common neighbor cell between cell 2 and cell 3, except cell 1
888 if (element32 == element34 && element32 != element)
889 {
890 //store pointer 3
891 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element32, 2);
892
893 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInInside(), 1,
894 0);
895 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInOutside(), 2,
896 1);
897 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInInside(), 3,
898 1);
899 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInOutside(), 2,
900 0);
901
902 // get global coordinate of neighbor cell 4 center
903 globalPos3 = element32.geometry().center();
904
905 globalPosFace23 = intersection2.geometry().center();
906 globalPosFace34 = intersection4.geometry().center();
907
908 Scalar faceVol23 = intersection2.geometry().volume() / 2.0;
909 Scalar faceVol34 = intersection4.geometry().volume() / 2.0;
910
911 // get outer normal vector scaled with half volume of face : for numbering of n see Aavatsmark, Eigestad
912 DimVector unitOuterNormal23 = intersection2.centerUnitOuterNormal();
913
914 DimVector unitOuterNormal43 = intersection4.centerUnitOuterNormal();
915
916 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 1, 0);
917 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 2, 1);
918 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 2, 0);
919 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 3, 1);
920 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 1, 0);
921 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 2, 1);
922 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 2, 0);
923 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 3, 1);
924
925 // get absolute permeability of neighbor cell 2
926 DimMatrix K3(
927 problem_.spatialParams().intrinsicPermeability(element32));
928
929 // compute normal vectors nu23, nu21; nu32, nu34; nu41, nu43;
930 DimVector nu23(0);
931 R.umv(globalPosFace12 - globalPos2, nu23);
932
933 DimVector nu21(0);
934 R.umv(globalPosFace23 - globalPos2, nu21);
935
936 DimVector nu32(0);
937 R.umv(globalPosFace34 - globalPos3, nu32);
938
939 DimVector nu34(0);
940 R.umv(globalPos3 - globalPosFace23, nu34);
941
942 DimVector nu41(0);
943 R.umv(globalPos4 - globalPosFace34, nu41);
944
945 DimVector nu43(0);
946 R.umv(globalPos4 - globalPosFace41, nu43);
947
948 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu23, K2, 1, 0);
949 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu21, K2, 1, 1);
950 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu34, K3, 2, 0);
951 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu32, K3, 2, 1);
952 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu41, K4, 3, 0);
953 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu43, K4, 3, 1);
954
955 // compute dF2, dF3, dF4 i.e., the area of quadrilateral made by normal vectors 'nu'
956 DimVector Rnu21(0);
957 R.umv(nu21, Rnu21);
958 interactionVolumes_[globalVertIdx1234].setDF(abs(nu23 * Rnu21), 1);
959
960 DimVector Rnu34(0);
961 R.umv(nu34, Rnu34);
962 interactionVolumes_[globalVertIdx1234].setDF(abs(nu32 * Rnu34), 2);
963
964 DimVector Rnu43(0);
965 R.umv(nu43, Rnu43);
966 interactionVolumes_[globalVertIdx1234].setDF(abs(nu41 * Rnu43), 3);
967
968 finished = true;
969
970 break;
971 }
972 }
973 }
974 if (finished)
975 {
976 break;
977 }
978 }
979 }
980
981 // 'isIt14' is on the boundary
982 else
983 {
984 problem_.boundaryTypes(bcType, intersection14);
985 PrimaryVariables boundValues(0.0);
986
987 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 3);
988 if (bcType.isNeumann(pressEqIdx))
989 {
990 problem_.neumann(boundValues, intersection14);
991 boundValues *= faceVol41;
992 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 3);
993 }
994 if (bcType.hasDirichlet())
995 {
996 problem_.dirichlet(boundValues, intersection14);
997 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 3);
998 }
999
1000 // get common geometry information for the following computation
1001 // get the information of the face 'isIt23' between cell2 and cell4 (locally numbered)
1002
1003 // center of face in global coordinates, i.e., the midpoint of edge 'isIt23'
1004 GlobalPosition globalPosFace23(0);
1005
1006 // get face volume
1007 Scalar faceVol23 = 0;
1008
1009 // get outer normal vector scaled with half volume of face 'isIt23'
1010 DimVector unitOuterNormal23(0);
1011
1012 bool finished = false;
1013
1014 for (const auto& intersection2
1015 : intersections(problem_.gridView(), element2))
1016 {
1017 if (intersection2.boundary())
1018 {
1019 for (int i = 0; i < intersection2.geometry().corners(); ++i)
1020 {
1021 int localVertIdx2corner = referenceElement.subEntity(intersection2.indexInInside(), dim - 1, i,
1022 dim);
1023
1024 int globalVertIdx2corner = problem_.variables().index(element2.template subEntity<dim>(localVertIdx2corner));
1025
1026 if (globalVertIdx2corner == globalVertIdx1234)
1027 {
1028 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection2.indexInInside(), 1,
1029 0);
1030
1031 globalPosFace23 = intersection2.geometry().center();
1032
1033 faceVol23 = intersection2.geometry().volume() / 2.0;
1034
1035 unitOuterNormal23 = intersection2.centerUnitOuterNormal();
1036
1037 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal23, 1, 0);
1038 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol23, 1, 0);
1039
1040 problem_.boundaryTypes(bcType, intersection2);
1041 boundValues = 0.0;
1042
1043 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 1);
1044 if (bcType.isNeumann(pressEqIdx))
1045 {
1046 problem_.neumann(boundValues, intersection2);
1047 boundValues *= faceVol23;
1048 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 1);
1049 }
1050 if (bcType.hasDirichlet())
1051 {
1052 problem_.dirichlet(boundValues, intersection2);
1053 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 1);
1054 }
1055
1056 interactionVolumes_[globalVertIdx1234].setOutsideFace(2);
1057
1058 innerBoundaryVolumeFaces_[eIdxGlobal1][intersection12.indexInInside()] = true;
1059 innerBoundaryVolumeFaces_[eIdxGlobal2][intersection12.indexInOutside()] = true;
1060
1061 // compute normal vectors nu23, nu21;
1062 DimVector nu23(0);
1063 R.umv(globalPosFace12 - globalPos2, nu23);
1064
1065 DimVector nu21(0);
1066 R.umv(globalPosFace23 - globalPos2, nu21);
1067
1068 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu23, K2, 1, 0);
1069 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu21, K2, 1, 1);
1070
1071 // compute dF2 i.e., the area of quadrilateral made by normal vectors 'nu'
1072 DimVector Rnu21(0);
1073 R.umv(nu21, Rnu21);
1074 interactionVolumes_[globalVertIdx1234].setDF(abs(nu23 * Rnu21), 1);
1075
1076 finished = true;
1077
1078 break;
1079 }
1080 }
1081 }
1082 if (finished)
1083 {
1084 break;
1085 }
1086 }
1087 if (!finished)
1088 {
1089 DUNE_THROW(
1090 Dune::NotImplemented,
1091 "fvmpfao2pfaboundpressure2p.hh, l. 997: boundary shape not available as interaction volume shape");
1092 }
1093 }
1094 }
1095
1096 // handle boundary face 'isIt12'
1097 else
1098 {
1099 problem_.boundaryTypes(bcType, *isIt12);
1100 PrimaryVariables boundValues(0.0);
1101
1102 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 0);
1103 if (bcType.isNeumann(pressEqIdx))
1104 {
1105 problem_.neumann(boundValues, *isIt12);
1106 boundValues *= faceVol12;
1107 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 0);
1108 }
1109 if (bcType.hasDirichlet())
1110 {
1111 problem_.dirichlet(boundValues, *isIt12);
1112 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 0);
1113 }
1114
1115 // 'isIt14' is on boundary
1116 if (intersection14.boundary())
1117 {
1118 problem_.boundaryTypes(bcType, intersection14);
1119 boundValues = 0.0;
1120
1121 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 3);
1122 if (bcType.isNeumann(pressEqIdx))
1123 {
1124 problem_.neumann(boundValues, intersection14);
1125 boundValues *= faceVol41;
1126 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 3);
1127 }
1128 if (bcType.hasDirichlet())
1129 {
1130 problem_.dirichlet(boundValues, intersection14);
1131 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 3);
1132 }
1133
1134 interactionVolumes_[globalVertIdx1234].setOutsideFace(1);
1135 interactionVolumes_[globalVertIdx1234].setOutsideFace(2);
1136 }
1137
1138 // 'isIt14' is inside
1139 else
1140 {
1141 // neighbor cell 3
1142 // access neighbor cell 3
1143 auto element4 = intersection14.outside();
1144 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection14.indexInOutside(), 3, 0);
1145
1146 //store pointer 4
1147 interactionVolumes_[globalVertIdx1234].setSubVolumeElement(element4, 3);
1148
1149 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal14, 3, 0);
1150 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol41, 3, 0);
1151
1152 // get global coordinate of neighbor cell 3 center
1153 const GlobalPosition& globalPos4 = element4.geometry().center();
1154
1155 int eIdxGlobal4 = problem_.variables().index(element4);
1156
1157 bool finished = false;
1158
1159 // get the information of the face 'isIt34' between cell3 and cell4 (locally numbered)
1160 for (const auto& intersection4
1161 : intersections(problem_.gridView(), element4))
1162 {
1163 if (intersection4.boundary())
1164 {
1165 for (int i = 0; i < intersection4.geometry().corners(); ++i)
1166 {
1167 int localVertIdx4corner = referenceElement.subEntity(intersection4.indexInInside(), dim - 1, i,
1168 dim);
1169
1170 int globalVertIdx4corner = problem_.variables().index(element4.template subEntity<dim>(localVertIdx4corner));
1171
1172 if (globalVertIdx4corner == globalVertIdx1234)
1173 {
1174 interactionVolumes_[globalVertIdx1234].setIndexOnElement(intersection4.indexInInside(), 3,
1175 1);
1176
1177 const GlobalPosition& globalPosFace34 = intersection4.geometry().center();
1178
1179 Scalar faceVol34 = intersection4.geometry().volume() / 2.0;
1180
1181 DimVector unitOuterNormal43 = intersection4.centerUnitOuterNormal();
1182
1183 interactionVolumes_[globalVertIdx1234].setNormal(unitOuterNormal43, 3, 1);
1184 interactionVolumes_[globalVertIdx1234].setFaceArea(faceVol34, 3, 1);
1185
1186 problem_.boundaryTypes(bcType, intersection4);
1187 boundValues = 0.0;
1188
1189 interactionVolumes_[globalVertIdx1234].setBoundary(bcType, 2);
1190 if (bcType.isNeumann(pressEqIdx))
1191 {
1192 problem_.neumann(boundValues, intersection4);
1193 boundValues *= faceVol34;
1194 interactionVolumes_[globalVertIdx1234].setNeumannCondition(boundValues, 2);
1195 }
1196 if (bcType.hasDirichlet())
1197 {
1198 problem_.dirichlet(boundValues, intersection4);
1199 interactionVolumes_[globalVertIdx1234].setDirichletCondition(boundValues, 2);
1200 }
1201
1202 interactionVolumes_[globalVertIdx1234].setOutsideFace(1);
1203
1204 innerBoundaryVolumeFaces_[eIdxGlobal1][intersection14.indexInInside()] = true;
1205 innerBoundaryVolumeFaces_[eIdxGlobal4][intersection14.indexInOutside()] = true;
1206
1207 // get absolute permeability of neighbor cell 2
1208 DimMatrix K4(
1209 problem_.spatialParams().intrinsicPermeability(element4));
1210
1211 // compute normal vectors nu41, nu43;
1212 DimVector nu41(0);
1213 R.umv(globalPos4 - globalPosFace34, nu41);
1214
1215 DimVector nu43(0);
1216 R.umv(globalPos4 - globalPosFace41, nu43);
1217
1218 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu41, K4, 3, 0);
1219 interactionVolumes_[globalVertIdx1234].setPermTimesNu(nu43, K4, 3, 1);
1220
1221 // compute dF1, dF3 i.e., the area of quadrilateral made by normal vectors 'nu'
1222 DimVector Rnu43(0);
1223 R.umv(nu43, Rnu43);
1224 interactionVolumes_[globalVertIdx1234].setDF(abs(nu41 * Rnu43), 3);
1225
1226 finished = true;
1227
1228 break;
1229 }
1230 }
1231 }
1232 if (finished)
1233 {
1234 break;
1235 }
1236 }
1237 if (!finished)
1238 {
1239 DUNE_THROW(
1240 Dune::NotImplemented,
1241 "fvmpfao2pfaboundpressure2p.hh, l. 1164: boundary shape not available as interaction volume shape");
1242 }
1243 }
1244 }
1245
1246 } // end all intersections
1247 } // end grid traversal
1248
1249 return;
1250}
1251
1252// only for 2-D general quadrilateral
1253//TODO doc me!
1254template<class TypeTag>
1255void FvMpfaO2dPressure2p<TypeTag>::assemble()
1256{
1257 // initialization: set global matrix this->A_ to zero
1258 this->A_ = 0;
1259 this->f_ = 0;
1260
1261 // run through all vertices
1262 for (const auto& vertex : vertices(problem_.gridView()))
1263 {
1264 int vIdxGlobal = problem_.variables().index(vertex);
1265
1266 InteractionVolume& interactionVolume = interactionVolumes_[vIdxGlobal];
1267
1268 if (interactionVolume.isInnerVolume())
1269 {
1270
1271 auto element1 = interactionVolume.getSubVolumeElement(0);
1272 auto element2 = interactionVolume.getSubVolumeElement(1);
1273 auto element3 = interactionVolume.getSubVolumeElement(2);
1274 auto element4 = interactionVolume.getSubVolumeElement(3);
1275
1276 // get global coordinate of cell centers
1277 const GlobalPosition& globalPos1 = element1.geometry().center();
1278 const GlobalPosition& globalPos2 = element2.geometry().center();
1279 const GlobalPosition& globalPos3 = element3.geometry().center();
1280 const GlobalPosition& globalPos4 = element4.geometry().center();
1281
1282 // cell volumes
1283 Scalar volume1 = element1.geometry().volume();
1284 Scalar volume2 = element2.geometry().volume();
1285 Scalar volume3 = element3.geometry().volume();
1286 Scalar volume4 = element4.geometry().volume();
1287
1288 // cell index
1289 int eIdxGlobal1 = problem_.variables().index(element1);
1290 int eIdxGlobal2 = problem_.variables().index(element2);
1291 int eIdxGlobal3 = problem_.variables().index(element3);
1292 int eIdxGlobal4 = problem_.variables().index(element4);
1293
1294 //get the cell Data
1295 CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
1296 CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
1297 CellData& cellData3 = problem_.variables().cellData(eIdxGlobal3);
1298 CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
1299
1300 // evaluate right hand side
1301 PrimaryVariables source(0.0);
1302 problem_.source(source, element1);
1303 this->f_[eIdxGlobal1] += volume1 / (4.0)
1304 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1305 problem_.source(source, element2);
1306 this->f_[eIdxGlobal2] += volume2 / (4.0)
1307 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1308 problem_.source(source, element3);
1309 this->f_[eIdxGlobal3] += volume3 / (4.0)
1310 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1311 problem_.source(source, element4);
1312 this->f_[eIdxGlobal4] += volume4 / (4.0)
1313 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1314
1315 this->f_[eIdxGlobal1] += evaluateErrorTerm_(cellData1) * volume1 / (4.0);
1316 this->f_[eIdxGlobal2] += evaluateErrorTerm_(cellData2) * volume2 / (4.0);
1317 this->f_[eIdxGlobal3] += evaluateErrorTerm_(cellData3) * volume3 / (4.0);
1318 this->f_[eIdxGlobal4] += evaluateErrorTerm_(cellData4) * volume4 / (4.0);
1319
1320 //get mobilities of the phases
1321 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
1322 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
1323
1324 //compute total mobility of cell 1
1325 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
1326
1327 //get mobilities of the phases
1328 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
1329 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
1330
1331 //compute total mobility of cell 1
1332 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
1333
1334 //get mobilities of the phases
1335 Dune::FieldVector<Scalar, numPhases> lambda3(cellData3.mobility(wPhaseIdx));
1336 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
1337
1338 //compute total mobility of cell 1
1339 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
1340
1341 //get mobilities of the phases
1342 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
1343 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
1344
1345 //compute total mobility of cell 1
1346 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
1347
1348 Scalar gn12nu14 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 0, 1);
1349 Scalar gn12nu12 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 0, 0);
1350 Scalar gn14nu14 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 1, 1);
1351 Scalar gn14nu12 = interactionVolume.getNtkrkNu_df(lambdaTotal1, 0, 1, 0);
1352 Scalar gn12nu23 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 1, 0);
1353 Scalar gn12nu21 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 1, 1);
1354 Scalar gn23nu23 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 0, 0);
1355 Scalar gn23nu21 = interactionVolume.getNtkrkNu_df(lambdaTotal2, 1, 0, 1);
1356 Scalar gn43nu32 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 0, 1);
1357 Scalar gn43nu34 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 0, 0);
1358 Scalar gn23nu32 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 1, 1);
1359 Scalar gn23nu34 = interactionVolume.getNtkrkNu_df(lambdaTotal3, 2, 1, 0);
1360 Scalar gn43nu41 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 1, 0);
1361 Scalar gn43nu43 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 1, 1);
1362 Scalar gn14nu41 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 0, 0);
1363 Scalar gn14nu43 = interactionVolume.getNtkrkNu_df(lambdaTotal4, 3, 0, 1);
1364
1365 // compute transmissibility matrix T = CA^{-1}B+F
1366 Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> C(0), F(0), A(0), B(0);
1367
1368 // evaluate matrix C, F, A, B
1369 C[0][0] = -gn12nu12;
1370 C[0][3] = -gn12nu14;
1371 C[1][0] = gn23nu21;
1372 C[1][1] = -gn23nu23;
1373 C[2][1] = gn43nu32;
1374 C[2][2] = gn43nu34;
1375 C[3][2] = -gn14nu43;
1376 C[3][3] = gn14nu41;
1377
1378 F[0][0] = gn12nu12 + gn12nu14;
1379 F[1][1] = -gn23nu21 + gn23nu23;
1380 F[2][2] = -gn43nu34 - gn43nu32;
1381 F[3][3] = gn14nu43 - gn14nu41;
1382
1383 A[0][0] = gn12nu12 + gn12nu21;
1384 A[0][1] = -gn12nu23;
1385 A[0][3] = gn12nu14;
1386 A[1][0] = -gn23nu21;
1387 A[1][1] = gn23nu23 + gn23nu32;
1388 A[1][2] = gn23nu34;
1389 A[2][1] = -gn43nu32;
1390 A[2][2] = -gn43nu34 - gn43nu43;
1391 A[2][3] = gn43nu41;
1392 A[3][0] = -gn14nu12;
1393 A[3][2] = gn14nu43;
1394 A[3][3] = -gn14nu41 - gn14nu14;
1395
1396 B[0][0] = gn12nu12 + gn12nu14;
1397 B[0][1] = gn12nu21 - gn12nu23;
1398 B[1][1] = -gn23nu21 + gn23nu23;
1399 B[1][2] = gn23nu34 + gn23nu32;
1400 B[2][2] = -gn43nu34 - gn43nu32;
1401 B[2][3] = -gn43nu43 + gn43nu41;
1402 B[3][0] = -gn14nu12 - gn14nu14;
1403 B[3][3] = gn14nu43 - gn14nu41;
1404
1405 // compute T
1406 A.invert();
1407 F += C.rightmultiply(B.leftmultiply(A));
1408 Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> T(F);
1409
1410 // std::cout<<"Tpress = "<<T<<"\n";
1411
1412 // assemble the global matrix this->A_ and right hand side f
1413 this->A_[eIdxGlobal1][eIdxGlobal1] += T[0][0] + T[3][0];
1414 this->A_[eIdxGlobal1][eIdxGlobal2] += T[0][1] + T[3][1];
1415 this->A_[eIdxGlobal1][eIdxGlobal3] += T[0][2] + T[3][2];
1416 this->A_[eIdxGlobal1][eIdxGlobal4] += T[0][3] + T[3][3];
1417
1418 this->A_[eIdxGlobal2][eIdxGlobal1] += -T[0][0] + T[1][0];
1419 this->A_[eIdxGlobal2][eIdxGlobal2] += -T[0][1] + T[1][1];
1420 this->A_[eIdxGlobal2][eIdxGlobal3] += -T[0][2] + T[1][2];
1421 this->A_[eIdxGlobal2][eIdxGlobal4] += -T[0][3] + T[1][3];
1422
1423 this->A_[eIdxGlobal3][eIdxGlobal1] -= T[1][0] + T[2][0];
1424 this->A_[eIdxGlobal3][eIdxGlobal2] -= T[1][1] + T[2][1];
1425 this->A_[eIdxGlobal3][eIdxGlobal3] -= T[1][2] + T[2][2];
1426 this->A_[eIdxGlobal3][eIdxGlobal4] -= T[1][3] + T[2][3];
1427
1428 this->A_[eIdxGlobal4][eIdxGlobal1] += T[2][0] - T[3][0];
1429 this->A_[eIdxGlobal4][eIdxGlobal2] += T[2][1] - T[3][1];
1430 this->A_[eIdxGlobal4][eIdxGlobal3] += T[2][2] - T[3][2];
1431 this->A_[eIdxGlobal4][eIdxGlobal4] += T[2][3] - T[3][3];
1432
1433 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
1434 {
1435 this->A_[eIdxGlobal1][eIdxGlobal1] += T[0][0];
1436 this->A_[eIdxGlobal1][eIdxGlobal2] += T[0][1];
1437 this->A_[eIdxGlobal1][eIdxGlobal3] += T[0][2];
1438 this->A_[eIdxGlobal1][eIdxGlobal4] += T[0][3];
1439 }
1440 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
1441 {
1442 this->A_[eIdxGlobal1][eIdxGlobal1] += T[3][0];
1443 this->A_[eIdxGlobal1][eIdxGlobal2] += T[3][1];
1444 this->A_[eIdxGlobal1][eIdxGlobal3] += T[3][2];
1445 this->A_[eIdxGlobal1][eIdxGlobal4] += T[3][3];
1446
1447 }
1448 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
1449 {
1450 this->A_[eIdxGlobal2][eIdxGlobal1] += T[1][0];
1451 this->A_[eIdxGlobal2][eIdxGlobal2] += T[1][1];
1452 this->A_[eIdxGlobal2][eIdxGlobal3] += T[1][2];
1453 this->A_[eIdxGlobal2][eIdxGlobal4] += T[1][3];
1454 }
1455 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
1456 {
1457 this->A_[eIdxGlobal2][eIdxGlobal1] += -T[0][0];
1458 this->A_[eIdxGlobal2][eIdxGlobal2] += -T[0][1];
1459 this->A_[eIdxGlobal2][eIdxGlobal3] += -T[0][2];
1460 this->A_[eIdxGlobal2][eIdxGlobal4] += -T[0][3];
1461 }
1462 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
1463 {
1464 this->A_[eIdxGlobal3][eIdxGlobal1] -= T[2][0];
1465 this->A_[eIdxGlobal3][eIdxGlobal2] -= T[2][1];
1466 this->A_[eIdxGlobal3][eIdxGlobal3] -= T[2][2];
1467 this->A_[eIdxGlobal3][eIdxGlobal4] -= T[2][3];
1468 }
1469 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
1470 {
1471 this->A_[eIdxGlobal3][eIdxGlobal1] -= T[1][0];
1472 this->A_[eIdxGlobal3][eIdxGlobal2] -= T[1][1];
1473 this->A_[eIdxGlobal3][eIdxGlobal3] -= T[1][2];
1474 this->A_[eIdxGlobal3][eIdxGlobal4] -= T[1][3];
1475 }
1476 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
1477 {
1478 this->A_[eIdxGlobal4][eIdxGlobal1] += -T[3][0];
1479 this->A_[eIdxGlobal4][eIdxGlobal2] += -T[3][1];
1480 this->A_[eIdxGlobal4][eIdxGlobal3] += -T[3][2];
1481 this->A_[eIdxGlobal4][eIdxGlobal4] += -T[3][3];
1482 }
1483 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
1484 {
1485 this->A_[eIdxGlobal4][eIdxGlobal1] += T[2][0];
1486 this->A_[eIdxGlobal4][eIdxGlobal2] += T[2][1];
1487 this->A_[eIdxGlobal4][eIdxGlobal3] += T[2][2];
1488 this->A_[eIdxGlobal4][eIdxGlobal4] += T[2][3];
1489 }
1490
1491 //add capillary pressure and gravity terms to right-hand-side
1492 //calculate capillary pressure velocity
1493 Dune::FieldVector<Scalar, 2 * dim> pc(0);
1494 pc[0] = cellData1.capillaryPressure();
1495 pc[1] = cellData2.capillaryPressure();
1496 pc[2] = cellData3.capillaryPressure();
1497 pc[3] = cellData4.capillaryPressure();
1498
1499 Dune::FieldVector<Scalar, 2 * dim> gravityDiff(0);
1500
1501// std::cout<<"maxPos = "<<problem_.bBoxMax()<<"\n";
1502
1503 gravityDiff[0] = (problem_.bBoxMax() - globalPos1) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1504 gravityDiff[1] = (problem_.bBoxMax() - globalPos2) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1505 gravityDiff[2] = (problem_.bBoxMax() - globalPos3) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1506 gravityDiff[3] = (problem_.bBoxMax() - globalPos4) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1507
1508 pc += gravityDiff;
1509
1510 if (pc[0] == 0 && pc[1] == 0 && pc[2] == 0 && pc[3] == 0)
1511 {
1512 continue;
1513 }
1514
1515 Dune::FieldVector<Scalar, 2 * dim> pcFlux(0);
1516
1517 T.mv(pc, pcFlux);
1518
1519 // std::cout<<"pcFlux = "<<pcFlux<<"\n";
1520
1521 Scalar pcPotential12 = pcFlux[0];
1522 Scalar pcPotential14 = pcFlux[3];
1523 Scalar pcPotential32 = -pcFlux[1];
1524 Scalar pcPotential34 = -pcFlux[2];
1525
1526 // std::cout<<"pcPotential12 = "<<pcPotential12<<"\n";
1527
1528 //compute mobilities of face 1
1529 Dune::FieldVector<Scalar, numPhases> lambda12Upw(0.0);
1530 lambda12Upw[wPhaseIdx] = (pcPotential12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
1531 lambda12Upw[nPhaseIdx] = (pcPotential12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
1532
1533 //compute mobilities of face 4
1534 Dune::FieldVector<Scalar, numPhases> lambda14Upw(0.0);
1535 lambda14Upw[wPhaseIdx] = (pcPotential14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
1536 lambda14Upw[nPhaseIdx] = (pcPotential14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
1537
1538 //compute mobilities of face 2
1539 Dune::FieldVector<Scalar, numPhases> lambda32Upw(0.0);
1540 lambda32Upw[wPhaseIdx] = (pcPotential32 >= 0) ? lambda3[wPhaseIdx] : lambda2[wPhaseIdx];
1541 lambda32Upw[nPhaseIdx] = (pcPotential32 >= 0) ? lambda3[nPhaseIdx] : lambda2[nPhaseIdx];
1542
1543 //compute mobilities of face 3
1544 Dune::FieldVector<Scalar, numPhases> lambda34Upw(0.0);
1545 lambda34Upw[wPhaseIdx] = (pcPotential34 >= 0) ? lambda3[wPhaseIdx] : lambda4[wPhaseIdx];
1546 lambda34Upw[nPhaseIdx] = (pcPotential34 >= 0) ? lambda3[nPhaseIdx] : lambda4[nPhaseIdx];
1547
1548 for (int i = 0; i < numPhases; i++)
1549 {
1550 Scalar lambdaT12 = lambda12Upw[wPhaseIdx] + lambda12Upw[nPhaseIdx];
1551 Scalar lambdaT14 = lambda14Upw[wPhaseIdx] + lambda14Upw[nPhaseIdx];
1552 Scalar lambdaT32 = lambda32Upw[wPhaseIdx] + lambda32Upw[nPhaseIdx];
1553 Scalar lambdaT34 = lambda34Upw[wPhaseIdx] + lambda34Upw[nPhaseIdx];
1554 Scalar fracFlow12 = (lambdaT12 > threshold_) ? lambda12Upw[i] / (lambdaT12) : 0.0;
1555 Scalar fracFlow14 = (lambdaT14 > threshold_) ? lambda14Upw[i] / (lambdaT14) : 0.0;
1556 Scalar fracFlow32 = (lambdaT32 > threshold_) ? lambda32Upw[i] / (lambdaT32) : 0.0;
1557 Scalar fracFlow34 = (lambdaT34 > threshold_) ? lambda34Upw[i] / (lambdaT34) : 0.0;
1558
1559 Dune::FieldVector<Scalar, 2 * dim> pcFluxReal(pcFlux);
1560
1561 pcFluxReal[0] *= fracFlow12;
1562 pcFluxReal[1] *= fracFlow32;
1563 pcFluxReal[2] *= fracFlow34;
1564 pcFluxReal[3] *= fracFlow14;
1565
1566// using std::isnan;
1567// if (isnan(pcFluxReal.two_norm()))
1568// std::cout<<"pcFlux = "<<pcFlux<<"\n";
1569
1570 switch (pressureType_)
1571 {
1572 case pw:
1573 {
1574 if (i == nPhaseIdx)
1575 {
1576 //add capillary pressure term to right hand side
1577 this->f_[eIdxGlobal1] -= (pcFluxReal[0] + pcFluxReal[3]);
1578 this->f_[eIdxGlobal2] -= (pcFluxReal[1] - pcFluxReal[0]);
1579 this->f_[eIdxGlobal3] -= (-pcFluxReal[2] - pcFluxReal[1]);
1580 this->f_[eIdxGlobal4] -= (-pcFluxReal[3] + pcFluxReal[2]);
1581
1582 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
1583 {
1584 this->f_[eIdxGlobal1] -= pcFluxReal[0];
1585 }
1586 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
1587 {
1588 this->f_[eIdxGlobal1] -= pcFluxReal[3];
1589 }
1590 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
1591 {
1592 this->f_[eIdxGlobal2] -= pcFluxReal[1];
1593 }
1594 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
1595 {
1596 this->f_[eIdxGlobal2] += pcFluxReal[0];
1597 }
1598 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
1599 {
1600 this->f_[eIdxGlobal3] += pcFluxReal[2];
1601 }
1602 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
1603 {
1604 this->f_[eIdxGlobal3] += pcFluxReal[1];
1605 }
1606 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
1607 {
1608 this->f_[eIdxGlobal4] += pcFluxReal[3];
1609 }
1610 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
1611 {
1612 this->f_[eIdxGlobal4] -= pcFluxReal[2];
1613 }
1614 }
1615 break;
1616 }
1617 case pn:
1618 {
1619 if (i == wPhaseIdx)
1620 {
1621 //add capillary pressure term to right hand side
1622 this->f_[eIdxGlobal1] += (pcFluxReal[0] + pcFluxReal[1]);
1623 this->f_[eIdxGlobal2] += (pcFluxReal[1] - pcFluxReal[0]);
1624 this->f_[eIdxGlobal3] += (-pcFluxReal[2] - pcFluxReal[1]);
1625 this->f_[eIdxGlobal4] += (-pcFluxReal[3] + pcFluxReal[2]);
1626
1627 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 0)])
1628 {
1629 this->f_[eIdxGlobal1] += pcFluxReal[0];
1630 }
1631 if (innerBoundaryVolumeFaces_[eIdxGlobal1][interactionVolume.getIndexOnElement(0, 1)])
1632 {
1633 this->f_[eIdxGlobal1] += pcFluxReal[3];
1634 }
1635 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 0)])
1636 {
1637 this->f_[eIdxGlobal2] += pcFluxReal[1];
1638 }
1639 if (innerBoundaryVolumeFaces_[eIdxGlobal2][interactionVolume.getIndexOnElement(1, 1)])
1640 {
1641 this->f_[eIdxGlobal2] -= pcFluxReal[0];
1642 }
1643 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 0)])
1644 {
1645 this->f_[eIdxGlobal3] -= pcFluxReal[2];
1646 }
1647 if (innerBoundaryVolumeFaces_[eIdxGlobal3][interactionVolume.getIndexOnElement(2, 1)])
1648 {
1649 this->f_[eIdxGlobal3] -= pcFluxReal[1];
1650 }
1651 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 0)])
1652 {
1653 this->f_[eIdxGlobal4] -= pcFluxReal[3];
1654 }
1655 if (innerBoundaryVolumeFaces_[eIdxGlobal4][interactionVolume.getIndexOnElement(3, 1)])
1656 {
1657 this->f_[eIdxGlobal4] += pcFluxReal[2];
1658 }
1659 }
1660 break;
1661 }
1662 }
1663 }
1664 }
1665
1666 // at least one face on boundary!
1667 else
1668 {
1669 for (int elemIdx = 0; elemIdx < 2 * dim; elemIdx++)
1670 {
1671 bool isOutside = false;
1672 for (int fIdx = 0; fIdx < dim; fIdx++)
1673 {
1674 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
1675 if (interactionVolume.isOutsideFace(intVolFaceIdx))
1676 {
1677 isOutside = true;
1678 break;
1679 }
1680 }
1681 if (isOutside)
1682 {
1683 continue;
1684 }
1685
1686 auto element = interactionVolume.getSubVolumeElement(elemIdx);
1687
1688 // get global coordinate of cell centers
1689 const GlobalPosition& globalPos = element.geometry().center();
1690
1691 // cell volumes
1692 Scalar volume = element.geometry().volume();
1693
1694 // cell index
1695 int eIdxGlobal = problem_.variables().index(element);
1696
1697 //get the cell Data
1698 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
1699
1700 //permeability vector at boundary
1701 DimMatrix permeability(problem_.spatialParams().intrinsicPermeability(element));
1702
1703 // evaluate right hand side
1704 PrimaryVariables source(0);
1705 problem_.source(source, element);
1706 this->f_[eIdxGlobal] += volume / (4.0)
1707 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
1708
1709 this->f_[eIdxGlobal] += evaluateErrorTerm_(cellData) * volume / (4.0);
1710
1711 //get mobilities of the phases
1712 Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
1713 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
1714
1715 Scalar pc = cellData.capillaryPressure();
1716
1717 Scalar gravityDiff = (problem_.bBoxMax() - globalPos) * gravity_
1718 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1719
1720 pc += gravityDiff; //minus because of gravity definition!
1721
1722 for (int fIdx = 0; fIdx < dim; fIdx++)
1723 {
1724 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
1725
1726 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
1727 {
1728
1729 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressEqIdx))
1730 {
1731 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
1732
1733 const auto referenceElement = ReferenceElements::general(element.type());
1734
1735 const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
1736
1737 const GlobalPosition& globalPosFace = element.geometry().global(localPos);
1738
1739 DimVector distVec(globalPosFace - globalPos);
1740 Scalar dist = distVec.two_norm();
1741 DimVector unitDistVec(distVec);
1742 unitDistVec /= dist;
1743
1744 Scalar faceArea = interactionVolume.getFaceArea(elemIdx, fIdx);
1745
1746 // get pc and lambda at the boundary
1747 Scalar satWBound = cellData.saturation(wPhaseIdx);
1748 //check boundary sat at face 1
1749 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
1750 {
1751 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
1752 switch (saturationType_)
1753 {
1754 case sw:
1755 {
1756 satWBound = satBound;
1757 break;
1758 }
1759 case sn:
1760 {
1761 satWBound = 1 - satBound;
1762 break;
1763 }
1764 }
1765
1766 }
1767
1768 Scalar pcBound = MaterialLaw::pc(
1769 problem_.spatialParams().materialLawParams(element), satWBound);
1770
1771 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
1772 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1773
1774 pcBound += gravityDiffBound;
1775
1776 Dune::FieldVector<Scalar, numPhases> lambdaBound(
1777 MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
1778 satWBound));
1779 lambdaBound[nPhaseIdx] = MaterialLaw::krn(
1780 problem_.spatialParams().materialLawParams(element), satWBound);
1781 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
1782 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
1783
1784 Scalar potentialBound = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx];
1785 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
1786
1787 //calculate potential gradients
1788 Scalar potentialDiffW = 0;
1789 Scalar potentialDiffNw = 0;
1790 switch (pressureType_)
1791 {
1792 case pw:
1793 {
1794 potentialBound += density_[wPhaseIdx]*gdeltaZ;
1795 potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound) / dist;
1796 potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound - pcBound) / dist;
1797 break;
1798 }
1799 case pn:
1800 {
1801 potentialBound += density_[nPhaseIdx]*gdeltaZ;
1802 potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound + pcBound) / dist;
1803 potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound) / dist;
1804 break;
1805 }
1806 }
1807
1808 Scalar lambdaTotal = (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
1809 lambdaTotal += (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
1810
1811 DimVector permTimesNormal(0);
1812 permeability.mv(unitDistVec, permTimesNormal);
1813
1814 //calculate current matrix entry
1815 Scalar entry = lambdaTotal * (unitDistVec * permTimesNormal) / dist * faceArea;
1816
1817 //calculate right hand side
1818 Scalar pcFlux = 0;
1819
1820 switch (pressureType_)
1821 {
1822 case pw:
1823 {
1824 // calculate capillary pressure gradient
1825 DimVector pcGradient = unitDistVec;
1826 pcGradient *= (pc - pcBound) / dist;
1827
1828 //add capillary pressure term to right hand side
1829 pcFlux = 0.5 * (lambda[nPhaseIdx] + lambdaBound[nPhaseIdx])
1830 * (permTimesNormal * pcGradient) * faceArea;
1831
1832 break;
1833 }
1834 case pn:
1835 {
1836 // calculate capillary pressure gradient
1837 DimVector pcGradient = unitDistVec;
1838 pcGradient *= (pc - pcBound) / dist;
1839
1840 //add capillary pressure term to right hand side
1841 pcFlux = 0.5 * (lambda[wPhaseIdx] + lambdaBound[wPhaseIdx])
1842 * (permTimesNormal * pcGradient) * faceArea;
1843
1844 break;
1845
1846 }
1847 }
1848
1849 // set diagonal entry and right hand side entry
1850 this->A_[eIdxGlobal][eIdxGlobal] += entry;
1851 this->f_[eIdxGlobal] += entry * potentialBound;
1852
1853 if (pc == 0 && pcBound == 0)
1854 {
1855 continue;
1856 }
1857
1858 for (int i = 0; i < numPhases; i++)
1859 {
1860 switch (pressureType_)
1861 {
1862 case pw:
1863 {
1864 if (i == nPhaseIdx)
1865 {
1866 //add capillary pressure term to right hand side
1867 this->f_[eIdxGlobal] -= pcFlux;
1868 }
1869 break;
1870 }
1871 case pn:
1872 {
1873 if (i == wPhaseIdx)
1874 {
1875 //add capillary pressure term to right hand side
1876 this->f_[eIdxGlobal] += pcFlux;
1877 }
1878 break;
1879 }
1880 }
1881 }
1882
1883 }
1884 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx))
1885 {
1886 Scalar J = interactionVolume.getNeumannValues(intVolFaceIdx)[wPhaseIdx]
1887 / density_[wPhaseIdx];
1888 J += interactionVolume.getNeumannValues(intVolFaceIdx)[nPhaseIdx] / density_[nPhaseIdx];
1889 this->f_[eIdxGlobal] -= J;
1890 }
1891 else
1892 {
1893 std::cout << "interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx)"
1894 << interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressEqIdx) << "\n";
1895 DUNE_THROW(Dune::NotImplemented,
1896 "No valid boundary condition type defined for pressure equation!");
1897 }
1898 }
1899 }
1900 }
1901
1902 } // end boundaries
1903
1904 } // end vertex iterator
1905
1906 // only do more if we have more than one process
1907 if (problem_.gridView().comm().size() > 1)
1908 {
1909 // set ghost and overlap element entries
1910 for (const auto& element : elements(problem_.gridView()))
1911 {
1912 if (element.partitionType() == Dune::InteriorEntity)
1913 continue;
1914
1915 // get the global index of the cell
1916 int eIdxGlobalI = problem_.variables().index(element);
1917
1918 this->A_[eIdxGlobalI] = 0.0;
1919 this->A_[eIdxGlobalI][eIdxGlobalI] = 1.0;
1920 this->f_[eIdxGlobalI] = this->pressure()[eIdxGlobalI];
1921 }
1922 }
1923
1924 return;
1925}
1926
1927
1933template<class TypeTag>
1934void FvMpfaO2dPressure2p<TypeTag>::updateMaterialLaws()
1935{
1936 // iterate through leaf grid an evaluate c0 at cell center
1937 for (const auto& element : elements(problem_.gridView()))
1938 {
1939 int eIdxGlobal = problem_.variables().index(element);
1940
1941 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
1942
1943 Scalar satW = cellData.saturation(wPhaseIdx);
1944
1945 Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
1946
1947 cellData.setCapillaryPressure(pc);
1948
1949 // initialize mobilities
1950 Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
1951 / viscosity_[wPhaseIdx];
1952 Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
1953 / viscosity_[nPhaseIdx];
1954
1955 // initialize mobilities
1956 cellData.setMobility(wPhaseIdx, mobilityW);
1957 cellData.setMobility(nPhaseIdx, mobilityNw);
1958
1959 //initialize fractional flow functions
1960 cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
1961 cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
1962 }
1963 return;
1964}
1965
1966} // end namespace Dumux
1967#endif
ointeractionvolume.hh
Class including the information of an interaction volume of a MPFA O-method that does not change with...
Dumux
Definition: adapt.hh:29
Dumux::GetProp
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
Dumux::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Dumux::IOName::viscosity
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
Dumux::IOName::permeability
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
Dumux::IOName::pressure
std::string pressure(int phaseIdx) noexcept
I/O name of pressure for multiphase systems.
Definition: name.hh:34
Dumux::IOName::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Dumux::FvMpfaO2dPressure2p
Finite volume MPFA O-method discretization of a two-phase flow pressure equation of the sequential IM...
Definition: omethod/2dpressure.hh:67
Dumux::FvMpfaO2dPressure2p::updateMaterialLaws
void updateMaterialLaws()
Constitutive functions are initialized and stored in the variables object.
Definition: omethod/2dpressure.hh:1934
Dumux::FvMpfaO2dPressure2p::interactionVolumes_
GlobalInteractionVolumeVector interactionVolumes_
Definition: omethod/2dpressure.hh:483
Dumux::FvMpfaO2dPressure2p::storePressureSolution
void storePressureSolution(const Element &element)
Stores the pressure solution of a cell.
Definition: omethod/2dpressure.hh:268
Dumux::FvMpfaO2dPressure2p::innerBoundaryVolumeFaces_
InnerBoundaryVolumeFaces innerBoundaryVolumeFaces_
Definition: omethod/2dpressure.hh:484
Dumux::FvMpfaO2dPressure2p::updateInteractionVolumeInfo
void updateInteractionVolumeInfo()
Updates interaction volumes.
Definition: omethod/2dpressure.hh:161
Dumux::FvMpfaO2dPressure2p::addOutputVtkFields
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: omethod/2dpressure.hh:327
Dumux::FvMpfaO2dPressure2p::update
void update()
Pressure update.
Definition: omethod/2dpressure.hh:213
Dumux::FvMpfaO2dPressure2p::initialize
void initialize()
Initializes the pressure model.
Definition: omethod/2dpressure.hh:177
Dumux::FvMpfaO2dPressure2p::storePressureSolution
void storePressureSolution()
Globally stores the pressure solution.
Definition: omethod/2dpressure.hh:254
Dumux::FvMpfaO2dPressure2p::FvMpfaO2dPressure2p
FvMpfaO2dPressure2p(Problem &problem)
Constructs a FvMpfaO2dPressure2p object.
Definition: omethod/2dpressure.hh:397
Dumux::FVMPFAOInteractionVolume
Class including the information of an interaction volume of a MPFA O-method that does not change with...
Definition: ointeractionvolume.hh:38
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::solve
void solve()
Solves the global system of equations to get the spatial distribution of the pressure.
Definition: sequential/cellcentered/pressure.hh:527
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