diff --git a/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh b/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
index 79ffb93e304e9eb4f7a074f3b100e8e5811feb75..e91accc2aa3c7f7586c9c59f5411eaa3940547a2 100644
--- a/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
+++ b/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
@@ -51,7 +51,7 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
@@ -59,6 +59,9 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(CellData)) CellData;
+
+
enum
{
dim = GridView::dimension, dimWorld = GridView::dimensionworld
@@ -99,6 +102,7 @@ public:
break;
}
int globalIdxI = problem_.variables().index(element);
+ CellData& CellDataI = problem_.variables().cellData(globalIdxI);
// get geometry type of face
//Dune::GeometryType faceGT = isIt->geometryInInside().type();
@@ -113,17 +117,19 @@ public:
if (preComput_)
{
- mobilityWI = problem_.variables().mobilityWetting(globalIdxI);
- mobilityNWI = problem_.variables().mobilityNonwetting(globalIdxI);
+ mobilityWI = CellDataI.mobility(wPhaseIdx);
+ mobilityNWI = CellDataI.mobility(nPhaseIdx);
}
else
{
FluidState fluidState;
- fluidState.update(satI, referencePressure, referencePressure, temperature);
+ fluidState.setPressure(wPhaseIdx, referencePressure);
+ fluidState.setPressure(nPhaseIdx, referencePressure);
+ fluidState.setTemperature(temperature);
mobilityWI = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satI);
- mobilityWI /= FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
+ mobilityWI /= FluidSystem::viscosity(fluidState, wPhaseIdx);
mobilityNWI = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(element), satI);
- mobilityNWI /= FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
+ mobilityNWI /= FluidSystem::viscosity(fluidState, nPhaseIdx);
}
FieldMatrix meanPermeability(0);
@@ -134,6 +140,7 @@ public:
ElementPointer neighborPointer = isIt->outside();
int globalIdxJ = problem_.variables().index(*neighborPointer);
+ CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
// neighbor cell center in global coordinates
const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
@@ -158,17 +165,20 @@ public:
//get lambda_bar = lambda_n*f_w
if(preComput_)
{
- mobilityWJ = problem_.variables().mobilityWetting(globalIdxJ);
- mobilityNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
+ mobilityWJ = cellDataJ.mobility(wPhaseIdx);
+ mobilityNWJ = cellDataJ.mobility(nPhaseIdx);
}
else
{
FluidState fluidState;
- fluidState.update(satJ, referencePressure, referencePressure, temperature);
+ fluidState.setPressure(wPhaseIdx, referencePressure);
+ fluidState.setPressure(nPhaseIdx, referencePressure);
+ fluidState.setTemperature(temperature);
+
mobilityWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*neighborPointer), satJ);
- mobilityWJ /= FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
+ mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
mobilityNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*neighborPointer), satJ);
- mobilityNWJ /= FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
+ mobilityNWJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
}
Scalar mobilityWMean = 0.5*(mobilityWI + mobilityWJ);
Scalar mobilityNWMean = 0.5*(mobilityNWI + mobilityNWJ);
@@ -185,11 +195,13 @@ public:
//calculate lambda_n*f_w at the boundary
FluidState fluidState;
- fluidState.update(satJ, referencePressure, referencePressure, temperature);
+ fluidState.setPressure(wPhaseIdx, referencePressure);
+ fluidState.setPressure(nPhaseIdx, referencePressure);
+ fluidState.setTemperature(temperature);
mobilityWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satJ);
- mobilityWJ /= FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
+ mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
mobilityNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(element), satJ);
- mobilityNWJ /= FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
+ mobilityNWJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
Scalar mobWMean = 0.5 * (mobilityWI + mobilityWJ);
Scalar mobNWMean = 0.5 * (mobilityNWI + mobilityNWJ);
diff --git a/dumux/decoupled/2p/transport/fv/evalcflflux_coats.hh b/dumux/decoupled/2p/transport/fv/evalcflflux_coats.hh
index fd294eee9f7e5466fe82e86f1e1295fdfa0b2757..520be40352f3a7670ffd38813774da8781b291d1 100644
--- a/dumux/decoupled/2p/transport/fv/evalcflflux_coats.hh
+++ b/dumux/decoupled/2p/transport/fv/evalcflflux_coats.hh
@@ -52,12 +52,14 @@ private:
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(CellData)) CellData;
+
enum
{
dim = GridView::dimension, dimWorld = GridView::dimensionworld
@@ -122,11 +124,13 @@ public:
// cell index
int globalIdxI = problem_.variables().index(element);
+ CellData& cellDataI = problem_.variables().cellData(globalIdxI);
+
int indexInInside = intersection.indexInInside();
- Scalar satI = problem_.variables().saturation()[globalIdxI];
- Scalar lambdaWI = problem_.variables().mobilityWetting(globalIdxI);
- Scalar lambdaNWI = problem_.variables().mobilityNonwetting(globalIdxI);
+ Scalar satI = cellDataI.saturation(wPhaseIdx);
+ Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
+ Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx);
Scalar dPcdSI = MaterialLaw::dpC_dSw(problem_.spatialParameters().materialLawParams(element), satI);
@@ -148,19 +152,18 @@ public:
int globalIdxJ = problem_.variables().index(neighbor);
+ CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
+
//calculate potential gradients
if (element.level() < neighbor.level())
{
hasHangingNode_ = true;
return;
}
- //get phase potentials
- Scalar potentialW = problem_.variables().potentialWetting(globalIdxI, indexInInside);
- Scalar potentialNW = problem_.variables().potentialNonwetting(globalIdxI, indexInInside);
- Scalar satJ = problem_.variables().saturation()[globalIdxJ];
- Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
+ Scalar satJ = cellDataJ.saturation(wPhaseIdx);
+ Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx);
+ Scalar lambdaNWJ = cellDataI.mobility(nPhaseIdx);
Scalar dPcdSJ = MaterialLaw::dpC_dSw(problem_.spatialParameters().materialLawParams(neighbor), satJ);
@@ -177,7 +180,7 @@ public:
Scalar transmissibility = (unitOuterNormal * permeability) * intersection.geometry().volume() / dist;
Scalar satUpw = 0;
- if (potentialW >= 0)
+ if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside))
{
satUpw = std::max(satI, 0.0);
}
@@ -200,7 +203,7 @@ public:
dLambdaWdS -= MaterialLaw::krw(problem_.spatialParameters().materialLawParams(neighbor), std::abs(satMinus)) / viscosityW;
dLambdaWdS /= (dS);
- if (potentialNW >= 0)
+ if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
{
satUpw = std::max(1 - satI, 0.0);
}
@@ -288,43 +291,33 @@ public:
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
- Scalar satBound = 0;
+ Scalar satWBound = cellDataI.saturation(wPhaseIdx);
if (bcType.isDirichlet(eqIdxSat))
{
PrimaryVariables bcValues;
problem_.dirichlet(bcValues, intersection);
- satBound = bcValues[eqIdxSat];
- }
- else
- {
- satBound = problem_.variables().saturation()[globalIdxI];
- }
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ satWBound = bcValues[eqIdxSat];
+ break;
+ }
+ case Sn:
+ {
+ satWBound = 1 - bcValues[eqIdxSat];
+ break;
+ }
+ default:
+ {
+ DUNE_THROW(Dune::RangeError, "saturation type not implemented");
+ break;
+ }
+ }
- //determine phase saturations from primary saturation variable
- Scalar satW;
- //Scalar satNW;
- switch (saturationType_)
- {
- case Sw:
- {
- satW = satBound;
- //satNW = 1 - satBound;
- break;
- }
- case Sn:
- {
- satW = 1 - satBound;
- //satNW = satBound;
- break;
- }
- default:
- {
- DUNE_THROW(Dune::RangeError, "saturation type not implemented");
- break;
- }
}
- Scalar dPcdSBound = MaterialLaw::dpC_dSw(problem_.spatialParameters().materialLawParams(element), satBound);
+ Scalar dPcdSBound = MaterialLaw::dpC_dSw(problem_.spatialParameters().materialLawParams(element), satWBound);
Scalar lambdaWBound = 0;
Scalar lambdaNWBound = 0;
@@ -332,28 +325,26 @@ public:
Scalar temperature = problem_.temperature(element);
Scalar referencePressure = problem_.referencePressure(element);
FluidState fluidState;
- Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
- Scalar viscosityNWBound =
- FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
- lambdaWBound = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satW) / viscosityWBound;
- lambdaNWBound = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(element), satW) / viscosityNWBound;
-
- Scalar potentialW = 0;
- Scalar potentialNW = 0;
+ fluidState.setPressure(wPhaseIdx, referencePressure);
+ fluidState.setPressure(nPhaseIdx, referencePressure);
+ fluidState.setTemperature(temperature);
- potentialW = problem_.variables().potentialWetting(globalIdxI, indexInInside);
- potentialNW = problem_.variables().potentialNonwetting(globalIdxI, indexInInside);
+ Scalar viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ Scalar viscosityNWBound =
+ FluidSystem::viscosity(fluidState, nPhaseIdx);
+ lambdaWBound = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satWBound) / viscosityWBound;
+ lambdaNWBound = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(element), satWBound) / viscosityNWBound;
Scalar transmissibility = (unitOuterNormal * permeability) * intersection.geometry().volume() / dist;
Scalar satUpw = 0;
- if (potentialW >= 0)
+ if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside))
{
satUpw = std::max(satI, 0.0);
}
else
{
- satUpw = std::max(satBound, 0.0);
+ satUpw = std::max(satWBound, 0.0);
}
Scalar dS = eps_;
@@ -370,13 +361,13 @@ public:
dLambdaWdS -= MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satMinus) / viscosityW;
dLambdaWdS /= (dS);
- if (potentialNW >= 0)
+ if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
{
satUpw = std::max(1 - satI, 0.0);
}
else
{
- satUpw = std::max(1 - satBound, 0.0);
+ satUpw = std::max(1 - satWBound, 0.0);
}
dS = eps_;
diff --git a/dumux/decoupled/2p/transport/fv/evalcflflux_default.hh b/dumux/decoupled/2p/transport/fv/evalcflflux_default.hh
index 8c429e2647ebbe41cfcbc2c4e0a9b934b1a09cd9..2b64a9c7b623e42ca5348be1a68ff20e7455980b 100644
--- a/dumux/decoupled/2p/transport/fv/evalcflflux_default.hh
+++ b/dumux/decoupled/2p/transport/fv/evalcflflux_default.hh
@@ -49,7 +49,7 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
enum
{
@@ -219,18 +219,9 @@ typename EvalCflFluxDefault<TypeTag>::Scalar EvalCflFluxDefault<TypeTag>::getCfl
Scalar volumeCorrectionFactorOutW = 0;
Scalar volumeCorrectionFactorOutNW = 0;
- Scalar sat = problem_.variables().saturation()[problem_.variables().index(element)];
-
- if (saturationType_ == Sw)
- {
- volumeCorrectionFactorOutW = std::max((sat - residualSatW), 1e-2);
- volumeCorrectionFactorOutNW = std::max((1 - sat - residualSatNW), 1e-2);
- }
- if (saturationType_ == Sn)
- {
- volumeCorrectionFactorOutW = std::max((1 - sat - residualSatW), 1e-2);
- volumeCorrectionFactorOutNW = std::max((sat - residualSatNW), 1e-2);
- }
+ Scalar satW = problem_.variables().cellData(problem_.variables().index(element)).saturation(wPhaseIdx);
+ volumeCorrectionFactorOutW = std::max((satW - residualSatW), 1e-2);
+ volumeCorrectionFactorOutNW = std::max((1 - satW - residualSatNW), 1e-2);
//make sure correction is in the right range. If not: force dt to be not min-dt!
if (volumeCorrectionFactorOutW <= 0)
diff --git a/dumux/decoupled/2p/transport/fv/evalcflflux_default_old.hh b/dumux/decoupled/2p/transport/fv/evalcflflux_default_old.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8c429e2647ebbe41cfcbc2c4e0a9b934b1a09cd9
--- /dev/null
+++ b/dumux/decoupled/2p/transport/fv/evalcflflux_default_old.hh
@@ -0,0 +1,268 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2010 by Markus Wolff *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+#ifndef DUMUX_EVALCFLFLUX_DEFAULT_HH
+#define DUMUX_EVALCFLFLUX_DEFAULT_HH
+
+/**
+ * @file
+ * @brief Fluxes to evaluate a CFL-Condition
+ * @author Markus Wolff
+ */
+
+#include "evalcflflux.hh"
+
+namespace Dumux
+{
+/*!\ingroup Saturation2p
+ * @brief Default implementation of cfl-fluxes to evaluate a CFL-Condition
+ *
+ * Compares the maximum of inflow and outflow to the element volume weighted by relative permeability and viscosity ratios.
+ *
+ * Template parameters are:
+
+ - TypeTag PropertyTag of the problem implementation
+ */
+template<class TypeTag>
+class EvalCflFluxDefault: public EvalCflFlux<TypeTag>
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ 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,
+ };
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::Intersection Intersection;
+
+public:
+
+ void addFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux, const Element& element, int phaseIdx = -1)
+ {
+ addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, flux, phaseIdx);
+ }
+
+ void addFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux, const Intersection& intersection, int phaseIdx = -1)
+ {
+ addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, flux, phaseIdx);
+ }
+
+ Scalar getCflFluxFunction(const GlobalPosition& globalPos, const Element& element);
+
+ EvalCflFluxDefault (Problem& problem)
+ : problem_(problem)
+ {
+ reset();
+ }
+
+private:
+
+ void addFlux(Scalar& lambdaW, Scalar& lambdaNW, Scalar& viscosityW, Scalar& viscosityNW, Scalar flux, int phaseIdx = -1)
+ {
+ Scalar krSum = lambdaW * viscosityW + lambdaNW * viscosityNW;
+ Scalar viscosityRatio = 1 - fabs(0.5 - viscosityNW / (viscosityW + viscosityNW));//1 - fabs(viscosityWI-viscosityNWI)/(viscosityWI+viscosityNWI);
+
+ switch (phaseIdx)
+ {
+ case wPhaseIdx:
+ {
+ //for time step criterion
+ if (flux >= 0)
+ {
+ fluxWettingOut_ += flux / (krSum * viscosityRatio);
+ }
+ if (flux < 0)
+ {
+ fluxIn_ -= flux / (krSum * viscosityRatio);
+ }
+
+ break;
+ }
+
+ //for time step criterion if the non-wetting phase velocity is used
+ case nPhaseIdx:
+ {
+ if (flux >= 0)
+ {
+ fluxNonwettingOut_ += flux / (krSum * viscosityRatio);
+ }
+ if (flux < 0)
+ {
+ fluxIn_ -= flux / (krSum * viscosityRatio);
+ }
+
+ break;
+ }
+ default:
+ {
+ if (flux >= 0)
+ {
+ fluxOut_ += flux / (krSum * viscosityRatio);
+ }
+ if (flux < 0)
+ {
+ fluxIn_ -= flux / (krSum * viscosityRatio);
+ }
+
+ break;
+ }
+ }
+ }
+
+ Scalar getCFLFluxIn(int phaseIdx = 0)
+ {
+ if (std::isnan(fluxIn_) || std::isinf(fluxIn_))
+ {
+ fluxIn_ = 1e-100;
+ }
+
+ return fluxIn_;
+ }
+
+ Scalar getCFLFluxOut(int phaseIdx = 0)
+ {
+ if (std::isnan(fluxWettingOut_) || std::isinf(fluxWettingOut_))
+ {
+ fluxWettingOut_ = 1e-100;
+ }
+
+ if (std::isnan(fluxNonwettingOut_) || std::isinf(fluxNonwettingOut_))
+ {
+ fluxNonwettingOut_ = 1e-100;
+ }
+ if (std::isnan(fluxOut_) || std::isinf(fluxOut_))
+ {
+ fluxOut_ = 1e-100;
+ }
+
+ if (phaseIdx == wPhaseIdx)
+ return fluxWettingOut_;
+ else if (phaseIdx == nPhaseIdx)
+ return fluxNonwettingOut_;
+ else
+ return fluxOut_;
+ }
+
+
+protected:
+
+ //! resets the accumulated CFL-fluxes to zero
+ void reset()
+ {
+ fluxWettingOut_ = 0;
+ fluxNonwettingOut_ = 0;
+ fluxIn_ = 0;
+ fluxOut_ = 0;
+ }
+
+private:
+ Problem& problem_;//problem data
+ Scalar fluxWettingOut_;
+ Scalar fluxNonwettingOut_;
+ Scalar fluxOut_;
+ Scalar fluxIn_;
+ static const int velocityType_ = GET_PROP_VALUE(TypeTag, PTAG(VelocityFormulation));
+ static const int saturationType_ = GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation));
+};
+
+template<class TypeTag>
+typename EvalCflFluxDefault<TypeTag>::Scalar EvalCflFluxDefault<TypeTag>::getCflFluxFunction(const GlobalPosition& globalPos, const Element& element)
+{
+ Scalar residualSatW = problem_.spatialParameters().materialLawParams(element).Swr();
+ Scalar residualSatNW = problem_.spatialParameters().materialLawParams(element).Snr();
+
+ // compute dt restriction
+ Scalar volumeCorrectionFactor = 1 - residualSatW - residualSatNW;
+ Scalar volumeCorrectionFactorOutW = 0;
+ Scalar volumeCorrectionFactorOutNW = 0;
+
+ Scalar sat = problem_.variables().saturation()[problem_.variables().index(element)];
+
+ if (saturationType_ == Sw)
+ {
+ volumeCorrectionFactorOutW = std::max((sat - residualSatW), 1e-2);
+ volumeCorrectionFactorOutNW = std::max((1 - sat - residualSatNW), 1e-2);
+ }
+ if (saturationType_ == Sn)
+ {
+ volumeCorrectionFactorOutW = std::max((1 - sat - residualSatW), 1e-2);
+ volumeCorrectionFactorOutNW = std::max((sat - residualSatNW), 1e-2);
+ }
+
+ //make sure correction is in the right range. If not: force dt to be not min-dt!
+ if (volumeCorrectionFactorOutW <= 0)
+ {
+ volumeCorrectionFactorOutW = 1e100;
+ }
+ if (volumeCorrectionFactorOutNW <= 0)
+ {
+ volumeCorrectionFactorOutNW = 1e100;
+ }
+
+ //correct volume
+ Scalar cFLFluxIn = volumeCorrectionFactor / getCFLFluxIn();
+ Scalar cFLFluxOut = 0;
+
+ if (velocityType_ == vt)
+ {
+ cFLFluxOut = volumeCorrectionFactor / getCFLFluxOut(-1);
+ }
+ else
+ {
+ cFLFluxOut = std::min(volumeCorrectionFactorOutW / getCFLFluxOut(wPhaseIdx), volumeCorrectionFactorOutNW / getCFLFluxOut(nPhaseIdx));
+ }
+
+ //determine timestep
+ Scalar cFLFluxFunction = std::min(cFLFluxIn, cFLFluxOut);
+
+ reset();
+
+ return cFLFluxFunction;
+}
+
+}
+
+#endif
diff --git a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
index fc9c62a6dbf3e789faef0d52ec1db340f8397c67..63a6d8178c11559586f4ad1cc6444e914d82ba3c 100644
--- a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
+++ b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
@@ -69,31 +69,30 @@ class FVSaturation2P
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Variables)) Variables;
- typedef Dune::GenericReferenceElements<Scalar, dim>
- ReferenceElementContainer;
- typedef Dune::GenericReferenceElements<Scalar, dim - 1>
- ReferenceElementFaceContainer;
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElementContainer;
+ typedef Dune::GenericReferenceElements<Scalar, dim - 1> ReferenceElementFaceContainer;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(DiffusivePart)) DiffusivePart;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConvectivePart))
- ConvectivePart;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(EvalCflFluxFunction))
- EvalCflFluxFunction;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConvectivePart)) ConvectivePart;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(EvalCflFluxFunction)) EvalCflFluxFunction;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters))
- SpatialParameters;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(WettingPhase)) WettingPhase;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NonwettingPhase)) NonwettingPhase;
typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(CellData)) CellData;
+
enum
{
pw = Indices::pressureW,
@@ -172,8 +171,7 @@ public:
* Additionally to the \a update vector, the recommended time step size \a dt is calculated
* employing a CFL condition.
*/
- void update(const Scalar t, Scalar& dt,
- RepresentationType& updateVec, bool impes);
+ void update(const Scalar t, Scalar& dt, RepresentationType& updateVec, bool impes);
//! Sets the initial solution \f$S_0\f$.
void initialize();
@@ -185,26 +183,71 @@ public:
*/
void updateMaterialLaws(RepresentationType& saturation, bool iterate);
+ void updateSaturationSolution()
+ {
+ int size = problem_.gridView().size(0);
+ for (int i = 0; i < size; i++)
+ {
+ Scalar sat = problem_.variables().primaryVariablesGlobal(satEqIdx)[i];
+ updateSaturationSolution(i, sat);
+ }
+ }
+
+ void updateSaturationSolution(int globalIdx, Scalar sat)
+ {
+ CellData& cellData = problem_.variables().cellData(globalIdx);
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ if (compressibility_)
+ {
+ cellData.fluidState().setSaturation(wPhaseIdx, sat);
+ cellData.fluidState().setSaturation(nPhaseIdx, 1 - sat);
+ }
+ else
+ {
+ cellData.setSaturation(wPhaseIdx, sat);
+ cellData.setSaturation(nPhaseIdx, 1 - sat);
+ }
+ break;
+ }
+ case Sn:
+ {
+ if (compressibility_)
+ {
+ cellData.fluidState().setSaturation(wPhaseIdx, 1 - sat);
+ cellData.fluidState().setSaturation(nPhaseIdx, sat);
+ }
+ else
+ {
+ cellData.setSaturation(wPhaseIdx,1 -sat);
+ cellData.setSaturation(nPhaseIdx, sat);
+ }
+ break;
+ }
+ }
+ }
+
//! \brief Write data files
/* \param name file name */
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
- typename Variables::ScalarSolutionType *saturation =
- writer.allocateManagedBuffer (
- problem_.gridView().size(0));
+ int size = problem_.gridView().size(0);
+ typename Variables::ScalarSolutionType *saturationW = writer.allocateManagedBuffer(size);
+ typename Variables::ScalarSolutionType *saturationN = writer.allocateManagedBuffer(size);
- *saturation = problem_.variables().saturation();
-
- if (saturationType_ == Sw)
- {
- writer.attachCellData(*saturation, "wetting saturation");
- }
- if (saturationType_ == Sn)
+ for (int i = 0; i < size; i++)
{
- writer.attachCellData(*saturation, "nonwetting saturation");
+ CellData& cellData = problem_.variables().cellData(i);
+ (*saturationW)[i] = cellData.saturation(wPhaseIdx);
+ (*saturationN)[i] = cellData.saturation(nPhaseIdx);
}
+ writer.attachCellData(*saturationW, "wetting saturation");
+ writer.attachCellData(*saturationN, "nonwetting saturation");
+
return;
}
@@ -215,14 +258,14 @@ public:
template<class Restarter>
void serialize(Restarter &res)
{
- problem_.variables().serialize<Restarter> (res);
+ problem_.variables().serialize<Restarter>(res);
}
//! Function needed for restart option.
template<class Restarter>
void deserialize(Restarter &res)
{
- problem_.variables().deserialize<Restarter> (res);
+ problem_.variables().deserialize<Restarter>(res);
}
//! Constructs a FVSaturation2P object
@@ -232,20 +275,18 @@ public:
*/
FVSaturation2P(Problem& problem) :
- problem_(problem), switchNormals_(false), threshold_(1e-6)
+ problem_(problem), switchNormals_(false), threshold_(1e-6)
{
if (compressibility_ && velocityType_ == vt)
{
- DUNE_THROW(
- Dune::NotImplemented,
+ DUNE_THROW(Dune::NotImplemented,
"Total velocity - global pressure - model cannot be used with compressible fluids!");
}
if (saturationType_ != Sw && saturationType_ != Sn)
{
DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
}
- if (pressureType_ != pw && pressureType_ != pn && pressureType_
- != pglobal)
+ if (pressureType_ != pw && pressureType_ != pn && pressureType_ != pglobal)
{
DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
}
@@ -272,15 +313,11 @@ private:
ConvectivePart* convectivePart_;
EvalCflFluxFunction* evalCflFluxFunction_;
- static const bool compressibility_ =
- GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility));
+ static const bool compressibility_ = GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility));
;
- static const int saturationType_ =
- GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation));
- static const int velocityType_ =
- GET_PROP_VALUE(TypeTag, PTAG(VelocityFormulation));
- static const int pressureType_ =
- GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation));
+ static const int saturationType_ = GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation));
+ static const int velocityType_ = GET_PROP_VALUE(TypeTag, PTAG(VelocityFormulation));
+ static const int pressureType_ = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation));
bool switchNormals_;
const Scalar threshold_;
@@ -298,8 +335,8 @@ void FVSaturation2P<TypeTag>::indicatorSaturation(RepresentationType &indicator,
{
// 1) calculate Indicator -> min, maxvalues
// Schleife über alle Leaf-Elemente
- for (ElementIterator it = problem_.gridView().template begin<0>();
- it!=problem_.gridView().template end<0>(); ++it)
+ for (ElementIterator it = problem_.gridView().template begin<0>(); it != problem_.gridView().template end<0>();
+ ++it)
{
// Bestimme maximale und minimale Sättigung
// Index des aktuellen Leaf-Elements
@@ -309,23 +346,23 @@ void FVSaturation2P<TypeTag>::indicatorSaturation(RepresentationType &indicator,
// Berechne Verfeinerungsindikator an allen Zellen
IntersectionIterator isend = problem_.gridView().iend(*it);
- for (IntersectionIterator is = problem_.gridView().ibegin(*it); is!= isend; ++is)
+ for (IntersectionIterator is = problem_.gridView().ibegin(*it); is != isend; ++is)
{
- const typename IntersectionIterator::Intersection &intersection =*is;
+ const typename IntersectionIterator::Intersection &intersection = *is;
// Steige aus, falls es sich nicht um einen Nachbarn handelt
- if ( !intersection.neighbor() )
- continue ;
+ if (!intersection.neighbor())
+ continue;
// Greife auf Nachbarn zu
- const Element &outside =*intersection.outside();
- int indexj = problem_.variables().index( outside );
+ const Element &outside = *intersection.outside();
+ int indexj = problem_.variables().index(outside);
// Jede Intersection nur von einer Seite betrachten
- if ( it.level() > outside.level() ||
- (it.level() == outside.level() && indexi<indexj) )
+ if (it.level() > outside.level() || (it.level() == outside.level() && indexi < indexj))
{
- Scalar localdelta = std::abs(problem_.variables().saturation()[indexi][0] - problem_.variables().saturation()[indexj][0]);
+ Scalar localdelta = std::abs(
+ problem_.variables().saturation()[indexi][0] - problem_.variables().saturation()[indexj][0]);
indicator[indexi][0] = std::max(indicator[indexi][0], localdelta);
indicator[indexj][0] = std::max(indicator[indexj][0], localdelta);
}
@@ -334,8 +371,7 @@ void FVSaturation2P<TypeTag>::indicatorSaturation(RepresentationType &indicator,
}
template<class TypeTag>
-void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
- RepresentationType& updateVec, bool impes = false)
+void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationType& updateVec, bool impes = false)
{
if (!impes)
{
@@ -353,10 +389,25 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
BoundaryTypes bcType;
+ ElementIterator eItBegin = problem_.gridView().template begin<0>();
+ Scalar densityW = 0;
+ Scalar densityNW = 0;
+ Scalar viscosityW = 0;
+ Scalar viscosityNW = 0;
+
+ if (!compressibility_)
+ {
+ Scalar temp = problem_.temperature(*eItBegin);
+ Scalar pRef = problem_.referencePressure(*eItBegin);
+ densityW = WettingPhase::density(temp, pRef);
+ densityNW = NonwettingPhase::density(temp, pRef);
+ viscosityW = WettingPhase::viscosity(temp, pRef);
+ viscosityNW = NonwettingPhase::viscosity(temp, pRef);
+ }
+
// compute update vector
- ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt
- != eItEnd; ++eIt)
+ ElementIterator eItEnd = problem_.gridView().template end<0>();
+ for (ElementIterator eIt = eItBegin; eIt != eItEnd; ++eIt)
{
#if HAVE_MPI
if (eIt->partitionType() != Dune::InteriorEntity)
@@ -373,27 +424,26 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
// cell index
int globalIdxI = problem_.variables().index(*eIt);
+ CellData& cellDataI = problem_.variables().cellData(globalIdxI);
+
//for benchmark only!
// problem_.variables().storeSrn(residualSatNW, globalIdxI);
//for benchmark only!
- Scalar porosity =
- problem_.spatialParameters().porosity(*eIt);
-
- Scalar viscosityWI = problem_.variables().viscosityWetting(globalIdxI);
- Scalar viscosityNWI = problem_.variables().viscosityNonwetting(
- globalIdxI);
+ Scalar porosity = problem_.spatialParameters().porosity(*eIt);
- Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
- Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
+ Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
+ Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx);
- Scalar lambdaWI = problem_.variables().mobilityWetting(globalIdxI);
- Scalar lambdaNWI = problem_.variables().mobilityNonwetting(globalIdxI);
+ if (compressibility_)
+ {
+ viscosityW = cellDataI.viscosity(wPhaseIdx);
+ viscosityNW = cellDataI.viscosity(nPhaseIdx);
+ }
// run through all intersections with neighbors and boundary
IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt
- != isItEnd; ++isIt)
+ for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
{
// local number of faces
int isIndex = isIt->indexInInside();
@@ -404,8 +454,10 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
Scalar faceArea = isIt->geometry().volume();
- Scalar factor = 0;
- Scalar factorSecondPhase = 0;
+ //get velocity*normalvector*facearea/(volume*porosity)
+ Scalar factorW = (cellDataI.fluxData().velocity(wPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
+ Scalar factorNW = (cellDataI.fluxData().velocity(nPhaseIdx, isIndex) * unitOuterNormal) * faceArea;
+ Scalar factorTotal = factorW + factorNW;
// handle interior face
if (isIt->neighbor())
@@ -414,82 +466,61 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
ElementPointer neighborPointer = isIt->outside();
int globalIdxJ = problem_.variables().index(*neighborPointer);
+ CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
+
// neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor =
- neighborPointer->geometry().center();
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
// distance vector between barycenters
GlobalPosition distVec = globalPosNeighbor - globalPos;
// compute distance between cell centers
- Scalar dist = distVec * unitOuterNormal;
-// Scalar dist = distVec.two_norm();
+ Scalar dist = distVec.two_norm();
bool takeNeighbor = (eIt->level() < neighborPointer->level());
//get phase potentials
- Scalar potentialW =
- (takeNeighbor) ? -problem_.variables().potentialWetting(globalIdxJ, isIt->indexInOutside()) :
- problem_.variables().potentialWetting(globalIdxI, isIndex);
- Scalar potentialNW =
- (takeNeighbor) ? -problem_.variables().potentialNonwetting(globalIdxJ, isIt->indexInOutside()) :
- problem_.variables().potentialNonwetting(globalIdxI, isIndex);
-
- Scalar densityWJ = problem_.variables().densityWetting(
- globalIdxJ);
- Scalar densityNWJ = problem_.variables().densityNonwetting(
- globalIdxJ);
-
- //get velocity*normalvector*facearea/(volume*porosity)
- factor = (problem_.variables().velocity()[globalIdxI][isIndex]
- * unitOuterNormal) * faceArea;
- factorSecondPhase
- = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex]
- * unitOuterNormal) * faceArea;
+ bool upwindWI =
+ (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(wPhaseIdx, isIt->indexInOutside()) :
+ cellDataI.fluxData().isUpwindCell(wPhaseIdx, isIndex);
+ bool upwindNWI =
+ (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(nPhaseIdx, isIt->indexInOutside()) :
+ cellDataI.fluxData().isUpwindCell(nPhaseIdx, isIndex);
Scalar lambdaW = 0;
Scalar lambdaNW = 0;
- //Scalar viscosityW = 0;
- //Scalar viscosityNW = 0;
//upwinding of lambda dependend on the phase potential gradients
- if (potentialW >= 0.)
+ if (upwindWI)
{
lambdaW = lambdaWI;
if (compressibility_)
{
- lambdaW /= densityWI;
- }//divide by density because lambda is saved as lambda*density
- //viscosityW = viscosityWI;
+ lambdaW /= cellDataI.density(wPhaseIdx);
+ } //divide by density because lambda is saved as lambda*density
}
else
{
- lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
+ lambdaW = cellDataJ.mobility(wPhaseIdx);
if (compressibility_)
{
- lambdaW /= densityWJ;
- }//divide by density because lambda is saved as lambda*density
- //viscosityW = problem_.variables().viscosityWetting(
- // globalIdxJ);
+ lambdaW /= cellDataJ.density(wPhaseIdx);
+ } //divide by density because lambda is saved as lambda*density
}
- if (potentialNW >= 0.)
+ if (upwindNWI)
{
lambdaNW = lambdaNWI;
if (compressibility_)
{
- lambdaNW /= densityNWI;
- }//divide by density because lambda is saved as lambda*density
- //viscosityNW = viscosityNWI;
+ lambdaNW /= cellDataI.density(nPhaseIdx);
+ } //divide by density because lambda is saved as lambda*density
}
else
{
- lambdaNW = problem_.variables().mobilityNonwetting(
- globalIdxJ);
+ lambdaNW = cellDataJ.mobility(nPhaseIdx);
if (compressibility_)
{
- lambdaNW /= densityNWJ;
- }//divide by density because lambda is saved as lambda*density
- //viscosityNW = problem_.variables().viscosityNonwetting(
- // globalIdxJ);
+ lambdaNW /= cellDataJ.density(nPhaseIdx);
+ } //divide by density because lambda is saved as lambda*density
}
switch (velocityType_)
@@ -497,111 +528,68 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
case vt:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorTotal, *isIt);
//determine phase saturations from primary saturation variable
- Scalar satWI = 0;
- Scalar satWJ = 0;
- switch (saturationType_)
- {
- case Sw:
- {
- satWI = problem_.variables().saturation()[globalIdxI];
- satWJ = problem_.variables().saturation()[globalIdxJ];
- break;
- }
- case Sn:
- {
- satWI = 1
- - problem_.variables().saturation()[globalIdxI];
- satWJ = 1
- - problem_.variables().saturation()[globalIdxJ];
- break;
- }
- }
+ Scalar satWI = cellDataI.saturation(wPhaseIdx);
+ Scalar satWJ = cellDataJ.saturation(wPhaseIdx);
- Scalar pcI = problem_.variables().capillaryPressure(
- globalIdxI);
- Scalar pcJ = problem_.variables().capillaryPressure(
- globalIdxJ);
+ Scalar pcI = cellDataI.capillaryPressure();
+ Scalar pcJ = cellDataJ.capillaryPressure();
// calculate the saturation gradient
GlobalPosition pcGradient = unitOuterNormal;
pcGradient *= (pcI - pcJ) / dist;
// get the diffusive part
- Scalar diffPart = diffusivePart()(*eIt, isIndex, satWI,
- satWJ, pcGradient) * unitOuterNormal * faceArea;
+ Scalar diffPart = diffusivePart()(*eIt, isIndex, satWI, satWJ, pcGradient) * unitOuterNormal
+ * faceArea;
-
- Scalar convPart = convectivePart()(*eIt, isIndex, satWI,
- satWJ) * unitOuterNormal * faceArea;
+ Scalar convPart = convectivePart()(*eIt, isIndex, satWI, satWJ) * unitOuterNormal * faceArea;
switch (saturationType_)
{
case Sw:
{
//vt*fw
- factor *= lambdaW / (lambdaW + lambdaNW);
+ factorTotal *= lambdaW / (lambdaW + lambdaNW);
break;
}
case Sn:
{
//vt*fn
- factor *= lambdaNW / (lambdaW + lambdaNW);
+ factorTotal *= lambdaNW / (lambdaW + lambdaNW);
diffPart *= -1;
convPart *= -1;
break;
}
}
- factor -= diffPart;
- factor += convPart;
-
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, 10 * diffPart, *isIt);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, 10 * convPart, *isIt);
-
- break;
- }
- case vw:
- {
- if (compressibility_)
- {
- factor /= densityWI;
- factorSecondPhase /= densityNWI;
- }
+ factorTotal -= diffPart;
+ factorTotal += convPart;
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt, wPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factorSecondPhase,
- *isIt, nPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, 10 * diffPart, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, 10 * convPart, *isIt);
break;
}
- case vn:
+ default:
{
if (compressibility_)
{
- factor /= densityNWI;
- factorSecondPhase /= densityWI;
+ factorW /= cellDataI.density(wPhaseIdx);
+ factorNW /= cellDataI.density(nPhaseIdx);
}
-
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt, nPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factorSecondPhase,
- *isIt, wPhaseIdx);
+ //add cflFlux for time-stepping
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorW, *isIt,
+ wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorNW, *isIt,
+ nPhaseIdx);
break;
}
}
- }//end intersection with neighbor element
+ } //end intersection with neighbor element
// handle boundary face
if (isIt->boundary())
@@ -616,8 +604,7 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
// distance vector between barycenters
GlobalPosition distVec = globalPosFace - globalPos;
// compute distance between cell centers
- Scalar dist = distVec * unitOuterNormal;
-// Scalar dist = distVec.two_norm();
+ Scalar dist = distVec.two_norm();
if (bcType.isDirichlet(satEqIdx))
{
@@ -625,126 +612,70 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
Scalar satBound = boundValues[saturationIdx];
- //get velocity*normalvector*facearea/(volume*porosity)
- factor
- = (problem_.variables().velocity()[globalIdxI][isIndex]
- * unitOuterNormal) * faceArea;
- factorSecondPhase
- = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex]
- * unitOuterNormal) * faceArea;
-
- Scalar pressBound = problem_.variables().pressure()[globalIdxI];
- Scalar temperature = problem_.temperature(*eIt);
-
//determine phase saturations from primary saturation variable
- Scalar satWI = 0;
+ Scalar satWI = cellDataI.saturation(wPhaseIdx);
Scalar satWBound = 0;
switch (saturationType_)
{
case Sw:
{
- satWI = problem_.variables().saturation()[globalIdxI];
satWBound = satBound;
break;
}
case Sn:
{
- satWI = 1
- - problem_.variables().saturation()[globalIdxI];
satWBound = 1 - satBound;
break;
}
}
- Scalar pcBound = MaterialLaw::pC(
- problem_.spatialParameters().materialLawParams(*eIt), satBound);
-
- //determine phase pressures from primary pressure variable
- Scalar pressW = 0;
- Scalar pressNW = 0;
- switch (pressureType_)
- {
- case pw:
- {
- pressW = pressBound;
- pressNW = pressBound + pcBound;
- break;
- }
- case pn:
- {
- pressW = pressBound - pcBound;
- pressNW = pressBound;
- break;
- }
- }
-
- FluidState fluidState;
- fluidState.update(satBound, pressW, pressNW, temperature);
-
- //get phase potentials
- Scalar potentialW = problem_.variables().potentialWetting(
- globalIdxI, isIndex);
- Scalar potentialNW =
- problem_.variables().potentialNonwetting(
- globalIdxI, isIndex);
+ Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt), satBound);
Scalar lambdaW = 0;
Scalar lambdaNW = 0;
//upwinding of lambda dependend on the phase potential gradients
- if (potentialW > 0.)
+ if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, isIndex))
{
lambdaW = lambdaWI;
if (compressibility_)
{
- lambdaW /= densityWI;
- }//divide by density because lambda is saved as lambda*density
+ lambdaW /= cellDataI.density(wPhaseIdx);
+ } //divide by density because lambda is saved as lambda*density
}
else
{
if (compressibility_)
{
- lambdaW
- = MaterialLaw::krw(
- problem_.spatialParameters().materialLawParams(*eIt), satWBound)
- / FluidSystem::phaseViscosity(
- wPhaseIdx, temperature,
- pressW, fluidState);
+ lambdaW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satWBound)
+ / FluidSystem::viscosity(cellDataI.fluidState(), wPhaseIdx);
}
else
{
- lambdaW
- = MaterialLaw::krw(
- problem_.spatialParameters().materialLawParams(*eIt), satWBound)
- / viscosityWI;
+ lambdaW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satWBound)
+ / viscosityW;
}
}
- if (potentialNW >= 0.)
+ if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, isIndex))
{
lambdaNW = lambdaNWI;
if (compressibility_)
{
- lambdaNW /= densityNWI;
- }//divide by density because lambda is saved as lambda*density
+ lambdaNW /= cellDataI.density(nPhaseIdx);
+ } //divide by density because lambda is saved as lambda*density
}
else
{
if (compressibility_)
{
- lambdaNW
- = MaterialLaw::krn(
- problem_.spatialParameters().materialLawParams(*eIt), satWBound)
- / FluidSystem::phaseViscosity(
- nPhaseIdx, temperature,
- pressNW, fluidState);
+ lambdaNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satWBound)
+ / FluidSystem::viscosity(cellDataI.fluidState(), nPhaseIdx);
}
else
{
- lambdaNW
- = MaterialLaw::krn(
- problem_.spatialParameters().materialLawParams(*eIt), satWBound)
- / viscosityNWI;
+ lambdaNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satWBound)
+ / viscosityNW;
}
}
// std::cout<<lambdaW<<" "<<lambdaNW<<std::endl;
@@ -753,136 +684,95 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
{
case vt:
{
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorTotal, *isIt);
- Scalar pcI = problem_.variables().capillaryPressure(
- globalIdxI);
+ Scalar pcI = cellDataI.capillaryPressure();
// calculate the saturation gradient
GlobalPosition pcGradient = unitOuterNormal;
pcGradient *= (pcI - pcBound) / dist;
// get the diffusive part -> give 1-sat because sat = S_n and lambda = lambda(S_w) and pc = pc(S_w)
- Scalar diffPart = diffusivePart()(*eIt, isIndex, satWI,
- satWBound, pcGradient) * unitOuterNormal
+ Scalar diffPart = diffusivePart()(*eIt, isIndex, satWI, satWBound, pcGradient) * unitOuterNormal
* faceArea;
- Scalar convPart = convectivePart()(*eIt, isIndex,
- satWI, satWBound) * unitOuterNormal * faceArea;
+ Scalar convPart = convectivePart()(*eIt, isIndex, satWI, satWBound) * unitOuterNormal
+ * faceArea;
switch (saturationType_)
{
case Sw:
{
//vt*fw
- factor *= lambdaW / (lambdaW + lambdaNW);
+ factorTotal *= lambdaW / (lambdaW + lambdaNW);
break;
}
case Sn:
{
//vt*fn
- factor *= lambdaNW / (lambdaW + lambdaNW);
- diffPart *= -1;
+ factorTotal *= lambdaNW / (lambdaW + lambdaNW);
+ diffPart *= -1; //add cflFlux for time-stepping
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorW, *isIt,
+ wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorNW,
+ *isIt, nPhaseIdx);
convPart *= -1;
break;
}
}
//vt*fw
- factor -= diffPart;
- factor += convPart;
-
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, 10 * diffPart, *isIt);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, 10 * convPart, *isIt);
-
- break;
- }
-
- case vw:
- {
- if (compressibility_)
- {
- factor /= densityWI;
- factorSecondPhase /= densityNWI;
- }
+ factorTotal -= diffPart;
+ factorTotal += convPart;
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt,
- wPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factorSecondPhase,
- *isIt, nPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, 10 * diffPart, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, 10 * convPart, *isIt);
break;
}
-
- //for time step criterion if the non-wetting phase velocity is used
- case vn:
+ default:
{
if (compressibility_)
{
- factor /= densityNWI;
- factorSecondPhase /= densityWI;
+ factorW /= cellDataI.density(wPhaseIdx);
+ factorNW /= cellDataI.density(nPhaseIdx);
}
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt,
- nPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factorSecondPhase,
- *isIt, wPhaseIdx);
+ //add cflFlux for time-stepping
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorW, *isIt,
+ wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorNW, *isIt,
+ nPhaseIdx);
break;
}
}
- }//end dirichlet boundary
+ } //end dirichlet boundary
if (bcType.isNeumann(satEqIdx))
{
problem_.neumann(boundValues, *isIt);
+ factorW = boundValues[wPhaseIdx];
+ factorNW = boundValues[nPhaseIdx];
+ factorW *= faceArea;
+ factorNW *= faceArea;
//get mobilities
Scalar lambdaW, lambdaNW;
lambdaW = lambdaWI;
- if (compressibility_)
- {
- lambdaW /= densityWI;
- }
-
lambdaNW = lambdaNWI;
if (compressibility_)
{
- lambdaNW /= densityNWI;
- }
-
- switch (saturationType_)
- {
- case Sw:
- {
- factor = boundValues[wPhaseIdx];
- factor /= densityWI;
- factor *= faceArea;
- factorSecondPhase = boundValues[nPhaseIdx];
- factorSecondPhase /= densityNWI;
- factorSecondPhase *= faceArea;
- break;
+ lambdaW /= cellDataI.density(wPhaseIdx);
+ lambdaNW /= cellDataI.density(nPhaseIdx);
+ factorW /= cellDataI.density(wPhaseIdx);
+ factorNW /= cellDataI.density(nPhaseIdx);
}
- case Sn:
+ else
{
- factor = boundValues[nPhaseIdx];
- factor /= densityNWI;
- factor *= faceArea;
- factorSecondPhase = boundValues[wPhaseIdx];
- factorSecondPhase /= densityWI;
- factorSecondPhase *= faceArea;
- break;
- }
+ factorW /= densityW;
+ factorNW /= densityNW;
}
switch (velocityType_)
@@ -890,60 +780,33 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
case vt:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor+factorSecondPhase, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorW + factorNW,
+ *isIt);
break;
}
- case vw:
+ default:
{
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt,
- wPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factorSecondPhase, *isIt,
- nPhaseIdx);
- break;
- }
- case vn:
- {
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt,
- nPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factorSecondPhase, *isIt,
- wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorW, *isIt,
+ wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorNW, *isIt,
+ nPhaseIdx);
+
break;
}
}
- }//end neumann boundary
+ } //end neumann boundary
if (bcType.isOutflow(satEqIdx))
{
//get mobilities
- Scalar lambdaW, lambdaNW;
-
- lambdaW = lambdaWI;
+ Scalar lambdaW = lambdaWI;
+ Scalar lambdaNW = lambdaNWI;
if (compressibility_)
{
- lambdaW /= densityWI;
+ lambdaW /= cellDataI.density(wPhaseIdx);
+ lambdaNW /= cellDataI.density(nPhaseIdx);
}
- lambdaNW = lambdaNWI;
- if (compressibility_)
- {
- lambdaNW /= densityNWI;
- }
-
- //get velocity*normalvector*facearea/(volume*porosity)
- factor
- = (problem_.variables().velocity()[globalIdxI][isIndex]
- * unitOuterNormal) * faceArea;
- factorSecondPhase
- = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex]
- * unitOuterNormal) * faceArea;
-
if (velocityType_ == vt)
{
switch (saturationType_)
@@ -951,131 +814,120 @@ void FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
case Sw:
{
//vt*fw
- factor *= lambdaW / (lambdaW + lambdaNW);
+ factorTotal *= lambdaW / (lambdaW + lambdaNW);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorTotal,
+ *isIt);
break;
}
case Sn:
{
//vt*fn
- factor *= lambdaNW / (lambdaW + lambdaNW);
+ factorTotal *= lambdaNW / (lambdaW + lambdaNW);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorTotal,
+ *isIt);
break;
}
}
}
-
- switch (velocityType_)
- {
- case vt:
- {
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt);
- break;
- }
- case vw:
- {
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt,
- wPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factorSecondPhase, *isIt,
- nPhaseIdx);
- break;
- }
- case vn:
+ else
{
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factor, *isIt,
- nPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
- viscosityWI, viscosityNWI, factorSecondPhase, *isIt,
- wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorW, *isIt,
+ wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, factorNW, *isIt,
+ nPhaseIdx);
+
break;
}
- }
- }//end outflow boundary
- }//end boundary
- // add to update vector
- updateVec[globalIdxI] -= factor / (volume * porosity);//-:v>0, if flow leaves the cell
- }// end all intersections
+ } //end outflow boundary
+ } //end boundary
+ // add to update vector
+
+ switch (velocityType_)
+ {
+ case vt:
+ updateVec[globalIdxI] -= factorTotal / (volume * porosity); //-:v>0, if flow leaves the cell
+ break;
+ default:
+ switch (saturationType_)
+ {
+ case Sw:
+ updateVec[globalIdxI] -= factorW / (volume * porosity); //-:v>0, if flow leaves the cell
+ break;
+ case Sn:
+ updateVec[globalIdxI] -= factorNW / (volume * porosity); //-:v>0, if flow leaves the cell
+ break;
+ }
+ break;
+ }
+ } // end all intersections
PrimaryVariables sourceVec(0.0);
problem_.source(sourceVec, *eIt);
- sourceVec[wPhaseIdx] /= densityWI;
- sourceVec[nPhaseIdx] /= densityNWI;
- //get mobilities
- Scalar lambdaW = 0;
- Scalar lambdaNW = 0;
+ if (compressibility_)
+ {
+ sourceVec[wPhaseIdx] /= cellDataI.density(wPhaseIdx);
+ sourceVec[nPhaseIdx] /= cellDataI.density(nPhaseIdx);
+ }
+ else
+ {
+ sourceVec[wPhaseIdx] /= densityW;
+ sourceVec[nPhaseIdx] /= densityNW;
+ }
- lambdaW = lambdaWI;
- if (compressibility_)
- {
- lambdaW /= densityWI;
- }
- lambdaNW = lambdaNWI;
- if (compressibility_)
- {
- lambdaNW /= densityNWI;
- }
+ //get mobilities
+ Scalar lambdaW = lambdaWI;
+ Scalar lambdaNW = lambdaNWI;
+ if (compressibility_)
+ {
+ lambdaW /= cellDataI.density(wPhaseIdx);
+ lambdaNW /= cellDataI.density(nPhaseIdx);
+ }
- switch (saturationType_)
- {
- case Sw:
- {
- if (sourceVec[wPhaseIdx] < 0 && problem_.variables().saturation()[globalIdxI] < threshold_)
- sourceVec[wPhaseIdx] = 0.0;
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ if (sourceVec[wPhaseIdx] < 0 && cellDataI.saturation(wPhaseIdx) < threshold_)
+ sourceVec[wPhaseIdx] = 0.0;
- updateVec[globalIdxI] += sourceVec[wPhaseIdx] / porosity;
- break;
- }
- case Sn:
- {
- if (sourceVec[nPhaseIdx] < 0 && problem_.variables().saturation()[globalIdxI] < threshold_)
- sourceVec[nPhaseIdx] = 0.0;
+ updateVec[globalIdxI] += sourceVec[wPhaseIdx] / porosity;
+ break;
+ }
+ case Sn:
+ {
+ if (sourceVec[nPhaseIdx] < 0 && cellDataI.saturation(nPhaseIdx) < threshold_)
+ sourceVec[nPhaseIdx] = 0.0;
- updateVec[globalIdxI] += sourceVec[nPhaseIdx] / porosity;
- break;
- }
- }
+ updateVec[globalIdxI] += sourceVec[nPhaseIdx] / porosity;
+ break;
+ }
+ }
- switch (velocityType_)
- {
- case vw:
- {
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI,
- viscosityNWI, sourceVec[wPhaseIdx] * volume, *eIt, wPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI,
- viscosityNWI, sourceVec[nPhaseIdx] * volume, *eIt, nPhaseIdx);
- break;
- }
- case vn:
- {
- //add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI,
- viscosityNWI, sourceVec[nPhaseIdx] * volume, *eIt, nPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI,
- viscosityNWI, sourceVec[wPhaseIdx] * volume, *eIt, wPhaseIdx);
- break;
- }
- case vt:
- {
+ switch (velocityType_)
+ {
+ case vt:
+ {
+ //add cflFlux for time-stepping
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW,
+ (sourceVec[wPhaseIdx] + sourceVec[nPhaseIdx]) * volume, *eIt);
+ break;
+ }
+ default:
+ {
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI,
- viscosityNWI, (sourceVec[wPhaseIdx]+sourceVec[nPhaseIdx]) * volume, *eIt);
- break;
- }
- }
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, sourceVec[wPhaseIdx] * volume,
+ *eIt, wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityW, viscosityNW, sourceVec[nPhaseIdx] * volume,
+ *eIt, nPhaseIdx);
+ break;
+ }
+ }
//calculate time step
- dt = std::min(dt, evalCflFluxFunction().getCflFluxFunction(globalPos,
- *eIt) * (porosity * volume));
+ dt = std::min(dt, evalCflFluxFunction().getCflFluxFunction(globalPos, *eIt) * (porosity * volume));
- problem_.variables().volumecorrection(globalIdxI)
- = updateVec[globalIdxI];
+ cellDataI.setUpdate(updateVec[globalIdxI]);
} // end grid traversal
}
@@ -1083,97 +935,80 @@ template<class TypeTag>
void FVSaturation2P<TypeTag>::initialize()
{
// iterate through leaf grid an evaluate c0 at cell center
- ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt
- != eItEnd; ++eIt)
+ ElementIterator eItEnd = problem_.gridView().template end<0>();
+ for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eItEnd; ++eIt)
{
PrimaryVariables initSol(0.0);
problem_.initial(initSol, *eIt);
+
+ int globalIdx = problem_.variables().index(*eIt);
+
// initialize cell concentration
- problem_.variables().saturation()[problem_.variables().index(*eIt)]
- = initSol[saturationIdx];
+ problem_.variables().primaryVariablesGlobal(satEqIdx)[globalIdx] = initSol[saturationIdx];
}
return;
}
template<class TypeTag>
-void FVSaturation2P<TypeTag>::updateMaterialLaws(
- RepresentationType& saturation = *(new RepresentationType(0)),
+void FVSaturation2P<TypeTag>::updateMaterialLaws(RepresentationType& saturation = *(new RepresentationType(0)),
bool iterate = false)
{
- FluidState fluidState;
+ ElementIterator eItBegin = problem_.gridView().template begin<0>();
+
+ Scalar temp = problem_.temperature(*eItBegin);
+ Scalar pRef = problem_.referencePressure(*eItBegin);
+ Scalar densityW = WettingPhase::density(temp, pRef);
+ Scalar densityNW = NonwettingPhase::density(temp, pRef);
+ Scalar viscosityW = WettingPhase::viscosity(temp, pRef);
+ Scalar viscosityNW = NonwettingPhase::viscosity(temp, pRef);
// iterate through leaf grid an evaluate c0 at cell center
- ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt
- != eItEnd; ++eIt)
+ ElementIterator eItEnd = problem_.gridView().template end<0>();
+ for (ElementIterator eIt = eItBegin; eIt != eItEnd; ++eIt)
{
int globalIdx = problem_.variables().index(*eIt);
- Scalar sat = 0;
- if (!iterate)
- {
- sat = problem_.variables().saturation()[globalIdx];
- }
- else
- {
- sat = saturation[globalIdx];
- }
+ CellData& cellData = problem_.variables().cellData(globalIdx);
+
+ Scalar temperature = problem_.temperature(*eIt);
//determine phase saturations from primary saturation variable
Scalar satW = 0;
- if (saturationType_ == Sw)
+ Scalar satNW = 0;
+ switch (saturationType_)
+ {
+ case Sw:
{
- satW = sat;
+ satW = problem_.variables().primaryVariablesGlobal(satEqIdx)[globalIdx];
+ satNW = 1 - satW;
+ break;
}
- if (saturationType_ == Sn)
+ case Sn:
{
- satW = 1 - sat;
+ satNW = problem_.variables().primaryVariablesGlobal(satEqIdx)[globalIdx];
+ satW = 1 - satNW;
+ break;
+ }
}
- Scalar temperature = problem_.temperature(*eIt);
- Scalar referencePressure = problem_.referencePressure(*eIt);
-
- fluidState.update(satW, referencePressure, referencePressure,
- temperature);
-
- problem_.variables().densityWetting(globalIdx)
- = FluidSystem::phaseDensity(wPhaseIdx, temperature,
- referencePressure, fluidState);
- problem_.variables().densityNonwetting(globalIdx)
- = FluidSystem::phaseDensity(nPhaseIdx, temperature,
- referencePressure, fluidState);
- problem_.variables().viscosityWetting(globalIdx)
- = FluidSystem::phaseViscosity(wPhaseIdx, temperature,
- referencePressure, fluidState);
- problem_.variables().viscosityNonwetting(globalIdx)
- = FluidSystem::phaseViscosity(nPhaseIdx, temperature,
- referencePressure, fluidState);
+ Scalar pc = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt), satW);
+
+ cellData.setSaturation(wPhaseIdx, satW);
+ cellData.setSaturation(nPhaseIdx, satNW);
+ cellData.setCapillaryPressure(pc);
// initialize mobilities
- problem_.variables().mobilityWetting(globalIdx)
- = MaterialLaw::krw(
- problem_.spatialParameters().materialLawParams(*eIt), satW)
- / problem_.variables().viscosityWetting(globalIdx);
- problem_.variables().mobilityNonwetting(globalIdx)
- = MaterialLaw::krn(
- problem_.spatialParameters().materialLawParams(*eIt), satW)
- / problem_.variables().viscosityNonwetting(globalIdx);
- problem_.variables().capillaryPressure(globalIdx)
- = MaterialLaw::pC(
- problem_.spatialParameters().materialLawParams(*eIt), satW);
-
- problem_.variables().fracFlowFuncWetting(globalIdx)
- = problem_.variables().mobilityWetting(globalIdx)
- / (problem_.variables().mobilityWetting(globalIdx)
- + problem_.variables().mobilityNonwetting(
- globalIdx));
- problem_.variables().fracFlowFuncNonwetting(globalIdx)
- = problem_.variables().mobilityNonwetting(globalIdx)
- / (problem_.variables().mobilityWetting(globalIdx)
- + problem_.variables().mobilityNonwetting(
- globalIdx));
+ Scalar mobilityW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosityW;
+ Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosityNW;
+
+ // initialize mobilities
+ cellData.mobility(wPhaseIdx) = mobilityW;
+ cellData.mobility(nPhaseIdx) = mobilityNW;
+
+ //initialize fractional flow functions
+ cellData.fracFlowFunc(wPhaseIdx) = mobilityW / (mobilityW + mobilityNW);
+ cellData.fracFlowFunc(nPhaseIdx) = mobilityNW / (mobilityW + mobilityNW);
}
return;
}
diff --git a/dumux/decoupled/2p/transport/fv/gravitypart.hh b/dumux/decoupled/2p/transport/fv/gravitypart.hh
index 58d5dabc007033741b6d59c2bc76db4384b535de..df47bfeaad7dea9f92caf05530e5476b70650894 100644
--- a/dumux/decoupled/2p/transport/fv/gravitypart.hh
+++ b/dumux/decoupled/2p/transport/fv/gravitypart.hh
@@ -52,13 +52,17 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(WettingPhase)) WettingPhase;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NonwettingPhase)) NonwettingPhase;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(CellData)) CellData;
enum
{
@@ -90,8 +94,9 @@ public:
{
Scalar temperature = problem_.temperature(element);
Scalar referencePressure = problem_.referencePressure(element);
- FluidState fluidState;
- fluidState.update(satI, referencePressure, referencePressure, temperature);//not for compressible flow -> thus constant
+
+ Scalar densityW = WettingPhase::density(temperature, referencePressure);
+ Scalar densityNW = NonwettingPhase::density(temperature, referencePressure);
IntersectionIterator isItEnd = problem_.gridView().iend(element);
IntersectionIterator isIt = problem_.gridView().ibegin(element);
@@ -101,38 +106,28 @@ public:
break;
}
int globalIdxI = problem_.variables().index(element);
+ CellData& cellDataI = problem_.variables().cellData(globalIdxI);
// get geometry type of face
//Dune::GeometryType faceGT = isIt->geometryInInside().type();
- Scalar potentialW = problem_.variables().potentialWetting(globalIdxI, indexInInside);
- Scalar potentialNW = problem_.variables().potentialNonwetting(globalIdxI, indexInInside);
-
//get lambda_bar = lambda_n*f_w
Scalar lambdaWI = 0;
Scalar lambdaNWI = 0;
Scalar lambdaWJ = 0;
Scalar lambdaNWJ = 0;
- Scalar densityWI = 0;
- Scalar densityNWI = 0;
- Scalar densityWJ = 0;
- Scalar densityNWJ = 0;
if (preComput_)
{
- lambdaWI=problem_.variables().mobilityWetting(globalIdxI);
- lambdaNWI=problem_.variables().mobilityNonwetting(globalIdxI);
- densityWI = problem_.variables().densityWetting(globalIdxI);
- densityNWI = problem_.variables().densityNonwetting(globalIdxI);
+ lambdaWI=cellDataI.mobility(wPhaseIdx);
+ lambdaNWI=cellDataI.mobility(nPhaseIdx);
}
else
{
lambdaWI = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satI);
- lambdaWI /= FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
+ lambdaWI /= FluidSystem::viscosity(fluidState, wPhaseIdx);
lambdaNWI = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(element), satI);
- lambdaNWI /= FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
- densityWI = FluidSystem::phaseDensity(wPhaseIdx, temperature, referencePressure, fluidState);
- densityNWI = FluidSystem::phaseDensity(nPhaseIdx, temperature, referencePressure, fluidState);
+ lambdaNWI /= FluidSystem::viscosity(fluidState, nPhaseIdx);
}
FieldMatrix meanPermeability(0);
@@ -143,6 +138,7 @@ public:
ElementPointer neighborPointer = isIt->outside();
int globalIdxJ = problem_.variables().index(*neighborPointer);
+ CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
// get permeability
problem_.spatialParameters().meanK(meanPermeability,
@@ -152,19 +148,15 @@ public:
//get lambda_bar = lambda_n*f_w
if (preComput_)
{
- lambdaWJ=problem_.variables().mobilityWetting(globalIdxJ);
- lambdaNWJ=problem_.variables().mobilityNonwetting(globalIdxJ);
- densityWJ = problem_.variables().densityWetting(globalIdxJ);
- densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
+ lambdaWJ=cellDataJ.mobility(wPhaseIdx);
+ lambdaNWJ=cellDataJ.mobility(nPhaseIdx);
}
else
{
lambdaWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*neighborPointer), satJ);
- lambdaWJ /= FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
+ lambdaWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
lambdaNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*neighborPointer), satJ);
- lambdaNWJ /= FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
- densityWJ = FluidSystem::phaseDensity(wPhaseIdx, temperature, referencePressure, fluidState);
- densityNWJ = FluidSystem::phaseDensity(nPhaseIdx, temperature, referencePressure, fluidState);
+ lambdaNWJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
}
}
else
@@ -175,25 +167,22 @@ public:
//calculate lambda_n*f_w at the boundary
lambdaWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satJ);
- lambdaWJ /= FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
+ lambdaWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
lambdaNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(element), satJ);
- lambdaNWJ /= FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
- densityWJ = FluidSystem::phaseDensity(wPhaseIdx, temperature, referencePressure, fluidState);
- densityNWJ = FluidSystem::phaseDensity(nPhaseIdx, temperature, referencePressure, fluidState);
+ lambdaNWJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
}
// set result to K*grad(pc)
FieldVector result(0);
meanPermeability.umv(problem_.gravity(), result);
+ Scalar potentialW = cellDataI.fluxData().potential(wPhaseIdx, indexInInside);
+ Scalar potentialNW = cellDataI.fluxData().potential(nPhaseIdx, indexInInside);
+
Scalar lambdaW = (potentialW >= 0) ? lambdaWI : lambdaWJ;
lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
Scalar lambdaNW = (potentialNW >= 0) ? lambdaNWI : lambdaNWJ;
lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
- Scalar densityW = (potentialW >= 0) ? densityWI : densityWJ;
- densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWJ) : densityW;
- Scalar densityNW = (potentialNW >= 0) ? densityNWI : densityNWJ;
- densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWJ) : densityNW;
// set result to f_w*lambda_n*K*grad(pc)
result *= lambdaW*lambdaNW/(lambdaW+lambdaNW);
diff --git a/dumux/decoupled/2p/transport/transportproblem2p.hh b/dumux/decoupled/2p/transport/transportproblem2p.hh
index d85e51117fe4d7967eb3508593c4d98bc4c8154d..53c4f0bb9606ea2cca37162acd173a0b93847c0a 100644
--- a/dumux/decoupled/2p/transport/transportproblem2p.hh
+++ b/dumux/decoupled/2p/transport/transportproblem2p.hh
@@ -29,8 +29,10 @@
#define DUMUX_TRANSPORTPROBLEM_2P_HH
#include <dumux/decoupled/common/onemodelproblem.hh>
-#include <dumux/decoupled/2p/variableclass2p.hh>
-#include <dumux/material/old_fluidsystems/2p_system.hh>
+#include <dumux/decoupled/common/variableclass.hh>
+#include <dumux/decoupled/2p/cellData2p.hh>
+#include <dumux/material/fluidsystems/2pimmisciblefluidsystem.hh>
+#include <dumux/material/fluidstates/isoimmisciblefluidstate.hh>
#include <dumux/decoupled/2p/2pproperties.hh>
@@ -65,10 +67,16 @@ class TransportProblem2P : public OneModelProblem<TypeTag>
typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
+
enum {
dim = Grid::dimension,
dimWorld = Grid::dimensionworld
};
+ enum
+ {
+ transportEqIdx = Indices::transportEqIdx
+ };
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
@@ -246,7 +254,7 @@ public:
void timeIntegration()
{
// allocate temporary vectors for the updates
- Solution k1 = this->asImp_().variables().saturation();
+ Solution k1 = asImp_().variables().primaryVariablesGlobal(transportEqIdx);
Scalar t = this->timeManager().time();
Scalar dt = 1e100;
@@ -259,7 +267,13 @@ public:
this->timeManager().setTimeStepSize(dt);
// explicit Euler: Sat <- Sat + dt*N(Sat)
- this->asImp_().variables().saturation() += (k1 *= dt);
+ int size = this->gridView().size(0);
+ Solution& newSol = asImp_().variables().primaryVariablesGlobal(transportEqIdx);
+ for (int i = 0; i < size; i++)
+ {
+ newSol[i] += (k1[i] *= dt);
+ this->model().updateSaturationSolution(i, newSol[i][0]);
+ }
}
// \}