diff --git a/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh b/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
index 51348de680dd98e988ad4b78b58514f87dc29e25..62e6e81f64af6fcb1a8bb7131b30d0a2f037df84 100644
--- a/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
+++ b/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
@@ -127,21 +127,22 @@ public:
*/
Scalar getCFLFluxFunction(const Element& element)
{
- Scalar cflFluxDefault = getCFLFluxFunctionDefault();
+ if (rejectForTimeStepping_)
+ return 1e100;
- // std::cout<<"cflFluxFunctionCoats_ = "<<cflFluxFunctionCoats_<<", cflFluxDefault = "<<cflFluxDefault<<"\n";
+ Scalar cflFluxDefault = getCFLFluxFunctionDefault();
if (std::isnan(cflFluxFunctionCoats_) || std::isinf(cflFluxFunctionCoats_))
{
return 0.99 / cflFluxDefault;
}
- else if (hasHangingNode_)
+ else if (cflFluxFunctionCoats_ <= 0)
{
return 0.99 / cflFluxDefault;
}
- else if (cflFluxFunctionCoats_ <= 0)
+ else if (cflFluxDefault > cflFluxFunctionCoats_)
{
- return 0.99 / cflFluxDefault;
+ return 0.99 / cflFluxDefault;
}
else
{
@@ -155,18 +156,19 @@ public:
*/
Scalar getDt(const Element& element)
{
- Scalar porosity = problem_.spatialParams().porosity(element);
- if (porosity > 1e-6)
- return (getCFLFluxFunction(element) * problem_.spatialParams().porosity(element) * element.geometry().volume());
- else
- return (getCFLFluxFunction(element) * element.geometry().volume());
+ if (rejectForTimeStepping_)
+ return 1e100;
+
+ Scalar porosity = std::max(problem_.spatialParams().porosity(element), porosityThreshold_);
+
+ return (getCFLFluxFunction(element) * porosity * element.geometry().volume());
}
//! \brief Resets the Timestep-estimator
void reset()
{
cflFluxFunctionCoats_ = 0;
- hasHangingNode_ = false;
+ rejectForTimeStepping_ = false;
fluxWettingOut_ = 0;
fluxNonwettingOut_ = 0;
fluxIn_ = 0;
@@ -183,6 +185,8 @@ public:
reset();
density_[wPhaseIdx] = 0.;
density_[nPhaseIdx] = 0.;
+
+ porosityThreshold_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, PorosityThreshold);
}
private:
@@ -239,13 +243,14 @@ private:
Scalar fluxNonwettingOut_;
Scalar fluxOut_;
Scalar fluxIn_;
- bool hasHangingNode_;
+ bool rejectForTimeStepping_;
Scalar density_[numPhases];
static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation);
static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation);
static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
const Scalar epsDerivative_;
const Scalar threshold_;
+ Scalar porosityThreshold_;
};
template<class TypeTag>
@@ -306,12 +311,10 @@ template<class TypeTag>
void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux,
const Intersection& intersection, int phaseIdx)
{
- Scalar lambdaT = (lambdaW + lambdaNW);
+ if (rejectForTimeStepping_)
+ return;
- if (lambdaT <= 0.0)
- {
- return;
- }
+ Scalar lambdaT = (lambdaW + lambdaNW);
ElementPointer element = intersection.inside();
@@ -323,6 +326,13 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNW,
CellData& cellDataI = problem_.variables().cellData(globalIdxI);
+ if (cellDataI.pressure(wPhaseIdx) < 0.0 || cellDataI.pressure(nPhaseIdx) < 0.0 )
+ {
+ rejectForTimeStepping_ = true;
+ cflFluxFunctionCoats_ = 0;
+ return;
+ }
+
int indexInInside = intersection.indexInInside();
Scalar satI = cellDataI.saturation(wPhaseIdx);
@@ -343,21 +353,32 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNW,
Dune::FieldVector < Scalar, dimWorld > distVec = globalPosNeighbor - globalPos;
// compute distance between cell centers
- // Scalar dist = distVec.two_norm();
- Scalar dist = std::abs(distVec*unitOuterNormal);
+ Scalar dist = distVec.two_norm();
int globalIdxJ = problem_.variables().index(*neighbor);
CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
+ if (cellDataJ.pressure(wPhaseIdx) < 0.0 || cellDataJ.pressure(nPhaseIdx) < 0.0 )
+ {
+ rejectForTimeStepping_ = true;
+ cflFluxFunctionCoats_ = 0;
+ return;
+ }
+
//calculate potential gradients
if (element.level() != neighbor.level())
{
- hasHangingNode_ = true;
+ rejectForTimeStepping_ = true;
cflFluxFunctionCoats_ = 0;
return;
}
+ if (lambdaT <= 0.0)
+ {
+ return;
+ }
+
Scalar satJ = cellDataJ.saturation(wPhaseIdx);
Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx);
Scalar lambdaNWJ = cellDataI.mobility(nPhaseIdx);
@@ -460,6 +481,11 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNW,
}
else
{
+ if (lambdaT <= 0.0)
+ {
+ return;
+ }
+
// center of face in global coordinates
GlobalPosition globalPosFace = intersection.geometry().center();
diff --git a/dumux/decoupled/2p/transport/fv/evalcflfluxdefault.hh b/dumux/decoupled/2p/transport/fv/evalcflfluxdefault.hh
index c0cb35b6b8706d0763667b508518546942c42f1d..78acf754750d2ec629cade4a9601326604cee283 100644
--- a/dumux/decoupled/2p/transport/fv/evalcflfluxdefault.hh
+++ b/dumux/decoupled/2p/transport/fv/evalcflfluxdefault.hh
@@ -104,11 +104,10 @@ public:
*/
Scalar getDt(const Element& element)
{
- Scalar porosity = problem_.spatialParams().porosity(element);
- if (porosity > 1e-6)
- return (getCFLFluxFunction(element) * problem_.spatialParams().porosity(element) * element.geometry().volume());
- else
- return (getCFLFluxFunction(element) * element.geometry().volume());
+ Scalar porosity = std::max(problem_.spatialParams().porosity(element), porosityThreshold_);
+
+ return (getCFLFluxFunction(element) * problem_.spatialParams().porosity(element) * element.geometry().volume());
+
}
//! resets the accumulated CFL-fluxes to zero
diff --git a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
index c2dd7e25ee25eaf206aaabd33c91b886390f39f6..a6d501677f0d8f237f7a6950e65e30a9ff7e4d20 100644
--- a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
+++ b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
@@ -128,16 +128,6 @@ class FVSaturation2P: public FVTransport<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dim> DimVector;
- Velocity& velocity()
- {
- return *velocity_;
- }
-
- Velocity& velocity() const
- {
- return *velocity_;
- }
-
CapillaryFlux& capillaryFlux()
{
return *capillaryFlux_;
@@ -159,6 +149,16 @@ class FVSaturation2P: public FVTransport<TypeTag>
}
public:
+ Velocity& velocity()
+ {
+ return *velocity_;
+ }
+
+ Velocity& velocity() const
+ {
+ return *velocity_;
+ }
+
// Function which calculates the flux update
void getFlux(Scalar& update, const Intersection& intersection, CellData& cellDataI);
@@ -407,6 +407,7 @@ public:
velocity_ = Dune::make_shared<Velocity>(problem);
vtkOutputLevel_ = GET_PARAM_FROM_GROUP(TypeTag, int, Vtk, OutputLevel);
+ porosityThreshold_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, PorosityThreshold);
}
private:
@@ -416,6 +417,7 @@ private:
Dune::shared_ptr<GravityFlux> gravityFlux_;
int vtkOutputLevel_;
+ Scalar porosityThreshold_;
static const bool compressibility_ = GET_PROP_VALUE(TypeTag, EnableCompressibility);
static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
@@ -450,7 +452,7 @@ void FVSaturation2P<TypeTag>::getFlux(Scalar& update, const Intersection& inters
// cell volume, assume linear map here
Scalar volume = elementI->geometry().volume();
- Scalar porosity = problem_.spatialParams().porosity(*elementI);
+ Scalar porosity = std::max(problem_.spatialParams().porosity(*elementI), porosityThreshold_);
if (compressibility_)
{
@@ -599,16 +601,16 @@ void FVSaturation2P<TypeTag>::getFlux(Scalar& update, const Intersection& inters
switch (velocityType_)
{
case vt:
- if (porosity > 0.0){update -= factorTotal / (volume * porosity);} //-:v>0, if flow leaves the cell
+ update -= factorTotal / (volume * porosity); //-:v>0, if flow leaves the cell
break;
default:
switch (saturationType_)
{
case Sw:
- if (porosity > 0.0){ update -= factorW / (volume * porosity);}//-:v>0, if flow leaves the cell
+ update -= factorW / (volume * porosity);//-:v>0, if flow leaves the cell
break;
case Sn:
- if (porosity > 0.0){update -= factorNW / (volume * porosity);} //-:v>0, if flow leaves the cell
+ update -= factorNW / (volume * porosity); //-:v>0, if flow leaves the cell
break;
}
break;
@@ -635,7 +637,7 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti
// cell volume, assume linear map here
Scalar volume = elementI->geometry().volume();
- Scalar porosity = problem_.spatialParams().porosity(*elementI);
+ Scalar porosity = std::max(problem_.spatialParams().porosity(*elementI), porosityThreshold_);
if (compressibility_)
{
@@ -895,16 +897,16 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti
switch (velocityType_)
{
case vt:
- if (porosity > 0.0){ update -= factorTotal / (volume * porosity);} //-:v>0, if flow leaves the cell
+ update -= factorTotal / (volume * porosity); //-:v>0, if flow leaves the cell
break;
default:
switch (saturationType_)
{
case Sw:
- if (porosity > 0.0){ update -= factorW / (volume * porosity);} //-:v>0, if flow leaves the cell
+ update -= factorW / (volume * porosity); //-:v>0, if flow leaves the cell
break;
case Sn:
- if (porosity > 0.0){ update -= factorNW / (volume * porosity);} //-:v>0, if flow leaves the cell
+ update -= factorNW / (volume * porosity); //-:v>0, if flow leaves the cell
break;
}
break;
@@ -923,7 +925,7 @@ void FVSaturation2P<TypeTag>::getSource(Scalar& update, const Element& element,
// cell volume, assume linear map here
Scalar volume = element.geometry().volume();
- Scalar porosity = problem_.spatialParams().porosity(element);
+ Scalar porosity = std::max(problem_.spatialParams().porosity(element), porosityThreshold_);
if (compressibility_)
{
@@ -961,7 +963,7 @@ void FVSaturation2P<TypeTag>::getSource(Scalar& update, const Element& element,
if (sourceVec[wPhaseIdx] < 0 && cellDataI.saturation(wPhaseIdx) < threshold_)
sourceVec[wPhaseIdx] = 0.0;
- if (porosity > 0.0){update += sourceVec[wPhaseIdx] / porosity;}
+ update += sourceVec[wPhaseIdx] / porosity;
break;
}
case Sn:
@@ -969,7 +971,7 @@ void FVSaturation2P<TypeTag>::getSource(Scalar& update, const Element& element,
if (sourceVec[nPhaseIdx] < 0 && cellDataI.saturation(nPhaseIdx) < threshold_)
sourceVec[nPhaseIdx] = 0.0;
- if (porosity > 0.0){ update += sourceVec[nPhaseIdx] / porosity;}
+ update += sourceVec[nPhaseIdx] / porosity;
break;
}
}
@@ -980,15 +982,15 @@ void FVSaturation2P<TypeTag>::getSource(Scalar& update, const Element& element,
{
//add cflFlux for time-stepping
this->evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosity_[wPhaseIdx], viscosity_[nPhaseIdx],
- (sourceVec[wPhaseIdx] + sourceVec[nPhaseIdx]) * volume, element);
+ (sourceVec[wPhaseIdx] + sourceVec[nPhaseIdx]) * -1 * volume, element);
break;
}
default:
{
//add cflFlux for time-stepping
- this->evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosity_[wPhaseIdx], viscosity_[nPhaseIdx], sourceVec[wPhaseIdx] * volume, element,
+ this->evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosity_[wPhaseIdx], viscosity_[nPhaseIdx], sourceVec[wPhaseIdx] * -1 * volume, element,
wPhaseIdx);
- this->evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosity_[wPhaseIdx], viscosity_[nPhaseIdx], sourceVec[nPhaseIdx] * volume, element,
+ this->evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosity_[wPhaseIdx], viscosity_[nPhaseIdx], sourceVec[nPhaseIdx] * -1 * volume, element,
nPhaseIdx);
break;
}
diff --git a/dumux/decoupled/common/impetproperties.hh b/dumux/decoupled/common/impetproperties.hh
index 7e94616c3d7b26c47ed1ef36458d025d2eeee94b..f0d2b66f3da00ef805d951f21f467c8f73512851 100644
--- a/dumux/decoupled/common/impetproperties.hh
+++ b/dumux/decoupled/common/impetproperties.hh
@@ -62,7 +62,8 @@ NEW_PROP_TAG(ImpetIterationFlag); //!< Flag to switch the iteration type of the
NEW_PROP_TAG(ImpetIterationNumber); //!< Number of iterations if IMPET iterations are enabled by the IterationFlags
NEW_PROP_TAG(ImpetMaximumDefect); //!< Maximum Defect if IMPET iterations are enabled by the IterationFlags
NEW_PROP_TAG(ImpetRelaxationFactor); //!< Used for IMPET iterations
-NEW_PROP_TAG(ImpetDtVariationRestrictionFactor);//!<Restricts the variation of a new time step compared to the old one
+NEW_PROP_TAG(ImpetDtVariationRestrictionFactor);
+NEW_PROP_TAG(ImpetPorosityThreshold);
//forward declaration!
NEW_PROP_TAG( Model );//! The model of the specific problem
@@ -85,6 +86,7 @@ SET_INT_PROP(IMPET, ImpetIterationNumber, 2);
SET_SCALAR_PROP(IMPET, ImpetMaximumDefect, 1e-5);
SET_SCALAR_PROP(IMPET, ImpetRelaxationFactor, 1.0);//!< 1 = new solution is new solution, 0 = old solution is new solution
SET_SCALAR_PROP(IMPET, ImpetDtVariationRestrictionFactor, std::numeric_limits<double>::max());
+SET_SCALAR_PROP(IMPET, ImpetPorosityThreshold, 1e-6);
}
}
diff --git a/dumux/decoupled/common/transportproperties.hh b/dumux/decoupled/common/transportproperties.hh
index 25c5dcea2d53ce90c1bff3a8defa4e8695f04b67..506ea1bb254ba338265df1be9bf11c2890e6c39a 100644
--- a/dumux/decoupled/common/transportproperties.hh
+++ b/dumux/decoupled/common/transportproperties.hh
@@ -50,7 +50,8 @@ NEW_PROP_TAG( TransportSolutionType);
NEW_PROP_TAG( EvalCflFluxFunction ); //!< Type of the evaluation of the CFL-condition
NEW_PROP_TAG( ImpetCFLFactor );
NEW_PROP_TAG( ImpetSwitchNormals );
-NEW_PROP_TAG(ImpetDtVariationRestrictionFactor);//!<Restricts the variation of a new time step compared to the old one
+NEW_PROP_TAG(ImpetPorosityThreshold);
+NEW_PROP_TAG(ImpetDtVariationRestrictionFactor);
SET_SCALAR_PROP(Transport, ImpetCFLFactor, 1.0);
@@ -72,6 +73,7 @@ SET_PROP(Transport, TransportSolutionType)
typedef typename SolutionType::ScalarSolution type;//!<type for vector of scalar properties
};
SET_SCALAR_PROP(Transport, ImpetDtVariationRestrictionFactor, std::numeric_limits<double>::max());
+SET_SCALAR_PROP(Transport, ImpetPorosityThreshold, 1e-6);
}
}