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