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