3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
evalcflfluxcoats.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_EVALCFLFLUX_COATS_HH
25#define DUMUX_EVALCFLFLUX_COATS_HH
26
27#include <dune/common/float_cmp.hh>
28#include <dumux/porousmediumflow/sequential/impetproperties.hh>
29#include "evalcflflux.hh"
30
31namespace Dumux {
38template<class TypeTag>
39class EvalCflFluxCoats: public EvalCflFlux<TypeTag>
40{
41private:
42 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
43 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
44 using Problem = GetPropType<TypeTag, Properties::Problem>;
45
46 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
47 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
48 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
49 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
50
51 using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
52 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
53 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
54
55 using CellData = GetPropType<TypeTag, Properties::CellData>;
56
57 enum
58 {
59 dim = GridView::dimension, dimWorld = GridView::dimensionworld
60 };
61 enum
62 {
63 wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
64 eqIdxPress = Indices::pressureEqIdx,
65 eqIdxSat = Indices::satEqIdx,
66 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
67 };
68
69 enum
70 {
71 pw = Indices::pressureW,
72 pn = Indices::pressureNw,
73 vt = Indices::velocityTotal,
74 sw = Indices::saturationW,
75 sn = Indices::saturationNw
76 };
77
78 using Element = typename GridView::Traits::template Codim<0>::Entity;
79 using Intersection = typename GridView::Intersection;
80
81 using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
82 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
83 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
84
85public:
87 void initialize()
88 {
89 const auto element = *problem_.gridView().template begin<0>();
90 FluidState fluidState;
91 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
92 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
93 fluidState.setTemperature(problem_.temperature(element));
94 fluidState.setSaturation(wPhaseIdx, 1.);
95 fluidState.setSaturation(nPhaseIdx, 0.);
96 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
97 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
98 }
99
105 void addFlux(Scalar& lambdaW, Scalar& lambdaNw, Scalar& viscosityW, Scalar& viscosityNw, Scalar flux,
106 const Element& element, int phaseIdx = -1)
107 {
108 addDefaultFlux(flux, phaseIdx);
109 }
110
116 void addFlux(Scalar& lambdaW, Scalar& lambdaNw, Scalar& viscosityW, Scalar& viscosityNw, Scalar flux,
117 const Intersection& intersection, int phaseIdx = -1)
118 {
119 addDefaultFlux(flux, phaseIdx);
120 addCoatsFlux(lambdaW, lambdaNw, viscosityW, viscosityNw, flux, intersection, phaseIdx);
121 }
122
128 Scalar getCflFluxFunction(const Element& element)
129 {
130 Scalar cflFluxDefault = getCflFluxFunctionDefault();
131
132 if (rejectForTimeStepping_)
133 return 0.99 / cflFluxDefault;
134
135 using std::isnan;
136 using std::isinf;
137 if (isnan(cflFluxFunctionCoatsOut_) || isinf(cflFluxFunctionCoatsOut_)){cflFluxFunctionCoatsOut_ = 0.0;}
138 if (isnan(cflFluxFunctionCoatsIn_) || isinf(cflFluxFunctionCoatsIn_)){cflFluxFunctionCoatsIn_ = 0.0;}
139
140 using std::max;
141 Scalar cflFluxFunctionCoats = max(cflFluxFunctionCoatsIn_, cflFluxFunctionCoatsOut_);
142
143 if (cflFluxFunctionCoats <= 0)
144 {
145 return 0.99 / cflFluxDefault;
146 }
147 else if (cflFluxDefault > cflFluxFunctionCoats)
148 {
149 return 0.99 / cflFluxDefault;
150 }
151 else
152 {
153 return 0.99 / cflFluxFunctionCoats;
154 }
155 }
156
162 Scalar getDt(const Element& element)
163 {
164 using std::max;
165 Scalar porosity = max(problem_.spatialParams().porosity(element), porosityThreshold_);
166 return (getCflFluxFunction(element) * porosity * element.geometry().volume());
167 }
168
170 void reset()
171 {
172 cflFluxFunctionCoatsIn_ = 0;
173 cflFluxFunctionCoatsOut_ = 0;
174 rejectForTimeStepping_ = false;
175 fluxWettingOut_ = 0;
176 fluxNonwettingOut_ = 0;
177 fluxIn_ = 0;
178 fluxOut_ = 0;
179 }
180
186 EvalCflFluxCoats(Problem& problem) :
187 problem_(problem), epsDerivative_(5e-3), threshold_(1e-8)
188 {
189 reset();
190 density_[wPhaseIdx] = 0.;
191 density_[nPhaseIdx] = 0.;
192
193 porosityThreshold_ = getParam<Scalar>("Impet.PorosityThreshold");
194 }
195
196private:
197 Scalar getCflFluxFunctionDefault()
198 {
199 using std::isnan;
200 using std::isinf;
201 using std::max;
202
203 if (isnan(fluxIn_) || isinf(fluxIn_))
204 {
205 fluxIn_ = 1e-100;
206 }
207
208 Scalar cFLFluxIn = fluxIn_;
209
210
211 Scalar cFLFluxOut = 0;
212
213 if (velocityType_ == vt)
214 {
215 if (isnan(fluxOut_) || isinf(fluxOut_))
216 {
217 fluxOut_ = 1e-100;
218 }
219
220 cFLFluxOut = fluxOut_;
221 }
222 else
223 {
224 if (isnan(fluxWettingOut_) || isinf(fluxWettingOut_))
225 {
226 fluxWettingOut_ = 1e-100;
227 }
228 if (isnan(fluxNonwettingOut_) || isinf(fluxNonwettingOut_))
229 {
230 fluxNonwettingOut_ = 1e-100;
231 }
232
233 cFLFluxOut = max(fluxWettingOut_, fluxNonwettingOut_);
234 }
235
236
237 //determine timestep
238 Scalar cFLFluxFunction = max(cFLFluxIn, cFLFluxOut);
239
240 return cFLFluxFunction;
241 }
242
243 void addDefaultFlux(Scalar flux,int phaseIdx);
244
245 void addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw, Scalar& viscosityW, Scalar& viscosityNw, Scalar flux,
246 const Intersection& intersection, int phaseIdx);
247
248 Problem& problem_;//problem data
249 Scalar cflFluxFunctionCoatsIn_;
250 Scalar cflFluxFunctionCoatsOut_;
251 Scalar fluxWettingOut_;
252 Scalar fluxNonwettingOut_;
253 Scalar fluxOut_;
254 Scalar fluxIn_;
255 bool rejectForTimeStepping_;
256 Scalar density_[numPhases];
257 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
258 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
259 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
260 const Scalar epsDerivative_;
261 const Scalar threshold_;
262 Scalar porosityThreshold_;
263};
264
265template<class TypeTag>
266void EvalCflFluxCoats<TypeTag>::addDefaultFlux(Scalar flux, int phaseIdx)
267{
268 switch (phaseIdx)
269 {
270 case wPhaseIdx:
271 {
272 //for time step criterion
273 if (flux >= 0)
274 {
275 fluxWettingOut_ += flux;
276 }
277 if (flux < 0)
278 {
279 fluxIn_ -= flux;
280 }
281
282 break;
283 }
284
285 //for time step criterion if the nonwetting phase velocity is used
286 case nPhaseIdx:
287 {
288 if (flux >= 0)
289 {
290 fluxNonwettingOut_ += flux;
291 }
292 if (flux < 0)
293 {
294 fluxIn_ -= flux;
295 }
296
297 break;
298 }
299 default:
300 {
301 if (flux >= 0)
302 {
303 fluxOut_ += flux;
304 }
305 if (flux < 0)
306 {
307 fluxIn_ -= flux;
308 }
309
310 break;
311 }
312 }
313}
314
320template<class TypeTag>
321void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
322 Scalar& viscosityW, Scalar& viscosityNw, Scalar flux,
323 const Intersection& intersection, int phaseIdx)
324{
325 if (rejectForTimeStepping_)
326 return;
327 else if (phaseIdx != wPhaseIdx)
328 return;
329
330 Scalar lambdaT = (lambdaW + lambdaNw);
331
332 if (lambdaT <= threshold_)
333 return;
334
335 auto element = intersection.inside();
336
337 //coordinates of cell center
338 const GlobalPosition& globalPos = element.geometry().center();
339
340 // cell index
341 int globalIdxI = problem_.variables().index(element);
342
343 CellData& cellDataI = problem_.variables().cellData(globalIdxI);
344
345 if (cellDataI.potential(wPhaseIdx) < 0.0 || cellDataI.potential(nPhaseIdx) < 0.0)
346 {
347 rejectForTimeStepping_ = true;
348 cflFluxFunctionCoatsIn_ = 0;
349 cflFluxFunctionCoatsOut_ = 0;
350 return;
351 }
352
353 int indexInInside = intersection.indexInInside();
354
355 Scalar satI = cellDataI.saturation(wPhaseIdx);
356 Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
357 Scalar lambdaNwI = cellDataI.mobility(nPhaseIdx);
358
359 const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(element.geometry().center());
360 const Scalar dpc_dsI = fluidMatrixInteraction.dpc_dsw(satI);
361
362 const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
363
364 if (intersection.neighbor())
365 {
366 auto neighbor = intersection.outside();
367
368 int globalIdxJ = problem_.variables().index(neighbor);
369
370 CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
371
372 if (cellDataJ.potential(wPhaseIdx) < 0.0 || cellDataJ.potential(nPhaseIdx) < 0.0 )
373 {
374 rejectForTimeStepping_ = true;
375 cflFluxFunctionCoatsIn_ = 0;
376 cflFluxFunctionCoatsOut_ = 0;
377 return;
378 }
379
380 if (element.level() != neighbor.level())
381 {
382 rejectForTimeStepping_ = true;
383 cflFluxFunctionCoatsIn_ = 0;
384 cflFluxFunctionCoatsOut_ = 0;
385 return;
386 }
387
388 bool takeNeighbor = (element.level() < neighbor.level());
389 //get phase potentials
390 bool upwindWI =
391 (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(wPhaseIdx, intersection.indexInOutside()) :
392 cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside);
393 bool upwindNwI =
394 (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(nPhaseIdx, intersection.indexInOutside()) :
395 cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside);
396
397 const GlobalPosition& globalPosNeighbor = neighbor.geometry().center();
398
399 // distance vector between barycenters
400 GlobalPosition distVec = globalPosNeighbor - globalPos;
401
402 // compute distance between cell centers
403 using std::abs;
404 Scalar dist = abs(distVec * unitOuterNormal);
405
406 Scalar satJ = cellDataJ.saturation(wPhaseIdx);
407 Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx);
408 Scalar lambdaNwJ = cellDataI.mobility(nPhaseIdx);
409
410 const auto fluidMatrixInteractionNeighbor = problem_.spatialParams().fluidMatrixInteractionAtPos(neighbor.geometry().center());
411 const Scalar dpc_dsJ = fluidMatrixInteractionNeighbor.dpc_dsw(satJ);
412
413 // compute vectorized permeabilities
414 DimVector permeability(0);
415 DimMatrix perm(0);
416 problem_.spatialParams().meanK(perm, problem_.spatialParams().intrinsicPermeability(element));
417 perm.mv(unitOuterNormal, permeability);
418
419 Scalar perm1 = permeability * unitOuterNormal;
420
421 permeability = 0;
422 problem_.spatialParams().meanK(perm, problem_.spatialParams().intrinsicPermeability(neighbor));
423 perm.mv(unitOuterNormal, permeability);
424
425 Scalar perm2 = permeability * unitOuterNormal;
426
427 Scalar meanPermeability = 2*perm1*perm2/(perm1 + perm2);
428
429 Scalar transmissibility = meanPermeability * intersection.geometry().volume() / dist;
430
431 Scalar satUpw = 0;
432 using std::max;
433 if (upwindWI)
434 {
435 satUpw = max(satI, 0.0);
436 }
437 else
438 {
439 satUpw = max(satJ, 0.0);
440 }
441
442 Scalar ds = epsDerivative_;
443
444 Scalar satPlus = satUpw + epsDerivative_;
445 Scalar satMinus = satUpw;
446 if (satMinus - epsDerivative_ > 0.0)
447 {
448 satMinus -= epsDerivative_;
449 ds += epsDerivative_;
450 }
451
452 Scalar dLambdaWDs = fluidMatrixInteractionNeighbor.krw(abs(satPlus)) / viscosityW;
453 dLambdaWDs -= fluidMatrixInteractionNeighbor.krw(abs(satMinus)) / viscosityW;
454 dLambdaWDs /= (ds);
455
456 if (upwindNwI)
457 {
458 satUpw = max(1 - satI, 0.0);
459 }
460 else
461 {
462 satUpw = max(1 - satJ, 0.0);
463 }
464
465 ds = epsDerivative_;
466
467 satPlus = satUpw + epsDerivative_;
468 satMinus = satUpw;
469 if (satMinus - epsDerivative_ > 0.0)
470 {
471 satMinus -= epsDerivative_;
472 ds += epsDerivative_;
473 }
474
475 Scalar dLambdaNwDs = fluidMatrixInteractionNeighbor.krn(satPlus) / viscosityNw;
476 dLambdaNwDs -= fluidMatrixInteractionNeighbor.krn(satMinus) / viscosityNw;
477 dLambdaNwDs /= (ds);
478
479 Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWJ);
480 Scalar lambdaNwCap = 0.5 * (lambdaNwI + lambdaNwJ);
481
482 Scalar potentialDiff = cellDataI.potential(wPhaseIdx) - cellDataJ.potential(wPhaseIdx);
483 Scalar cflFlux = transmissibility * lambdaNw * dLambdaWDs * abs(potentialDiff) / lambdaT;
484
485 potentialDiff = cellDataI.potential(nPhaseIdx) - cellDataJ.potential(nPhaseIdx);
486 cflFlux -= transmissibility * lambdaW * dLambdaNwDs * abs(potentialDiff) / lambdaT;
487
488 cflFlux -= transmissibility * lambdaWCap * lambdaNwCap * (dpc_dsI + dpc_dsJ) / lambdaT;
489
490 if ((upwindWI && lambdaW > threshold_)|| (upwindNwI && lambdaW < threshold_))
491 {
492 cflFluxFunctionCoatsOut_ += cflFlux;
493 }
494 else
495 {
496 cflFluxFunctionCoatsIn_ += cflFlux;
497 }
498 }
499 else
500 {
501 // center of face in global coordinates
502 GlobalPosition globalPosFace = intersection.geometry().center();
503
504 //get boundary type
505 BoundaryTypes bcType;
506 problem_.boundaryTypes(bcType, intersection);
507
508 // distance vector between barycenters
509 Dune::FieldVector < Scalar, dimWorld > distVec = globalPosFace - globalPos;
510
511 // compute distance between cell centers
512 Scalar dist = distVec.two_norm();
513
514 //permeability vector at boundary
515 // compute vectorized permeabilities
516
517 Dune::FieldVector<Scalar, dim> permeability(0);
518 DimMatrix perm(0);
519 problem_.spatialParams().meanK(perm, problem_.spatialParams().intrinsicPermeability(element));
520 perm.mv(unitOuterNormal, permeability);
521
522 Scalar faceArea = intersection.geometry().volume();
523
524 Scalar transmissibility = (unitOuterNormal * permeability) * faceArea / dist;
525
526 Scalar satWBound = cellDataI.saturation(wPhaseIdx);
527 if (bcType.isDirichlet(eqIdxSat))
528 {
529 PrimaryVariables bcValues;
530 problem_.dirichlet(bcValues, intersection);
531 switch (saturationType_)
532 {
533 case sw:
534 {
535 satWBound = bcValues[eqIdxSat];
536 break;
537 }
538 case sn:
539 {
540 satWBound = 1 - bcValues[eqIdxSat];
541 break;
542 }
543 default:
544 {
545 DUNE_THROW(Dune::RangeError, "saturation type not implemented");
546 }
547 }
548
549 }
550
551 Scalar potWBound = cellDataI.potential(wPhaseIdx);
552 Scalar potNwBound = cellDataI.potential(nPhaseIdx);
553 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * problem_.gravity();
554 if (bcType.isDirichlet(eqIdxPress))
555 {
556 PrimaryVariables bcValues;
557 problem_.dirichlet(bcValues, intersection);
558 switch (pressureType_)
559 {
560 case pw:
561 {
562 potWBound = bcValues[eqIdxPress] + density_[wPhaseIdx] * gdeltaZ;
563 potNwBound = bcValues[eqIdxPress] + fluidMatrixInteraction.pc(satWBound)
564 + density_[nPhaseIdx] * gdeltaZ;
565 break;
566 }
567 case pn:
568 {
569 potWBound = bcValues[eqIdxPress] - fluidMatrixInteraction.pc(satWBound)
570 + density_[wPhaseIdx] * gdeltaZ;
571 potNwBound = bcValues[eqIdxPress] + density_[nPhaseIdx] * gdeltaZ;
572 break;
573 }
574 default:
575 {
576 DUNE_THROW(Dune::RangeError, "pressure type not implemented");
577 }
578 }
579 }
580 else if (bcType.isNeumann(eqIdxPress) && bcType.isDirichlet(eqIdxSat))
581 {
582 PrimaryVariables bcValues;
583 problem_.neumann(bcValues, intersection);
584
585 bcValues[wPhaseIdx] /= density_[wPhaseIdx];
586 bcValues[nPhaseIdx] /= density_[nPhaseIdx];
587
588 bcValues[wPhaseIdx] *= faceArea;
589 bcValues[nPhaseIdx] *= faceArea;
590
591 bool hasPotWBound = false;
592 if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(lambdaW, 0.0, 1.0e-30) && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(bcValues[wPhaseIdx], 0.0, 1.0e-30))
593 {
594 potWBound -= bcValues[wPhaseIdx] / (transmissibility * lambdaW);
595 hasPotWBound = true;
596 }
597 bool hasPotNwBound = false;
598 if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(lambdaNw, 0.0, 1.0e-30) && Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(bcValues[nPhaseIdx], 0.0, 1.0e-30))
599 {
600 potNwBound -= bcValues[nPhaseIdx] / (transmissibility * lambdaNw);
601 hasPotNwBound = true;
602 }
603
604 if (hasPotWBound && !hasPotNwBound)
605 {
606 potNwBound = potWBound + fluidMatrixInteraction.pc(satWBound)
607 + (density_[nPhaseIdx] - density_[wPhaseIdx]) * gdeltaZ;
608 }
609 else if (!hasPotWBound && hasPotNwBound)
610 {
611 potWBound = potNwBound - fluidMatrixInteraction.pc(satWBound)
612 + (density_[nPhaseIdx] - density_[wPhaseIdx]) * gdeltaZ;
613 }
614 }
615 else if (bcType.isNeumann(eqIdxPress))
616 {
617 PrimaryVariables bcValues;
618 problem_.neumann(bcValues, intersection);
619
620 bcValues[wPhaseIdx] /= density_[wPhaseIdx];
621 bcValues[nPhaseIdx] /= density_[nPhaseIdx];
622
623 bcValues[wPhaseIdx] *= faceArea;
624 bcValues[nPhaseIdx] *= faceArea;
625
626 using std::abs;
627
628 if (bcValues[wPhaseIdx] > 0)
629 {
630 cflFluxFunctionCoatsOut_ += abs(bcValues[wPhaseIdx]);
631 }
632 else
633 {
634 cflFluxFunctionCoatsIn_ += abs(bcValues[wPhaseIdx]);
635 }
636 if (bcValues[nPhaseIdx] > 0)
637 {
638 cflFluxFunctionCoatsOut_ += abs(bcValues[nPhaseIdx]);
639 }
640 else
641 {
642 cflFluxFunctionCoatsIn_ += abs(bcValues[nPhaseIdx]);
643 }
644
645 return;
646 }
647 else
648 {
649 rejectForTimeStepping_ = true;
650 cflFluxFunctionCoatsIn_ = 0;
651 cflFluxFunctionCoatsOut_ = 0;
652 return;
653 }
654
655 const Scalar dpc_dsBound = fluidMatrixInteraction.dpc_dsw(satWBound);
656
657 Scalar lambdaWBound = 0;
658 Scalar lambdaNwBound = 0;
659
660 Scalar temperature = problem_.temperature(element);
661 Scalar referencePressure = problem_.referencePressure(element);
662 FluidState fluidState;
663 fluidState.setPressure(wPhaseIdx, referencePressure);
664 fluidState.setPressure(nPhaseIdx, referencePressure);
665 fluidState.setTemperature(temperature);
666
667 Scalar viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx);
668 Scalar viscosityNwBound =
669 FluidSystem::viscosity(fluidState, nPhaseIdx);
670 lambdaWBound = fluidMatrixInteraction.krw(satWBound) / viscosityWBound;
671 lambdaNwBound = fluidMatrixInteraction.krn(satWBound) / viscosityNwBound;
672
673 Scalar satUpw = 0;
674 using std::max;
675 if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside))
676 {
677 satUpw = max(satI, 0.0);
678 }
679 else
680 {
681 satUpw = max(satWBound, 0.0);
682 }
683
684 Scalar ds = epsDerivative_;
685
686 Scalar satPlus = satUpw + epsDerivative_;
687 Scalar satMinus = satUpw;
688 if (satMinus - epsDerivative_ > 0.0)
689 {
690 satMinus -= epsDerivative_;
691 ds += epsDerivative_;
692 }
693
694 Scalar dLambdaWDs = fluidMatrixInteraction.krw(satPlus) / viscosityW;
695 dLambdaWDs -= fluidMatrixInteraction.krw(satMinus) / viscosityW;
696 dLambdaWDs /= (ds);
697
698 if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
699 {
700 satUpw = max(1 - satI, 0.0);
701 }
702 else
703 {
704 satUpw = max(1 - satWBound, 0.0);
705 }
706
707 ds = epsDerivative_;
708
709 satPlus = satUpw + epsDerivative_;
710 satMinus = satUpw;
711 if (satMinus - epsDerivative_ > 0.0)
712 {
713 satMinus -= epsDerivative_;
714 ds += epsDerivative_;
715 }
716
717 Scalar dLambdaNwDs = fluidMatrixInteraction.krn(satPlus) / viscosityNw;
718 dLambdaNwDs -= fluidMatrixInteraction.krn(satMinus) / viscosityNw;
719 dLambdaNwDs /= (ds);
720
721 Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWBound);
722 Scalar lambdaNwCap = 0.5 * (lambdaNwI + lambdaNwBound);
723
724 Scalar potDiff = cellDataI.potential(wPhaseIdx) - potWBound;
725 using std::abs;
726 Scalar cflFlux = transmissibility * lambdaNw * dLambdaWDs * abs(potDiff) / lambdaT;
727
728 cflFlux -= transmissibility * lambdaWCap * lambdaNwCap * (dpc_dsI + dpc_dsBound) / lambdaT;
729
730 potDiff = cellDataI.potential(nPhaseIdx) - potNwBound;
731 cflFlux -= transmissibility * lambdaW * dLambdaNwDs * abs(potDiff) / lambdaT;
732
733 if ((cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside) && lambdaW > threshold_) ||
734 (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside) && lambdaW < threshold_))
735 {
736 cflFluxFunctionCoatsOut_ += cflFlux;
737 }
738 else
739 {
740 cflFluxFunctionCoatsIn_ += cflFlux;
741 }
742 }
743}
744
745} // end namespace Dumux
746
747#endif
evalcflflux.hh
Base class for implementations of different kinds of fluxes to evaluate a CFL-Condition.
impetproperties.hh
Base file for properties related to sequential IMPET algorithms.
Dumux
Definition: adapt.hh:29
Dumux::GetProp
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property
Definition: propertysystem.hh:141
Dumux::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
Dumux::IOName::temperature
std::string temperature() noexcept
I/O name of temperature for equilibrium models.
Definition: name.hh:51
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::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Dumux::IOName::porosity
std::string porosity() noexcept
I/O name of porosity.
Definition: name.hh:139
Dumux::EvalCflFlux
Base class for implementations of different kinds of fluxes to evaluate a CFL-Condition.
Definition: evalcflflux.hh:52
Dumux::EvalCflFluxCoats
Cfl-flux-function to evaluate a Cfl-Condition after Coats 2003.
Definition: evalcflfluxcoats.hh:40
Dumux::EvalCflFluxCoats::getDt
Scalar getDt(const Element &element)
Returns the Cfl time-step.
Definition: evalcflfluxcoats.hh:162
Dumux::EvalCflFluxCoats::reset
void reset()
Resets the Timestep-estimator.
Definition: evalcflfluxcoats.hh:170
Dumux::EvalCflFluxCoats::addFlux
void addFlux(Scalar &lambdaW, Scalar &lambdaNw, Scalar &viscosityW, Scalar &viscosityNw, Scalar flux, const Element &element, int phaseIdx=-1)
adds a flux to the cfl-criterion evaluation
Definition: evalcflfluxcoats.hh:105
Dumux::EvalCflFluxCoats::getCflFluxFunction
Scalar getCflFluxFunction(const Element &element)
Returns the Cfl flux-function.
Definition: evalcflfluxcoats.hh:128
Dumux::EvalCflFluxCoats::EvalCflFluxCoats
EvalCflFluxCoats(Problem &problem)
Constructs an EvalCflFluxDefault object.
Definition: evalcflfluxcoats.hh:186
Dumux::EvalCflFluxCoats::addFlux
void addFlux(Scalar &lambdaW, Scalar &lambdaNw, Scalar &viscosityW, Scalar &viscosityNw, Scalar flux, const Intersection &intersection, int phaseIdx=-1)
adds a flux to the cfl-criterion evaluation
Definition: evalcflfluxcoats.hh:116
Dumux::EvalCflFluxCoats::initialize
void initialize()
Initializes the cfl-flux-model.
Definition: evalcflfluxcoats.hh:87
Generated by doxygen 1.9.3