diff --git a/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh b/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
index cf972f8adcaab0b43287d80974b533a0dd3f3abb..51348de680dd98e988ad4b78b58514f87dc29e25 100644
--- a/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
+++ b/dumux/decoupled/2p/transport/fv/evalcflfluxcoats.hh
@@ -24,7 +24,7 @@
* @brief CFL-flux-function to evaluate a CFL-Condition after Coats 2003
*/
-#include "evalcflfluxdefault.hh"
+#include "evalcflflux.hh"
namespace Dumux
{
@@ -34,10 +34,10 @@ namespace Dumux
* tparam TypeTag The problem TypeTag
*/
template<class TypeTag>
-class EvalCflFluxCoats: public EvalCflFluxDefault<TypeTag>
+class EvalCflFluxCoats: public EvalCflFlux<TypeTag>
{
private:
- typedef EvalCflFluxDefault<TypeTag> ParentType;typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
@@ -54,50 +54,60 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
enum
- {
- dim = GridView::dimension, dimWorld = GridView::dimensionworld
- };
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
enum
- {
- wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
- eqIdxPress = Indices::pressEqIdx,
- eqIdxSat = Indices::satEqIdx
- };
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ eqIdxPress = Indices::pressEqIdx,
+ eqIdxSat = Indices::satEqIdx,
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
enum
- {
- pw = Indices::pressureW,
- pn = Indices::pressureNW,
- pglobal = Indices::pressureGlobal,
- vw = Indices::velocityW,
- vn = Indices::velocityNW,
- vt = Indices::velocityTotal,
- Sw = Indices::saturationW,
- Sn = Indices::saturationNW
- };
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ vw = Indices::velocityW,
+ vn = Indices::velocityNW,
+ vt = Indices::velocityTotal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW
+ };
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
public:
+ //! \brief Initializes the cfl-flux-model
+ void initialize()
+ {
+ ElementIterator element = problem_.gridView().template begin<0>();
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, problem_.referencePressure(*element));
+ fluidState.setPressure(nPhaseIdx, problem_.referencePressure(*element));
+ fluidState.setTemperature(problem_.temperature(*element));
+ fluidState.setSaturation(wPhaseIdx, 1.);
+ fluidState.setSaturation(nPhaseIdx, 0.);
+ density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
+ density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
+ }
/*! \brief adds a flux to the cfl-criterion evaluation
*
* \copydetails EvalCflFlux::addFlux(Scalar&,Scalar&,Scalar&,Scalar&,Scalar,const Element&,int)
*/
void addFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux,
- const Element& element, int phaseIdx = -1)
+ const Element& element, int phaseIdx = -1)
{
-// if (hasHangingNode_)
-// {
-// return;
-// }
-
- ParentType::addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, flux, element, phaseIdx);
+ addDefaultFlux(flux, phaseIdx);
}
/*! \brief adds a flux to the cfl-criterion evaluation
@@ -105,404 +115,546 @@ public:
* \copydetails EvalCflFlux::addFlux(Scalar&,Scalar&,Scalar&,Scalar&,Scalar,const Intersection&,int)
*/
void addFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux,
- const Intersection& intersection, int phaseIdx = -1)
+ const Intersection& intersection, int phaseIdx = -1)
{
- Scalar parentMob = 0.5;
- Scalar parentVis = 1;
- // ParentType::addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, flux, intersection, phaseIdx);
- ParentType::addFlux(parentMob, parentMob, parentVis, parentVis, flux, intersection, phaseIdx);
-
- if (hasHangingNode_)
- {
- return;
- }
-
- ElementPointer element = intersection.inside();
-
- //coordinates of cell center
- const GlobalPosition& globalPos = element->geometry().center();
-
- // cell index
- int globalIdxI = problem_.variables().index(*element);
+ addDefaultFlux(flux, phaseIdx);
+ addCoatsFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, flux, intersection, phaseIdx);
+ }
- CellData& cellDataI = problem_.variables().cellData(globalIdxI);
+ /*! \brief Returns the CFL flux-function
+ *
+ * \copydetails EvalCflFlux::getCFLFluxFunction(const Element&)
+ */
+ Scalar getCFLFluxFunction(const Element& element)
+ {
+ Scalar cflFluxDefault = getCFLFluxFunctionDefault();
- int indexInInside = intersection.indexInInside();
+ // std::cout<<"cflFluxFunctionCoats_ = "<<cflFluxFunctionCoats_<<", cflFluxDefault = "<<cflFluxDefault<<"\n";
- Scalar satI = cellDataI.saturation(wPhaseIdx);
- Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
- Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx);
+ if (std::isnan(cflFluxFunctionCoats_) || std::isinf(cflFluxFunctionCoats_))
+ {
+ return 0.99 / cflFluxDefault;
+ }
+ else if (hasHangingNode_)
+ {
+ return 0.99 / cflFluxDefault;
+ }
+ else if (cflFluxFunctionCoats_ <= 0)
+ {
+ return 0.99 / cflFluxDefault;
+ }
+ else
+ {
+ return 1.0 / cflFluxFunctionCoats_;
+ }
+ }
- Scalar dPcdSI = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*element), satI);
+ /*! \brief Returns the CFL time-step
+ *
+ * \copydetails EvalCflFlux::getDt(const Element&)
+ */
+ 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());
+ }
+ //! \brief Resets the Timestep-estimator
+ void reset()
+ {
+ cflFluxFunctionCoats_ = 0;
+ hasHangingNode_ = false;
+ fluxWettingOut_ = 0;
+ fluxNonwettingOut_ = 0;
+ fluxIn_ = 0;
+ fluxOut_ = 0;
+ }
- const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
+ /*! \brief Constructs an EvalCflFluxDefault object
+ *
+ * \param problem A problem type object
+ */
+ EvalCflFluxCoats(Problem& problem) :
+ problem_(problem), epsDerivative_(5e-3), threshold_(1e-12)
+ {
+ reset();
+ density_[wPhaseIdx] = 0.;
+ density_[nPhaseIdx] = 0.;
+ }
- if (intersection.neighbor())
+private:
+ Scalar getCFLFluxFunctionDefault()
+ {
+ if (std::isnan(fluxIn_) || std::isinf(fluxIn_))
{
- ElementPointer neighbor = intersection.outside();
-
- const GlobalPosition& globalPosNeighbor = neighbor->geometry().center();
+ fluxIn_ = 1e-100;
+ }
- // distance vector between barycenters
- Dune::FieldVector < Scalar, dimWorld > distVec = globalPosNeighbor - globalPos;
+ Scalar cFLFluxIn = fluxIn_;
- // compute distance between cell centers
-// Scalar dist = distVec.two_norm();
- Scalar dist = std::abs(distVec*unitOuterNormal);
- int globalIdxJ = problem_.variables().index(*neighbor);
+ Scalar cFLFluxOut = 0;
- CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
+ if (velocityType_ == vt)
+ {
+ if (std::isnan(fluxOut_) || std::isinf(fluxOut_))
+ {
+ fluxOut_ = 1e-100;
+ }
- //calculate potential gradients
- if (element.level() < neighbor.level())
+ cFLFluxOut = fluxOut_;
+ }
+ else
+ {
+ if (std::isnan(fluxWettingOut_) || std::isinf(fluxWettingOut_))
{
- hasHangingNode_ = true;
- cflFluxFunction_ = 0;
- return;
+ fluxWettingOut_ = 1e-100;
+ }
+ if (std::isnan(fluxNonwettingOut_) || std::isinf(fluxNonwettingOut_))
+ {
+ fluxNonwettingOut_ = 1e-100;
}
- Scalar satJ = cellDataJ.saturation(wPhaseIdx);
- Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx);
- Scalar lambdaNWJ = cellDataI.mobility(nPhaseIdx);
+ cFLFluxOut = std::max(fluxWettingOut_, fluxNonwettingOut_);
+ }
- Scalar dPcdSJ = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*neighbor), satJ);
- // compute vectorized permeabilities
- DimMatrix meanPermeability(0);
+ //determine timestep
+ Scalar cFLFluxFunction = std::max(cFLFluxIn, cFLFluxOut);
- problem_.spatialParams().meanK(meanPermeability,
- problem_.spatialParams().intrinsicPermeability(*element),
- problem_.spatialParams().intrinsicPermeability(*neighbor));
+ return cFLFluxFunction;
+ }
+
+ void addDefaultFlux(Scalar flux,int phaseIdx);
- Dune::FieldVector < Scalar, dim > permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
+ void addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux,
+ const Intersection& intersection, int phaseIdx);
- Scalar transmissibility = (unitOuterNormal * permeability) * intersection.geometry().volume() / dist;
+ Problem& problem_;//problem data
+ Scalar cflFluxFunctionCoats_;
+ Scalar fluxWettingOut_;
+ Scalar fluxNonwettingOut_;
+ Scalar fluxOut_;
+ Scalar fluxIn_;
+ bool hasHangingNode_;
+ 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 satUpw = 0;
- if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside))
+template<class TypeTag>
+void EvalCflFluxCoats<TypeTag>::addDefaultFlux(Scalar flux, int phaseIdx)
+{
+ switch (phaseIdx)
+ {
+ case wPhaseIdx:
+ {
+ //for time step criterion
+ if (flux >= 0)
{
- satUpw = std::max(satI, 0.0);
+ fluxWettingOut_ += flux;
}
- else
+ if (flux < 0)
{
- satUpw = std::max(satJ, 0.0);
+ fluxIn_ -= flux;
}
- Scalar dS = eps_;
+ break;
+ }
- Scalar satPlus = satUpw + eps_;
- Scalar satMinus = satUpw;
- if (satMinus - eps_ > 0.0)
+ //for time step criterion if the non-wetting phase velocity is used
+ case nPhaseIdx:
+ {
+ if (flux >= 0)
{
- satMinus -= eps_;
- dS += eps_;
+ fluxNonwettingOut_ += flux;
+ }
+ if (flux < 0)
+ {
+ fluxIn_ -= flux;
}
- Scalar dLambdaWdS = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighbor), std::abs(satPlus)) / viscosityW;
- dLambdaWdS -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighbor), std::abs(satMinus)) / viscosityW;
- dLambdaWdS /= (dS);
-
- if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
+ break;
+ }
+ default:
+ {
+ if (flux >= 0)
{
- satUpw = std::max(1 - satI, 0.0);
+ fluxOut_ += flux;
}
- else
+ if (flux < 0)
{
- satUpw = std::max(1 - satJ, 0.0);
+ fluxIn_ -= flux;
}
- dS = eps_;
+ break;
+ }
+ }
+}
- satPlus = satUpw + eps_;
- satMinus = satUpw;
- if (satMinus - eps_ > 0.0)
- {
- satMinus -= eps_;
- dS += eps_;
- }
+/*! \brief adds a flux to the cfl-criterion evaluation
+ *
+ * \copydetails EvalCflFlux::addFlux(Scalar&,Scalar&,Scalar&,Scalar&,Scalar,const Intersection&,int)
+ */
+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);
- Scalar dLambdaNWdS = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighbor), satPlus) / viscosityNW;
- dLambdaNWdS -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighbor), satMinus) / viscosityNW;
- dLambdaNWdS /= (dS);
+ if (lambdaT <= 0.0)
+ {
+ return;
+ }
- Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWJ);
- Scalar lambdaNWCap = 0.5 * (lambdaNWI + lambdaNWJ);
+ ElementPointer element = intersection.inside();
- if (phaseIdx == wPhaseIdx)
- {
- if (flux != 0)
- {
- cflFluxFunction_ += lambdaNW * dLambdaWdS * std::abs(flux) / (lambdaW * (lambdaW + lambdaNW));
- }
+ //coordinates of cell center
+ const GlobalPosition& globalPos = element->geometry().center();
- cflFluxFunction_ -= transmissibility * lambdaWCap * lambdaNWCap * (dPcdSI + dPcdSJ) / (lambdaW
- + lambdaNW);
- }
- else if (phaseIdx == nPhaseIdx)
- {
- if (flux != 0)
- {
- cflFluxFunction_ -= lambdaW * dLambdaNWdS * std::abs(flux) / (lambdaNW * (lambdaW + lambdaNW));
- }
- }
- else
+ // cell index
+ int globalIdxI = problem_.variables().index(*element);
+
+ CellData& cellDataI = problem_.variables().cellData(globalIdxI);
+
+ int indexInInside = intersection.indexInInside();
+
+ Scalar satI = cellDataI.saturation(wPhaseIdx);
+ Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
+ Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx);
+
+ Scalar dPcdSI = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*element), satI);
+
+ const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
+
+ if (intersection.neighbor())
+ {
+ ElementPointer neighbor = intersection.outside();
+
+ const GlobalPosition& globalPosNeighbor = neighbor->geometry().center();
+
+ // distance vector between barycenters
+ Dune::FieldVector < Scalar, dimWorld > distVec = globalPosNeighbor - globalPos;
+
+ // compute distance between cell centers
+ // Scalar dist = distVec.two_norm();
+ Scalar dist = std::abs(distVec*unitOuterNormal);
+
+ int globalIdxJ = problem_.variables().index(*neighbor);
+
+ CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
+
+ //calculate potential gradients
+ if (element.level() != neighbor.level())
+ {
+ hasHangingNode_ = true;
+ cflFluxFunctionCoats_ = 0;
+ return;
+ }
+
+ Scalar satJ = cellDataJ.saturation(wPhaseIdx);
+ Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx);
+ Scalar lambdaNWJ = cellDataI.mobility(nPhaseIdx);
+
+ Scalar dPcdSJ = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*neighbor), satJ);
+
+ // compute vectorized permeabilities
+ DimMatrix meanPermeability(0);
+
+ problem_.spatialParams().meanK(meanPermeability,
+ problem_.spatialParams().intrinsicPermeability(*element),
+ problem_.spatialParams().intrinsicPermeability(*neighbor));
+
+ Dune::FieldVector < Scalar, dim > permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
+
+ Scalar transmissibility = (unitOuterNormal * permeability) * intersection.geometry().volume() / dist;
+
+ Scalar satUpw = 0;
+ if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside))
+ {
+ satUpw = std::max(satI, 0.0);
+ }
+ else
+ {
+ satUpw = std::max(satJ, 0.0);
+ }
+
+ Scalar dS = epsDerivative_;
+
+ Scalar satPlus = satUpw + epsDerivative_;
+ Scalar satMinus = satUpw;
+ if (satMinus - epsDerivative_ > 0.0)
+ {
+ satMinus -= epsDerivative_;
+ dS += epsDerivative_;
+ }
+
+ Scalar dLambdaWdS = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighbor), std::abs(satPlus)) / viscosityW;
+ dLambdaWdS -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighbor), std::abs(satMinus)) / viscosityW;
+ dLambdaWdS /= (dS);
+
+ if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
+ {
+ satUpw = std::max(1 - satI, 0.0);
+ }
+ else
+ {
+ satUpw = std::max(1 - satJ, 0.0);
+ }
+
+ dS = epsDerivative_;
+
+ satPlus = satUpw + epsDerivative_;
+ satMinus = satUpw;
+ if (satMinus - epsDerivative_ > 0.0)
+ {
+ satMinus -= epsDerivative_;
+ dS += epsDerivative_;
+ }
+
+ Scalar dLambdaNWdS = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighbor), satPlus) / viscosityNW;
+ dLambdaNWdS -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighbor), satMinus) / viscosityNW;
+ dLambdaNWdS /= (dS);
+
+ Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWJ);
+ Scalar lambdaNWCap = 0.5 * (lambdaNWI + lambdaNWJ);
+
+ if (phaseIdx == wPhaseIdx)
+ {
+ Scalar pressDiff = cellDataI.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx);
+ cflFluxFunctionCoats_+= transmissibility * lambdaNW * dLambdaWdS * std::abs(pressDiff) / lambdaT;
+
+ cflFluxFunctionCoats_ -= transmissibility * lambdaWCap * lambdaNWCap * (dPcdSI + dPcdSJ) / lambdaT;
+ }
+ else if (phaseIdx == nPhaseIdx)
+ {
+ Scalar pressDiff = cellDataI.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx);
+ cflFluxFunctionCoats_ -= transmissibility * lambdaW * dLambdaNWdS * std::abs(pressDiff) / lambdaT;
+ }
+ else
+ {
+ if (flux != 0)
{
- if (flux != 0)
+ switch (saturationType_)
{
- switch (saturationType_)
- {
- case Sw:
+ case Sw:
{
- cflFluxFunction_ += dLambdaWdS / (dLambdaWdS + dLambdaNWdS) * std::abs(flux);
- break;
+ cflFluxFunctionCoats_ += dLambdaWdS / (dLambdaWdS + dLambdaNWdS) * std::abs(flux);
+ break;
}
- case Sn:
+ case Sn:
{
- cflFluxFunction_ += dLambdaNWdS / (dLambdaWdS + dLambdaNWdS) * std::abs(flux);
- break;
- }
+ cflFluxFunctionCoats_ += dLambdaNWdS / (dLambdaWdS + dLambdaNWdS) * std::abs(flux);
+ break;
}
}
}
}
- else
- {
- // center of face in global coordinates
- GlobalPosition globalPosFace = intersection.geometry().center();
+ }
+ else
+ {
+ // center of face in global coordinates
+ GlobalPosition globalPosFace = intersection.geometry().center();
- //get boundary type
- BoundaryTypes bcType;
- problem_.boundaryTypes(bcType, intersection);
+ //get boundary type
+ BoundaryTypes bcType;
+ problem_.boundaryTypes(bcType, intersection);
- // distance vector between barycenters
- Dune::FieldVector < Scalar, dimWorld > distVec = globalPosFace - globalPos;
+ // distance vector between barycenters
+ Dune::FieldVector < Scalar, dimWorld > distVec = globalPosFace - globalPos;
- // compute distance between cell centers
- Scalar dist = distVec.two_norm();
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
- //permeability vector at boundary
- // compute vectorized permeabilities
- DimMatrix meanPermeability(0);
+ //permeability vector at boundary
+ // compute vectorized permeabilities
+ DimMatrix meanPermeability(0);
- problem_.spatialParams().meanK(meanPermeability,
- problem_.spatialParams().intrinsicPermeability(*element));
+ problem_.spatialParams().meanK(meanPermeability,
+ problem_.spatialParams().intrinsicPermeability(*element));
- Dune::FieldVector<Scalar, dim> permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
- Scalar satWBound = cellDataI.saturation(wPhaseIdx);
- if (bcType.isDirichlet(eqIdxSat))
+ Scalar transmissibility = (unitOuterNormal * permeability) * intersection.geometry().volume() / dist;
+
+ Scalar satWBound = cellDataI.saturation(wPhaseIdx);
+ if (bcType.isDirichlet(eqIdxSat))
+ {
+ PrimaryVariables bcValues;
+ problem_.dirichlet(bcValues, intersection);
+ switch (saturationType_)
{
- PrimaryVariables bcValues;
- problem_.dirichlet(bcValues, intersection);
- switch (saturationType_)
- {
- case Sw:
+ case Sw:
{
satWBound = bcValues[eqIdxSat];
break;
}
- case Sn:
+ case Sn:
{
satWBound = 1 - bcValues[eqIdxSat];
break;
}
- default:
+ default:
{
DUNE_THROW(Dune::RangeError, "saturation type not implemented");
break;
}
- }
-
- }
-
- Scalar dPcdSBound = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*element), satWBound);
-
- Scalar lambdaWBound = 0;
- Scalar lambdaNWBound = 0;
-
- Scalar temperature = problem_.temperature(*element);
- Scalar referencePressure = problem_.referencePressure(*element);
- FluidState fluidState;
- fluidState.setPressure(wPhaseIdx, referencePressure);
- fluidState.setPressure(nPhaseIdx, referencePressure);
- fluidState.setTemperature(temperature);
-
- Scalar viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx);
- Scalar viscosityNWBound =
- FluidSystem::viscosity(fluidState, nPhaseIdx);
- lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satWBound) / viscosityWBound;
- lambdaNWBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satWBound) / viscosityNWBound;
-
- Scalar transmissibility = (unitOuterNormal * permeability) * intersection.geometry().volume() / dist;
-
- Scalar satUpw = 0;
- if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside))
- {
- satUpw = std::max(satI, 0.0);
- }
- else
- {
- satUpw = std::max(satWBound, 0.0);
}
- Scalar dS = eps_;
+ }
- Scalar satPlus = satUpw + eps_;
- Scalar satMinus = satUpw;
- if (satMinus - eps_ > 0.0)
+ Scalar potWBound = cellDataI.pressure(wPhaseIdx);
+ Scalar potNWBound = cellDataI.pressure(nPhaseIdx);
+ Scalar gdeltaZ = (problem_.bboxMax()-globalPosFace) * problem_.gravity();
+ if (bcType.isDirichlet(eqIdxPress))
+ {
+ PrimaryVariables bcValues;
+ problem_.dirichlet(bcValues, intersection);
+ switch (pressureType_)
{
- satMinus -= eps_;
- dS += eps_;
+ case pw:
+ {
+ potWBound = bcValues[eqIdxPress] + density_[wPhaseIdx] * gdeltaZ;
+ potNWBound = bcValues[eqIdxPress] + MaterialLaw::pC(problem_.spatialParams().materialLawParams(*element), satWBound) + density_[nPhaseIdx] * gdeltaZ;
+ break;
+ }
+ case pn:
+ {
+ potWBound = bcValues[eqIdxPress] - MaterialLaw::pC(problem_.spatialParams().materialLawParams(*element), satWBound) + density_[wPhaseIdx] * gdeltaZ;
+ potNWBound = bcValues[eqIdxPress] + density_[nPhaseIdx] * gdeltaZ;
+ break;
+ }
+ default:
+ {
+ DUNE_THROW(Dune::RangeError, "pressure type not implemented");
+ break;
+ }
}
+ }
+ else if (bcType.isNeumann(eqIdxPress))
+ {
+ PrimaryVariables bcValues;
+ problem_.neumann(bcValues, intersection);
- Scalar dLambdaWdS = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satPlus) / viscosityW;
- dLambdaWdS -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satMinus) / viscosityW;
- dLambdaWdS /= (dS);
-
- if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
+ if (lambdaW != 0 && bcValues[wPhaseIdx] != 0)
{
- satUpw = std::max(1 - satI, 0.0);
+ potWBound += bcValues[wPhaseIdx] / (transmissibility * lambdaW);
}
- else
+ if (lambdaNW != 0 && bcValues[nPhaseIdx] != 0)
{
- satUpw = std::max(1 - satWBound, 0.0);
+ potNWBound += bcValues[nPhaseIdx] / (transmissibility * lambdaNW);
}
+ }
- dS = eps_;
+ Scalar dPcdSBound = MaterialLaw::dpC_dSw(problem_.spatialParams().materialLawParams(*element), satWBound);
- satPlus = satUpw + eps_;
- satMinus = satUpw;
- if (satMinus - eps_ > 0.0)
- {
- satMinus -= eps_;
- dS += eps_;
- }
+ Scalar lambdaWBound = 0;
+ Scalar lambdaNWBound = 0;
- Scalar dLambdaNWdS = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satPlus) / viscosityNW;
- dLambdaNWdS -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satMinus) / viscosityNW;
- dLambdaNWdS /= (dS);
+ Scalar temperature = problem_.temperature(*element);
+ Scalar referencePressure = problem_.referencePressure(*element);
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, referencePressure);
+ fluidState.setPressure(nPhaseIdx, referencePressure);
+ fluidState.setTemperature(temperature);
- Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWBound);
- Scalar lambdaNWCap = 0.5 * (lambdaNWI + lambdaNWBound);
+ Scalar viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ Scalar viscosityNWBound =
+ FluidSystem::viscosity(fluidState, nPhaseIdx);
+ lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satWBound) / viscosityWBound;
+ lambdaNWBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satWBound) / viscosityNWBound;
- if (phaseIdx == wPhaseIdx)
- {
- if (flux != 0)
- {
- cflFluxFunction_ += lambdaNW * dLambdaWdS * std::abs(flux) / (lambdaW * (lambdaW + lambdaNW));
- }
-
- cflFluxFunction_ -= transmissibility * lambdaWCap * lambdaNWCap * (dPcdSI + dPcdSBound) / (lambdaW
- + lambdaNW);
- }
- else if (phaseIdx == nPhaseIdx)
- {
- if (flux != 0)
- {
- cflFluxFunction_ -= lambdaW * dLambdaNWdS * std::abs(flux) / (lambdaNW * (lambdaW + lambdaNW));
- }
- }
- else
- {
- if (flux != 0)
- {
- switch (saturationType_)
- {
- case Sw:
- {
- cflFluxFunction_ += dLambdaWdS / (dLambdaWdS + dLambdaNWdS) * std::abs(flux);
- break;
- }
- case Sn:
- {
- cflFluxFunction_ += dLambdaNWdS / (dLambdaWdS + dLambdaNWdS) * std::abs(flux);
- break;
- }
- }
- }
- }
+ Scalar satUpw = 0;
+ if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside))
+ {
+ satUpw = std::max(satI, 0.0);
+ }
+ else
+ {
+ satUpw = std::max(satWBound, 0.0);
}
- }
- /*! \brief Returns the CFL flux-function
- *
- * \copydetails EvalCflFlux::getCFLFluxFunction(const Element&)
- */
- Scalar getCFLFluxFunction(const Element& element)
- {
- Scalar cflFluxDefault = 1 / ParentType::getCFLFluxFunction(element);
+ Scalar dS = epsDerivative_;
- if (std::isnan(cflFluxFunction_) || std::isinf(cflFluxFunction_) || cflFluxFunction_ > 100 * cflFluxDefault)
+ Scalar satPlus = satUpw + epsDerivative_;
+ Scalar satMinus = satUpw;
+ if (satMinus - epsDerivative_ > 0.0)
{
- cflFluxFunction_ = 0;
+ satMinus -= epsDerivative_;
+ dS += epsDerivative_;
}
-// if (cflFluxDefault > 1)
-// std::cout << "cflFluxDefault = " << cflFluxDefault << "\n";
-//
-// if (cflFluxFunction_ > 1)
-// std::cout << "cflFluxFunction_ = " << cflFluxFunction_ << "\n";
-
- Scalar returnValue = std::max(cflFluxFunction_, cflFluxDefault);
- reset();
+ Scalar dLambdaWdS = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satPlus) / viscosityW;
+ dLambdaWdS -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satMinus) / viscosityW;
+ dLambdaWdS /= (dS);
- if (returnValue > 0 && !hasHangingNode_)
+ if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
{
- return (returnValue == cflFluxDefault) ? 0.95 / returnValue : 1 / returnValue;
+ satUpw = std::max(1 - satI, 0.0);
}
else
{
- return 0.95 / cflFluxDefault;
+ satUpw = std::max(1 - satWBound, 0.0);
}
- }
- /*! \brief Returns the CFL time-step
- *
- * \copydetails EvalCflFlux::getDt(const Element&)
- */
- 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());
- }
+ dS = epsDerivative_;
- /*! \brief Constructs an EvalCflFluxDefault object
- *
- * \param problem A problem type object
- */
- EvalCflFluxCoats(Problem& problem) :
- ParentType(problem), problem_(problem), eps_(5e-3)
- {
- reset();
- }
+ satPlus = satUpw + epsDerivative_;
+ satMinus = satUpw;
+ if (satMinus - epsDerivative_ > 0.0)
+ {
+ satMinus -= epsDerivative_;
+ dS += epsDerivative_;
+ }
-protected:
- //! resets the accumulated CFL-fluxes to zero
- void reset()
- {
- ParentType::reset();
- cflFluxFunction_ = 0;
- hasHangingNode_ = false;
- }
+ Scalar dLambdaNWdS = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satPlus) / viscosityNW;
+ dLambdaNWdS -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satMinus) / viscosityNW;
+ dLambdaNWdS /= (dS);
-private:
+ Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWBound);
+ Scalar lambdaNWCap = 0.5 * (lambdaNWI + lambdaNWBound);
+ if (phaseIdx == wPhaseIdx)
+ {
+ Scalar potDiff = cellDataI.pressure(wPhaseIdx) - potWBound;
+ cflFluxFunctionCoats_ += transmissibility * lambdaNW * dLambdaWdS * std::abs(potDiff) / lambdaT;
- Problem& problem_;//problem data
- Scalar cflFluxFunction_;
- bool hasHangingNode_;
- 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 eps_;
-};
+ cflFluxFunctionCoats_ -= transmissibility * lambdaWCap * lambdaNWCap * (dPcdSI + dPcdSBound) / lambdaT;
+ }
+ else if (phaseIdx == nPhaseIdx)
+ {
+ Scalar potDiff = cellDataI.pressure(nPhaseIdx) - potNWBound;
+ cflFluxFunctionCoats_ -= transmissibility * lambdaW * dLambdaNWdS * std::abs(potDiff) / lambdaT;
+ }
+ else
+ {
+ if (flux != 0)
+ {
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ cflFluxFunctionCoats_ += dLambdaWdS / (dLambdaWdS + dLambdaNWdS) * std::abs(flux);
+ break;
+ }
+ case Sn:
+ {
+ cflFluxFunctionCoats_ += dLambdaNWdS / (dLambdaWdS + dLambdaNWdS) * std::abs(flux);
+ break;
+ }
+ }
+ }
+ }
+ }
+}
}