3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
saturation.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_FVSATURATION2P_HH
25#define DUMUX_FVSATURATION2P_HH
26
29
30namespace Dumux {
31
74template<class TypeTag>
75class FVSaturation2P: public FVTransport<TypeTag>
76{
78 using Implementation = typename GET_PROP_TYPE(TypeTag, TransportModel);
79
80 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
81
82 enum
83 {
84 dim = GridView::dimension, dimWorld = GridView::dimensionworld
85 };
86
87 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
88 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
89
90 using Velocity = typename GET_PROP_TYPE(TypeTag, Velocity);
91 using CapillaryFlux = typename GET_PROP_TYPE(TypeTag, CapillaryFlux);
92 using GravityFlux = typename GET_PROP_TYPE(TypeTag, GravityFlux);
93
94 using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
95 using MaterialLaw = typename SpatialParams::MaterialLaw;
96
97 using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
98
99 using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
100 using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
101
102 using SolutionTypes = typename GET_PROP(TypeTag, SolutionTypes);
103
104 using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
105 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
106
107 using CellData = typename GET_PROP_TYPE(TypeTag, CellData);
108
109 enum
110 {
111 pw = Indices::pressureW,
112 pn = Indices::pressureNw,
113 pGlobal = Indices::pressureGlobal,
114 vw = Indices::velocityW,
115 vn = Indices::velocityNw,
116 vt = Indices::velocityTotal,
117 sw = Indices::saturationW,
118 sn = Indices::saturationNw
119 };
120 enum
121 {
122 wPhaseIdx = Indices::wPhaseIdx,
123 nPhaseIdx = Indices::nPhaseIdx,
124 saturationIdx = Indices::saturationIdx,
125 satEqIdx = Indices::satEqIdx,
126 numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
127 };
128
129 using TransportSolutionType = typename GET_PROP_TYPE(TypeTag, TransportSolutionType);
130
131 using Element = typename GridView::Traits::template Codim<0>::Entity;
132 using Intersection = typename GridView::Intersection;
133
134 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
135 using DimVector = Dune::FieldVector<Scalar, dim>;
136
137protected:
138 CapillaryFlux& capillaryFlux()
139 {
140 return *capillaryFlux_;
141 }
142
143 const CapillaryFlux& capillaryFlux() const
144 {
145 return *capillaryFlux_;
146 }
147
148 GravityFlux& gravityFlux()
149 {
150 return *gravityFlux_;
151 }
152
153 const GravityFlux& gravityFlux() const
154 {
155 return *gravityFlux_;
156 }
157
158public:
159 Velocity& velocity()
160 {
161 return *velocity_;
162 }
163
164 Velocity& velocity() const
165 {
166 return *velocity_;
167 }
168
170 void getFlux(Scalar& update, const Intersection& intersection, CellData& cellDataI);
171
173 void getFluxOnBoundary(Scalar& update, const Intersection& intersection, CellData& cellDataI);
174
176 void getSource(Scalar& update, const Element& element, CellData& cellDataI);
177
179 void initialize();
180
182 void updateMaterialLaws();
183
184
191 void getTransportedQuantity(TransportSolutionType& transportedQuantity)
192 {
193 int size = problem_.gridView().size(0);
194 transportedQuantity.resize(size);
195 for (int i = 0; i < size; i++)
196 {
197 switch (saturationType_)
198 {
199 case sw:
200 {
201 transportedQuantity[i] = problem_.variables().cellData(i).saturation(wPhaseIdx);
202 break;
203 }
204 case sn:
205 {
206 transportedQuantity[i] = problem_.variables().cellData(i).saturation(nPhaseIdx);
207 break;
208 }
209 }
210 }
211 }
212
219 void setTransportedQuantity(TransportSolutionType& transportedQuantity)
220 {
221 int size = problem_.gridView().size(0);
222 for (int i = 0; i < size; i++)
223 {
224 CellData& cellData = problem_.variables().cellData(i);
225 switch (saturationType_)
226 {
227 case sw:
228 {
229 Scalar sat = transportedQuantity[i];
230
231 cellData.setSaturation(wPhaseIdx, sat);
232 cellData.setSaturation(nPhaseIdx, 1 - sat);
233 break;
234 }
235 case sn:
236 {
237 Scalar sat = transportedQuantity[i];
238
239 cellData.setSaturation(wPhaseIdx,1 -sat);
240 cellData.setSaturation(nPhaseIdx, sat);
241 break;
242 }
243 }
244 }
245 }
246
253 template<class DataEntry>
254 bool inPhysicalRange(DataEntry& entry)
255 {
256 return (entry > -1e-6 && entry < 1.0 + 1e-6);
257 }
258
264 void updateTransportedQuantity(TransportSolutionType& updateVec)
265 {
266 if (this->enableLocalTimeStepping())
267 this->innerUpdate(updateVec);
268 else
269 asImp_().updateSaturationSolution(updateVec);
270 }
271
278 void updateTransportedQuantity(TransportSolutionType& updateVec, Scalar dt)
279 {
280 asImp_().updateSaturationSolution(updateVec, dt);
281 }
282
288 void updateSaturationSolution(TransportSolutionType& updateVec)
289 {
290 Scalar dt = problem_.timeManager().timeStepSize();
291 int size = problem_.gridView().size(0);
292 for (int i = 0; i < size; i++)
293 {
294 asImp_().updateSaturationSolution(i, updateVec[i][0], dt);
295 }
296 }
297
304 void updateSaturationSolution(TransportSolutionType& updateVec, Scalar dt)
305 {
306 int size = problem_.gridView().size(0);
307 for (int i = 0; i < size; i++)
308 {
309 asImp_().updateSaturationSolution(i, updateVec[i][0], dt);
310 }
311 }
312
322 void updateSaturationSolution(int eIdxGlobal, Scalar update, Scalar dt)
323 {
324 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
325
326 switch (saturationType_)
327 {
328 case sw:
329 {
330 Scalar sat = cellData.saturation(wPhaseIdx) + dt*update;
331
332 cellData.setSaturation(wPhaseIdx, sat);
333 cellData.setSaturation(nPhaseIdx, 1 - sat);
334 break;
335 }
336 case sn:
337 {
338 Scalar sat = cellData.saturation(nPhaseIdx) + dt*update;
339
340 cellData.setSaturation(wPhaseIdx,1 -sat);
341 cellData.setSaturation(nPhaseIdx, sat);
342 break;
343 }
344 }
345 }
346
359 template<class MultiWriter>
360 void addOutputVtkFields(MultiWriter &writer)
361 {
362 int size = problem_.gridView().size(0);
363 TransportSolutionType *saturation = writer.allocateManagedBuffer(size);
364 TransportSolutionType *saturationSecond = 0;
365
366 if (vtkOutputLevel_ > 0)
367 {
368 saturationSecond = writer.allocateManagedBuffer(size);
369 }
370
371 for (int i = 0; i < size; i++)
372 {
373 CellData& cellData = problem_.variables().cellData(i);
374
375 if (saturationType_ == sw)
376 {
377 (*saturation)[i] = cellData.saturation(wPhaseIdx);
378 if (vtkOutputLevel_ > 0)
379 {
380 (*saturationSecond)[i] = cellData.saturation(nPhaseIdx);
381 }
382 }
383 else if (saturationType_ == sn)
384 {
385 (*saturation)[i] = cellData.saturation(nPhaseIdx);
386 if (vtkOutputLevel_ > 0)
387 {
388 (*saturationSecond)[i] = cellData.saturation(wPhaseIdx);
389 }
390 }
391 }
392
393 if (saturationType_ == sw)
394 {
395 writer.attachCellData(*saturation, "wetting saturation");
396 if (vtkOutputLevel_ > 0)
397 {
398 writer.attachCellData(*saturationSecond, "nonwetting saturation");
399 }
400 }
401 else if (saturationType_ == sn)
402 {
403 writer.attachCellData(*saturation, "nonwetting saturation");
404 if (vtkOutputLevel_ > 0)
405 {
406 writer.attachCellData(*saturationSecond, "wetting saturation");
407 }
408 }
409
410 if (velocity().calculateVelocityInTransport())
411 velocity().addOutputVtkFields(writer);
412
413 return;
414 }
415
421 void serializeEntity(std::ostream &outstream, const Element &element)
422 {
423 int eIdxGlobal = problem_.variables().index(element);
424
425 Scalar sat = 0.0;
426 switch (saturationType_)
427 {
428 case sw:
429 sat = problem_.variables().cellData(eIdxGlobal).saturation(wPhaseIdx);
430 break;
431 case sn:
432 sat = problem_.variables().cellData(eIdxGlobal).saturation(nPhaseIdx);
433 break;
434 }
435
436 outstream << sat;
437 }
438
444 void deserializeEntity(std::istream &instream, const Element &element)
445 {
446 int eIdxGlobal = problem_.variables().index(element);
447
448 Scalar sat = 0.;
449 instream >> sat;
450
451 switch (saturationType_)
452 {
453 case sw:
454 problem_.variables().cellData(eIdxGlobal).setSaturation(wPhaseIdx, sat);
455 problem_.variables().cellData(eIdxGlobal).setSaturation(nPhaseIdx, 1-sat);
456 break;
457 case sn:
458 problem_.variables().cellData(eIdxGlobal).setSaturation(nPhaseIdx, sat);
459 problem_.variables().cellData(eIdxGlobal).setSaturation(wPhaseIdx, 1-sat);
460 break;
461 }
462 }
463
464
470 FVSaturation2P(Problem& problem) :
471 ParentType(problem), problem_(problem), threshold_(1e-6),
472 switchNormals_(getParam<bool>("Impet.SwitchNormals"))
473 {
474 if (compressibility_ && velocityType_ == vt)
475 {
476 DUNE_THROW(Dune::NotImplemented,
477 "Total velocity - global pressure - model cannot be used with compressible fluids!");
478 }
479 if (saturationType_ != sw && saturationType_ != sn)
480 {
481 DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
482 }
483 if (pressureType_ != pw && pressureType_ != pn && pressureType_ != pGlobal)
484 {
485 DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
486 }
487 if (velocityType_ != vw && velocityType_ != vn && velocityType_ != vt)
488 {
489 DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
490 }
491
492 capillaryFlux_ = std::make_shared<CapillaryFlux>(problem);
493 gravityFlux_ = std::make_shared<GravityFlux>(problem);
494 velocity_ = std::make_shared<Velocity>(problem);
495
496 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
497 porosityThreshold_ = getParam<Scalar>("Impet.PorosityThreshold");
498 }
499
500private:
502 Implementation &asImp_()
503 { return *static_cast<Implementation *>(this); }
504
506 const Implementation &asImp_() const
507 { return *static_cast<const Implementation *>(this); }
508
509 Problem& problem_;
510 std::shared_ptr<Velocity> velocity_;
511 std::shared_ptr<CapillaryFlux> capillaryFlux_;
512 std::shared_ptr<GravityFlux> gravityFlux_;
513
514 int vtkOutputLevel_;
515 Scalar porosityThreshold_;
516
517 static const bool compressibility_ = GET_PROP_VALUE(TypeTag, EnableCompressibility);
518 static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
519 static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation);
520 static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation);
521
522 Scalar density_[numPhases];
523 Scalar viscosity_[numPhases];
524
525 const Scalar threshold_;
526 const bool switchNormals_;
527};
528
540template<class TypeTag>
541void FVSaturation2P<TypeTag>::getFlux(Scalar& update, const Intersection& intersection, CellData& cellDataI)
542{
543 auto elementI = intersection.inside();
544 auto elementJ = intersection.outside();
545
546 const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(elementJ));
547
548 // get global coordinates of cell centers
549 const GlobalPosition& globalPosI = elementI.geometry().center();
550 const GlobalPosition& globalPosJ = elementJ.geometry().center();
551
552 // cell volume, assume linear map here
553 Scalar volume = elementI.geometry().volume();
554 using std::max;
555 Scalar porosity = max(problem_.spatialParams().porosity(elementI), porosityThreshold_);
556
557 if (compressibility_)
558 {
559 viscosity_[wPhaseIdx] = cellDataI.viscosity(wPhaseIdx);
560 viscosity_[nPhaseIdx] = cellDataI.viscosity(nPhaseIdx);
561 }
562
563 // local number of faces
564 int isIndex = intersection.indexInInside();
565
566 GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();
567 if (switchNormals_)
568 unitOuterNormal *= -1.0;
569
570 Scalar faceArea = intersection.geometry().volume();
571
572 if (velocity().calculateVelocityInTransport() && !cellDataI.fluxData().haveVelocity(isIndex))
573 velocity().calculateVelocity(intersection, cellDataI);
574
575 //get velocity*normalvector*facearea/(volume*porosity)
576 Scalar factorW = (cellDataI.fluxData().velocity(wPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
577 Scalar factorNw = (cellDataI.fluxData().velocity(nPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
578 Scalar factorTotal = factorW + factorNw;
579
580 // distance vector between barycenters
581 GlobalPosition distVec = globalPosJ - globalPosI;
582 // compute distance between cell centers
583 Scalar dist = distVec.two_norm();
584
585 bool takeNeighbor = (elementI.level() < elementJ.level());
586 //get phase potentials
587 bool upwindWI =
588 (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(wPhaseIdx, intersection.indexInOutside()) :
589 cellDataI.fluxData().isUpwindCell(wPhaseIdx, isIndex);
590 bool upwindNwI =
591 (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(nPhaseIdx, intersection.indexInOutside()) :
592 cellDataI.fluxData().isUpwindCell(nPhaseIdx, isIndex);
593
594 Scalar lambdaW = 0;
595 Scalar lambdaNw = 0;
596
597 //upwinding of lambda dependend on the phase potential gradients
598 if (upwindWI)
599 {
600 lambdaW = cellDataI.mobility(wPhaseIdx);
601 if (compressibility_)
602 {
603 lambdaW /= cellDataI.density(wPhaseIdx);
604 } //divide by density because lambda is saved as lambda*density
605 }
606 else
607 {
608 lambdaW = cellDataJ.mobility(wPhaseIdx);
609 if (compressibility_)
610 {
611 lambdaW /= cellDataJ.density(wPhaseIdx);
612 } //divide by density because lambda is saved as lambda*density
613 }
614
615 if (upwindNwI)
616 {
617 lambdaNw = cellDataI.mobility(nPhaseIdx);
618 if (compressibility_)
619 {
620 lambdaNw /= cellDataI.density(nPhaseIdx);
621 } //divide by density because lambda is saved as lambda*density
622 }
623 else
624 {
625 lambdaNw = cellDataJ.mobility(nPhaseIdx);
626 if (compressibility_)
627 {
628 lambdaNw /= cellDataJ.density(nPhaseIdx);
629 } //divide by density because lambda is saved as lambda*density
630 }
631
632 switch (velocityType_)
633 {
634 case vt:
635 {
636 //add cflFlux for time-stepping
637 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
638 viscosity_[nPhaseIdx], factorTotal, intersection);
639
640 //determine phase saturations from primary saturation variable
641 Scalar satWI = cellDataI.saturation(wPhaseIdx);
642 Scalar satWJ = cellDataJ.saturation(wPhaseIdx);
643
644 Scalar pcI = cellDataI.capillaryPressure();
645 Scalar pcJ = cellDataJ.capillaryPressure();
646
647 // calculate the saturation gradient
648 GlobalPosition pcGradient = unitOuterNormal;
649 pcGradient *= (pcI - pcJ) / dist;
650
651 // get the diffusive part
652 DimVector flux(0.);
653 capillaryFlux().getFlux(flux, intersection, satWI, satWJ, pcGradient);
654 Scalar capillaryFlux = (flux * unitOuterNormal * faceArea);
655
656 flux = 0.0;
657 gravityFlux().getFlux(flux, intersection, satWI, satWJ);
658 Scalar gravityFlux = (flux * unitOuterNormal * faceArea);
659
660 switch (saturationType_)
661 {
662 case sw:
663 {
664 //vt*fw
665 factorTotal *= lambdaW / (lambdaW + lambdaNw);
666 break;
667 }
668 case sn:
669 {
670 //vt*fn
671 factorTotal *= lambdaNw / (lambdaW + lambdaNw);
672 capillaryFlux *= -1;
673 gravityFlux *= -1;
674 break;
675 }
676 }
677 factorTotal -= capillaryFlux;
678 factorTotal += gravityFlux;
679
680 //add cflFlux for time-stepping
681 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
682 viscosity_[nPhaseIdx], 10 * capillaryFlux, intersection);
683 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
684 viscosity_[nPhaseIdx], 10 * gravityFlux, intersection);
685
686 break;
687 }
688 default:
689 {
690 if (compressibility_)
691 {
692 factorW /= cellDataI.density(wPhaseIdx);
693 factorNw /= cellDataI.density(nPhaseIdx);
694 }
695
696 //add cflFlux for time-stepping
697 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
698 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
699 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
700 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
701
702 break;
703 }
704 }
705
706 switch (velocityType_)
707 {
708 case vt:
709 update -= factorTotal / (volume * porosity); //-:v>0, if flow leaves the cell
710 break;
711 default:
712 switch (saturationType_)
713 {
714 case sw:
715 update -= factorW / (volume * porosity);//-:v>0, if flow leaves the cell
716 break;
717 case sn:
718 update -= factorNw / (volume * porosity); //-:v>0, if flow leaves the cell
719 break;
720 }
721 break;
722 }
723}
724
733template<class TypeTag>
734void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersection& intersection, CellData& cellDataI)
735{
736 auto elementI = intersection.inside();
737
738 // get global coordinates of cell centers
739 const GlobalPosition& globalPosI = elementI.geometry().center();
740
741 // center of face in global coordinates
742 const GlobalPosition& globalPosJ = intersection.geometry().center();
743
744 // cell volume, assume linear map here
745 Scalar volume = elementI.geometry().volume();
746 using std::max;
747 Scalar porosity = max(problem_.spatialParams().porosity(elementI), porosityThreshold_);
748
749 if (compressibility_)
750 {
751 viscosity_[wPhaseIdx] = cellDataI.viscosity(wPhaseIdx);
752 viscosity_[nPhaseIdx] = cellDataI.viscosity(nPhaseIdx);
753 }
754
755 // local number of faces
756 int isIndex = intersection.indexInInside();
757
758 GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();
759 if (switchNormals_)
760 unitOuterNormal *= -1.0;
761
762 Scalar faceArea = intersection.geometry().volume();
763
764 if (velocity().calculateVelocityInTransport())
765 velocity().calculateVelocityOnBoundary(intersection, cellDataI);
766
767 //get velocity*normalvector*facearea/(volume*porosity)
768 Scalar factorW = (cellDataI.fluxData().velocity(wPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
769 Scalar factorNw = (cellDataI.fluxData().velocity(nPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
770 Scalar factorTotal = factorW + factorNw;
771
772 // distance vector between barycenters
773 GlobalPosition distVec = globalPosJ - globalPosI;
774
775 // compute distance between cell centers
776 Scalar dist = distVec.two_norm();
777
778 //get boundary type
779 BoundaryTypes bcType;
780 problem_.boundaryTypes(bcType, intersection);
781 PrimaryVariables boundValues(0.0);
782
783 if (bcType.isDirichlet(satEqIdx))
784 {
785 problem_.dirichlet(boundValues, intersection);
786
787 Scalar satBound = boundValues[saturationIdx];
788
789 //determine phase saturations from primary saturation variable
790 Scalar satWI = cellDataI.saturation(wPhaseIdx);
791 Scalar satWBound = 0;
792 switch (saturationType_)
793 {
794 case sw:
795 {
796 satWBound = satBound;
797 break;
798 }
799 case sn:
800 {
801 satWBound = 1 - satBound;
802 break;
803 }
804 }
805
806 Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(elementI), satWBound);
807
808 Scalar lambdaW = 0;
809 Scalar lambdaNw = 0;
810
811 //upwinding of lambda dependend on the phase potential gradients
812 if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, isIndex))
813 {
814 lambdaW = cellDataI.mobility(wPhaseIdx);
815 if (compressibility_)
816 {
817 lambdaW /= cellDataI.density(wPhaseIdx);
818 } //divide by density because lambda is saved as lambda*density
819 }
820 else
821 {
822 if (compressibility_)
823 {
824 lambdaW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(elementI), satWBound)
825 / FluidSystem::viscosity(cellDataI.fluidState(), wPhaseIdx);
826 }
827 else
828 {
829 lambdaW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(elementI), satWBound)
830 / viscosity_[wPhaseIdx];
831 }
832 }
833
834 if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, isIndex))
835 {
836 lambdaNw = cellDataI.mobility(nPhaseIdx);
837 if (compressibility_)
838 {
839 lambdaNw /= cellDataI.density(nPhaseIdx);
840 } //divide by density because lambda is saved as lambda*density
841 }
842 else
843 {
844 if (compressibility_)
845 {
846 lambdaNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(elementI), satWBound)
847 / FluidSystem::viscosity(cellDataI.fluidState(), nPhaseIdx);
848 }
849 else
850 {
851 lambdaNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(elementI), satWBound)
852 / viscosity_[nPhaseIdx];
853 }
854 }
855
856 switch (velocityType_)
857 {
858 case vt:
859 {
860 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
861 viscosity_[nPhaseIdx], factorTotal, intersection);
862
863 Scalar pcI = cellDataI.capillaryPressure();
864
865 // calculate the saturation gradient
866 GlobalPosition pcGradient = unitOuterNormal;
867 pcGradient *= (pcI - pcBound) / dist;
868
869 // get the diffusive part -> give 1-sat because sat = S_n and lambda = lambda(S_w) and pc = pc(S_w)
870 DimVector flux(0.);
871 capillaryFlux().getFlux(flux, intersection, satWI, satWBound, pcGradient);
872 Scalar capillaryFlux = flux * unitOuterNormal * faceArea;
873
874 flux = 0.0;
875 gravityFlux().getFlux(flux, intersection, satWI, satWBound);
876 Scalar gravityFlux = flux * unitOuterNormal * faceArea;
877
878 switch (saturationType_)
879 {
880 case sw:
881 {
882 //vt*fw
883 factorTotal *= lambdaW / (lambdaW + lambdaNw);
884 break;
885 }
886 case sn:
887 {
888 //vt*fn
889 factorTotal *= lambdaNw / (lambdaW + lambdaNw);
890 capillaryFlux *= -1; //add cflFlux for time-stepping
891 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
892 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
893 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
894 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
895 gravityFlux *= -1;
896 break;
897 }
898 }
899 //vt*fw
900 factorTotal -= capillaryFlux;
901 factorTotal += gravityFlux;
902
903 //add cflFlux for time-stepping
904 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
905 viscosity_[nPhaseIdx], 10 * capillaryFlux, intersection);
906 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
907 viscosity_[nPhaseIdx], 10 * gravityFlux, intersection);
908
909 break;
910 }
911 default:
912 {
913 if (compressibility_)
914 {
915 factorW /= cellDataI.density(wPhaseIdx);
916 factorNw /= cellDataI.density(nPhaseIdx);
917 }
918
919 //add cflFlux for time-stepping
920 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
921 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
922 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
923 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
924
925 break;
926 }
927 }
928 } //end dirichlet boundary
929
930 if (bcType.isNeumann(satEqIdx))
931 {
932 problem_.neumann(boundValues, intersection);
933 factorW = boundValues[wPhaseIdx];
934 factorNw = boundValues[nPhaseIdx];
935 factorW *= faceArea;
936 factorNw *= faceArea;
937
938 //get mobilities
939 Scalar lambdaW, lambdaNw;
940
941 lambdaW = cellDataI.mobility(wPhaseIdx);
942 lambdaNw = cellDataI.mobility(nPhaseIdx);
943 if (compressibility_)
944 {
945 lambdaW /= cellDataI.density(wPhaseIdx);
946 lambdaNw /= cellDataI.density(nPhaseIdx);
947 factorW /= cellDataI.density(wPhaseIdx);
948 factorNw /= cellDataI.density(nPhaseIdx);
949 }
950 else
951 {
952 factorW /= density_[wPhaseIdx];
953 factorNw /= density_[nPhaseIdx];
954 }
955
956 switch (velocityType_)
957 {
958 case vt:
959 {
960 //add cflFlux for time-stepping
961 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
962 viscosity_[nPhaseIdx], factorW + factorNw, intersection);
963 break;
964 }
965 default:
966 {
967 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
968 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
969 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
970 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
971
972 break;
973 }
974 }
975
976 } //end neumann boundary
977 if (bcType.isOutflow(satEqIdx))
978 {
979 //get mobilities
980 Scalar lambdaW = cellDataI.mobility(wPhaseIdx);
981 Scalar lambdaNw = cellDataI.mobility(nPhaseIdx);
982 if (compressibility_)
983 {
984 lambdaW /= cellDataI.density(wPhaseIdx);
985 lambdaNw /= cellDataI.density(nPhaseIdx);
986 }
987
988 if (velocityType_ == vt)
989 {
990 switch (saturationType_)
991 {
992 case sw:
993 {
994 //vt*fw
995 factorTotal *= lambdaW / (lambdaW + lambdaNw);
996 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
997 viscosity_[nPhaseIdx], factorTotal, intersection);
998 break;
999 }
1000 case sn:
1001 {
1002 //vt*fn
1003 factorTotal *= lambdaNw / (lambdaW + lambdaNw);
1004 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1005 viscosity_[nPhaseIdx], factorTotal, intersection);
1006 break;
1007 }
1008 }
1009 }
1010 else
1011 {
1012 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1013 viscosity_[nPhaseIdx], factorW, intersection, wPhaseIdx);
1014 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1015 viscosity_[nPhaseIdx], factorNw, intersection, nPhaseIdx);
1016 }
1017 }
1018 switch (velocityType_)
1019 {
1020 case vt:
1021 update -= factorTotal / (volume * porosity); //-:v>0, if flow leaves the cell
1022 break;
1023 default:
1024 switch (saturationType_)
1025 {
1026 case sw:
1027 update -= factorW / (volume * porosity); //-:v>0, if flow leaves the cell
1028 break;
1029 case sn:
1030 update -= factorNw / (volume * porosity); //-:v>0, if flow leaves the cell
1031 break;
1032 }
1033 break;
1034 }
1035}
1036
1044template<class TypeTag>
1045void FVSaturation2P<TypeTag>::getSource(Scalar& update, const Element& element, CellData& cellDataI)
1046{
1047 // cell volume, assume linear map here
1048 Scalar volume = element.geometry().volume();
1049
1050 using std::max;
1051 Scalar porosity = max(problem_.spatialParams().porosity(element), porosityThreshold_);
1052
1053 if (compressibility_)
1054 {
1055 viscosity_[wPhaseIdx] = cellDataI.viscosity(wPhaseIdx);
1056 viscosity_[nPhaseIdx] = cellDataI.viscosity(nPhaseIdx);
1057 }
1058
1059 PrimaryVariables sourceVec(0.0);
1060 problem_.source(sourceVec, element);
1061
1062 if (compressibility_)
1063 {
1064 sourceVec[wPhaseIdx] /= cellDataI.density(wPhaseIdx);
1065 sourceVec[nPhaseIdx] /= cellDataI.density(nPhaseIdx);
1066 }
1067 else
1068 {
1069 sourceVec[wPhaseIdx] /= density_[wPhaseIdx];
1070 sourceVec[nPhaseIdx] /= density_[nPhaseIdx];
1071 }
1072
1073 //get mobilities
1074 Scalar lambdaW = cellDataI.mobility(wPhaseIdx);
1075 Scalar lambdaNw = cellDataI.mobility(nPhaseIdx);
1076 if (compressibility_)
1077 {
1078 lambdaW /= cellDataI.density(wPhaseIdx);
1079 lambdaNw /= cellDataI.density(nPhaseIdx);
1080 }
1081
1082 switch (saturationType_)
1083 {
1084 case sw:
1085 {
1086 if (sourceVec[wPhaseIdx] < 0 && cellDataI.saturation(wPhaseIdx) < threshold_)
1087 sourceVec[wPhaseIdx] = 0.0;
1088
1089 update += sourceVec[wPhaseIdx] / porosity;
1090 break;
1091 }
1092 case sn:
1093 {
1094 if (sourceVec[nPhaseIdx] < 0 && cellDataI.saturation(nPhaseIdx) < threshold_)
1095 sourceVec[nPhaseIdx] = 0.0;
1096
1097 update += sourceVec[nPhaseIdx] / porosity;
1098 break;
1099 }
1100 }
1101
1102 switch (velocityType_)
1103 {
1104 case vt:
1105 {
1106 //add cflFlux for time-stepping
1107 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx], viscosity_[nPhaseIdx],
1108 (sourceVec[wPhaseIdx] + sourceVec[nPhaseIdx]) * -1 * volume, element);
1109 break;
1110 }
1111 default:
1112 {
1113 //add cflFlux for time-stepping
1114 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1115 viscosity_[nPhaseIdx], sourceVec[wPhaseIdx] * -1 * volume, element,
1116 wPhaseIdx);
1117 this->evalCflFluxFunction().addFlux(lambdaW, lambdaNw, viscosity_[wPhaseIdx],
1118 viscosity_[nPhaseIdx], sourceVec[nPhaseIdx] * -1 * volume, element,
1119 nPhaseIdx);
1120 break;
1121 }
1122 }
1123}
1124
1126template<class TypeTag>
1128{
1129 ParentType::initialize();
1130
1131 if (!compressibility_)
1132 {
1133 const auto element = *problem_.gridView().template begin<0>();
1134 FluidState fluidState;
1135 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
1136 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
1137 fluidState.setTemperature(problem_.temperature(element));
1138 fluidState.setSaturation(wPhaseIdx, 1.);
1139 fluidState.setSaturation(nPhaseIdx, 0.);
1140 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
1141 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
1142 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
1143 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
1144 }
1145
1146 // iterate through leaf grid an evaluate c0 at cell center
1147 for (const auto& element : elements(problem_.gridView()))
1148 {
1149 PrimaryVariables initSol(0.0);
1150 problem_.initial(initSol, element);
1151
1152 int eIdxGlobal = problem_.variables().index(element);
1153
1154 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
1155
1156 switch (saturationType_)
1157 {
1158 case sw:
1159 {
1160 cellData.setSaturation(wPhaseIdx, initSol[saturationIdx]);
1161 cellData.setSaturation(nPhaseIdx, 1 - initSol[saturationIdx]);
1162 break;
1163 }
1164 case sn:
1165 {
1166 cellData.setSaturation(wPhaseIdx,1 -initSol[saturationIdx]);
1167 cellData.setSaturation(nPhaseIdx, initSol[saturationIdx]);
1168
1169 break;
1170 }
1171 }
1172 }
1173
1174 velocity_->initialize();
1175 capillaryFlux_->initialize();
1176 gravityFlux_->initialize();
1177
1178 return;
1179}
1180
1186template<class TypeTag>
1188{
1189 // iterate through leaf grid an evaluate c0 at cell center
1190 for (const auto& element : elements(problem_.gridView()))
1191 {
1192 int eIdxGlobal = problem_.variables().index(element);
1193
1194 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
1195
1196 //determine phase saturations from primary saturation variable
1197 Scalar satW = cellData.saturation(wPhaseIdx);
1198
1199 Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
1200
1201 cellData.setCapillaryPressure(pc);
1202
1203 // initialize mobilities
1204 Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[wPhaseIdx];
1205 Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[nPhaseIdx];
1206
1207 // initialize mobilities
1208 cellData.setMobility(wPhaseIdx, mobilityW);
1209 cellData.setMobility(nPhaseIdx, mobilityNw);
1210
1211 //initialize fractional flow functions
1212 cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
1213 cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
1214 }
1215 return;
1216}
1217
1218} // end namespace Dumux
1219#endif
#define GET_PROP_VALUE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:282
#define GET_PROP(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:281
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
Finite Volume discretization of a transport equation.
T getParam(Args &&... args)
A free function to get a parameter from the parameter tree singleton.
Definition: parameters.hh:428
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag EnableCompressibility
Returns whether compressibility is allowed.
Definition: porousmediumflow/2p/sequential/properties.hh:55
Property tag TransportModel
The type of the discretization of a transport model.
Definition: porousmediumflow/sequential/properties.hh:66
Property tag SaturationFormulation
The formulation of the saturation model.
Definition: porousmediumflow/2p/sequential/properties.hh:53
Property tag NumPhases
Number of phases in the system.
Definition: porousmediumflow/sequential/properties.hh:69
Property tag VelocityFormulation
The type of velocity reconstructed for the transport model.
Definition: porousmediumflow/2p/sequential/properties.hh:54
std::string saturation(int phaseIdx) noexcept
I/O name of saturation for multiphase systems.
Definition: name.hh:43
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
std::string porosity() noexcept
I/O name of porosity.
Definition: name.hh:139
The finite volume discretization of a saturation transport equation.
Definition: saturation.hh:76
const CapillaryFlux & capillaryFlux() const
Definition: saturation.hh:143
void getFlux(Scalar &update, const Intersection &intersection, CellData &cellDataI)
Function which calculates the flux update.
Definition: saturation.hh:541
Velocity & velocity()
Definition: saturation.hh:159
GravityFlux & gravityFlux()
Definition: saturation.hh:148
void updateTransportedQuantity(TransportSolutionType &updateVec)
Updates the primary transport variable.
Definition: saturation.hh:264
void updateSaturationSolution(int eIdxGlobal, Scalar update, Scalar dt)
Updates the saturation solution of a cell.
Definition: saturation.hh:322
void getFluxOnBoundary(Scalar &update, const Intersection &intersection, CellData &cellDataI)
Function which calculates the boundary flux update.
Definition: saturation.hh:734
void updateSaturationSolution(TransportSolutionType &updateVec)
Globally updates the saturation solution.
Definition: saturation.hh:288
CapillaryFlux & capillaryFlux()
Definition: saturation.hh:138
Velocity & velocity() const
Definition: saturation.hh:164
void updateSaturationSolution(TransportSolutionType &updateVec, Scalar dt)
Globally updates the saturation solution.
Definition: saturation.hh:304
void initialize()
Sets the initial solution.
Definition: saturation.hh:1127
void getTransportedQuantity(TransportSolutionType &transportedQuantity)
Writes the current values of the primary transport variable into the transportedQuantity-vector (come...
Definition: saturation.hh:191
const GravityFlux & gravityFlux() const
Definition: saturation.hh:153
void serializeEntity(std::ostream &outstream, const Element &element)
Function for serialization of the primary transport variable.
Definition: saturation.hh:421
void deserializeEntity(std::istream &instream, const Element &element)
Function for deserialization of the primary transport variable.
Definition: saturation.hh:444
void updateMaterialLaws()
Update the values of the material laws and constitutive relations.
Definition: saturation.hh:1187
void addOutputVtkFields(MultiWriter &writer)
Adds saturation output to the output file.
Definition: saturation.hh:360
FVSaturation2P(Problem &problem)
Constructs a FVSaturation2P object.
Definition: saturation.hh:470
void getSource(Scalar &update, const Element &element, CellData &cellDataI)
Function which calculates the source update.
Definition: saturation.hh:1045
void setTransportedQuantity(TransportSolutionType &transportedQuantity)
Writes the current values of the primary transport variable into the variable container.
Definition: saturation.hh:219
void updateTransportedQuantity(TransportSolutionType &updateVec, Scalar dt)
Updates the primary transport variable.
Definition: saturation.hh:278
bool inPhysicalRange(DataEntry &entry)
Check if saturation is in physical range.
Definition: saturation.hh:254
The finite volume discretization of a transport equation.
Definition: transport.hh:52
bool enableLocalTimeStepping()
Definition: transport.hh:202
void update(const Scalar t, Scalar &dt, TransportSolutionType &updateVec, bool impet=false)
Calculate the update vector.
Definition: transport.hh:270
void innerUpdate(TransportSolutionType &updateVec)
Definition: transport.hh:564
Specifies the properties for immiscible 2p transport.