From 4fa636e44c809469e03966e15e25c03cdc631688 Mon Sep 17 00:00:00 2001
From: Markus Wolff <markus.wolff@twt-gmbh.de>
Date: Wed, 15 Aug 2012 09:35:56 +0000
Subject: [PATCH] Added initalization function if it did not exist and moved
initialization of constant viscosities and densities from the constructor to
the initialization function
authorized by Benjamin
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8872 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../decoupled/2p/diffusion/fv/fvpressure2p.hh | 34 +++--
.../decoupled/2p/diffusion/fv/fvvelocity2p.hh | 143 ++++++++++--------
.../2p/diffusion/fvmpfa/fvmpfaopressure2p.hh | 25 +--
.../2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh | 24 +--
.../diffusion/mimetic/mimeticgroundwater.hh | 103 ++++++-------
.../2p/diffusion/mimetic/mimeticpressure2p.hh | 26 ++--
.../2p/transport/fv/capillarydiffusion.hh | 1 +
.../2p/transport/fv/convectivepart.hh | 3 +
.../2p/transport/fv/diffusivepart.hh | 4 +
.../decoupled/2p/transport/fv/evalcflflux.hh | 3 +
.../2p/transport/fv/fvsaturation2p.hh | 36 +++--
.../decoupled/2p/transport/fv/gravitypart.hh | 28 ++--
dumux/decoupled/common/fv/fvpressure.hh | 10 ++
dumux/decoupled/common/fv/fvtransport.hh | 6 +-
dumux/decoupled/common/fv/fvvelocity.hh | 5 +
.../decoupled/common/fv/fvvelocitydefault.hh | 3 +
16 files changed, 261 insertions(+), 193 deletions(-)
diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
index 3feaa84025..01bdd17de5 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
@@ -174,6 +174,22 @@ public:
void initialize(bool solveTwice = true)
{
ParentType::initialize();
+
+ if (!compressibility_)
+ {
+ 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);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+ }
+
updateMaterialLaws();
this->assemble(true);
@@ -401,20 +417,10 @@ public:
ErrorTermLowerBound_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, ErrorTermLowerBound);
ErrorTermUpperBound_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, ErrorTermUpperBound);
- if (!compressibility_)
- {
- 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);
- viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
- viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
- }
+ density_[wPhaseIdx] = 0.;
+ density_[nPhaseIdx] = 0.;
+ viscosity_[wPhaseIdx] = 0.;
+ viscosity_[nPhaseIdx] = 0.;
}
private:
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index 5a1cd7804d..5b061207e9 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -57,7 +57,7 @@ namespace Dumux
template<class TypeTag>
class FVVelocity2P
{
- 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;
@@ -115,7 +115,7 @@ public:
* \param problem A Problem class object
*/
FVVelocity2P(Problem& problem) :
- problem_(problem), gravity_(problem.gravity())
+ problem_(problem), gravity_(problem.gravity())
{
if (GET_PROP_VALUE(TypeTag, EnableCompressibility) && velocityType_ == vt)
{
@@ -127,6 +127,15 @@ public:
DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
}
+ density_[wPhaseIdx] = 0.;
+ density_[nPhaseIdx] = 0.;
+ viscosity_[wPhaseIdx] = 0.;
+ viscosity_[nPhaseIdx] = 0.;
+ }
+
+ //! For initialization
+ void initialize()
+ {
if (!compressibility_)
{
ElementIterator element = problem_.gridView().template begin<0> ();
@@ -172,7 +181,7 @@ public:
Dune::BlockVector < Dune::FieldVector<Scalar, dim> > &velocity = *(writer.template allocateManagedBuffer<Scalar,
dim>(problem_.gridView().size(0)));
Dune::BlockVector < Dune::FieldVector<Scalar, dim> > &velocitySecondPhase =
- *(writer.template allocateManagedBuffer<Scalar, dim>(problem_.gridView().size(0)));
+ *(writer.template allocateManagedBuffer<Scalar, dim>(problem_.gridView().size(0)));
// compute update vector
ElementIterator eItEnd = problem_.gridView().template end<0>();
@@ -192,20 +201,20 @@ public:
int isIndex = isIt->indexInInside();
fluxW[isIndex] += isIt->geometry().volume()
- * (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
+ * (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNW[isIndex] += isIt->geometry().volume()
- * (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
+ * (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
}
Dune::FieldVector < Scalar, dim > refVelocity(0);
- refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
+ refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
if (dim > 1)
- refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
+ refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
if (dim == 3)
- refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
+ refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
const Dune::FieldVector<Scalar, dim>& localPos =
- ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
+ ReferenceElementContainer::general(eIt->geometry().type()).position(0, 0);
// get the transposed Jacobian of the element mapping
const DimMatrix& jacobianInv = eIt->geometry().jacobianInverseTransposed(localPos);
@@ -224,7 +233,7 @@ public:
if (dim > 1)
refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
if (dim == 3)
- refVelocity[2] = 0.5 * (fluxNW[5] - fluxNW[4]);
+ refVelocity[2] = 0.5 * (fluxNW[5] - fluxNW[4]);
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
@@ -237,23 +246,23 @@ public:
//switch velocities
switch (velocityType_)
{
- case vw:
- {
- writer.attachCellData(velocity, "wetting-velocity", dim);
- writer.attachCellData(velocitySecondPhase, "non-wetting-velocity", dim);
- break;
- }
- case vn:
- {
- writer.attachCellData(velocity, "non-wetting-velocity", dim);
- writer.attachCellData(velocitySecondPhase, "wetting-velocity", dim);
- break;
- }
- case vt:
- {
- writer.attachCellData(velocity, "total velocity", dim);
- break;
- }
+ case vw:
+ {
+ writer.attachCellData(velocity, "wetting-velocity", dim);
+ writer.attachCellData(velocitySecondPhase, "non-wetting-velocity", dim);
+ break;
+ }
+ case vn:
+ {
+ writer.attachCellData(velocity, "non-wetting-velocity", dim);
+ writer.attachCellData(velocitySecondPhase, "wetting-velocity", dim);
+ break;
+ }
+ case vt:
+ {
+ writer.attachCellData(velocity, "total velocity", dim);
+ break;
+ }
}
return;
@@ -265,19 +274,19 @@ private:
Scalar density_[numPhases];
Scalar viscosity_[numPhases];
- static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation); //!< gives kind of velocity used (\f$ 0 = v_w\f$, \f$ 1 = v_n\f$, \f$ 2 = v_t\f$)
+ static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation);//!< gives kind of velocity used (\f$ 0 = v_w\f$, \f$ 1 = v_n\f$, \f$ 2 = v_t\f$)
static const bool compressibility_ = GET_PROP_VALUE(TypeTag, EnableCompressibility);
- static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$p_w\f$, \f$p_n\f$, \f$p_{global}\f$)
- static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation); //!< gives kind of saturation used (\f$S_w\f$, \f$S_n\f$)
+ static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation);//!< gives kind of pressure used (\f$p_w\f$, \f$p_n\f$, \f$p_{global}\f$)
+ static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);//!< gives kind of saturation used (\f$S_w\f$, \f$S_n\f$)
};
/*! \brief Calculates the velocity at a cell-cell interface.
-*
-* Calculates the velocity at a cell-cell interface from a given pressure field.
-*
-* \param intersection Intersection of two grid cells
-* \param cellData Object containing all model relevant cell data
-*/
+ *
+ * Calculates the velocity at a cell-cell interface from a given pressure field.
+ *
+ * \param intersection Intersection of two grid cells
+ * \param cellData Object containing all model relevant cell data
+ */
template<class TypeTag>
void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
{
@@ -325,7 +334,7 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*elementI),
problem_.spatialParams().intrinsicPermeability(*elementJ));
- Dune::FieldVector < Scalar, dim > permeability(0);
+ Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
//calculate potential gradients
@@ -393,14 +402,14 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
Scalar scalarPerm = permeability.two_norm();
//calculate the gravity term
- Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
- Dune::FieldVector < Scalar, dimWorld > velocityNW(unitOuterNormal);
+ Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
+ Dune::FieldVector<Scalar, dimWorld> velocityNW(unitOuterNormal);
//calculate unit distVec
distVec /= dist;
Scalar areaScaling = (unitOuterNormal * distVec);
//this treatment of g allows to account for gravity flux through faces where the face normal has no z component (e.g. parallelepiped grids)
- Scalar gravityTermW = (gravity_* distVec) * density_[wPhaseIdx] * areaScaling;
+ Scalar gravityTermW = (gravity_ * distVec) * density_[wPhaseIdx] * areaScaling;
Scalar gravityTermNW = (gravity_ * distVec) * density_[nPhaseIdx] * areaScaling;
//calculate velocity depending on the pressure used -> use pc = pn - pw
@@ -408,22 +417,26 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
{
case pw:
{
- velocityW *= lambdaW * scalarPerm * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermW);
- velocityNW *= lambdaNW * scalarPerm * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermNW)
+ velocityW *= lambdaW * scalarPerm
+ * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermW);
+ velocityNW *= lambdaNW * scalarPerm
+ * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermNW)
+ 0.5 * (lambdaNWI + lambdaNWJ) * scalarPerm * (pcI - pcJ) / dist;
break;
}
case pn:
{
- velocityW *= lambdaW * scalarPerm * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermW)
+ velocityW *= lambdaW * scalarPerm
+ * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermW)
- 0.5 * (lambdaWI + lambdaWJ) * scalarPerm * (pcI - pcJ) / dist;
- velocityNW *= lambdaNW * scalarPerm * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermNW);
+ velocityNW *= lambdaNW * scalarPerm
+ * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermNW);
break;
}
case pglobal:
{
- velocityW *= (lambdaW + lambdaNW) * scalarPerm * (cellData.globalPressure() - cellDataJ.globalPressure()) / dist +
- scalarPerm * (lambdaW * gravityTermW + lambdaNW * gravityTermNW);
+ velocityW *= (lambdaW + lambdaNW) * scalarPerm * (cellData.globalPressure() - cellDataJ.globalPressure()) / dist
+ + scalarPerm * (lambdaW * gravityTermW + lambdaNW * gravityTermNW);
velocityNW = 0;
break;
}
@@ -443,12 +456,12 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
}
/*! \brief Calculates the velocity at a boundary.
-*
-* Calculates the velocity at a boundary from a given pressure field.
-*
-* \param intersection Boundary intersection
-* \param cellData Object containing all model relevant cell data
-*/
+ *
+ * Calculates the velocity at a boundary from a given pressure field.
+ *
+ * \param intersection Boundary intersection
+ * \param cellData Object containing all model relevant cell data
+ */
template<class TypeTag>
void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
{
@@ -496,7 +509,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*element));
- Dune::FieldVector < Scalar, dim > permeability(0);
+ Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
//determine saturation at the boundary -> if no saturation is known directly at the boundary use the cell saturation
@@ -580,9 +593,11 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
if (compressibility_)
{
density_[wPhaseIdx] = (potentialW > 0.) ? cellData.density(wPhaseIdx) : densityWBound;
- density_[wPhaseIdx] = (potentialW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
+ density_[wPhaseIdx] =
+ (potentialW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
density_[nPhaseIdx] = (potentialNW > 0.) ? cellData.density(nPhaseIdx) : densityNWBound;
- density_[nPhaseIdx] = (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx];
+ density_[nPhaseIdx] =
+ (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx];
}
//calculate potential gradient
@@ -613,16 +628,18 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
if (compressibility_)
{
density_[wPhaseIdx] = (potentialW > 0.) ? cellData.density(wPhaseIdx) : densityWBound;
- density_[wPhaseIdx] = (potentialW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
+ density_[wPhaseIdx] =
+ (potentialW == 0) ? 0.5 * (cellData.density(wPhaseIdx) + densityWBound) : density_[wPhaseIdx];
density_[nPhaseIdx] = (potentialNW > 0.) ? cellData.density(nPhaseIdx) : densityNWBound;
- density_[nPhaseIdx] = (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx];
+ density_[nPhaseIdx] =
+ (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx];
}
Scalar scalarPerm = permeability.two_norm();
//calculate the gravity term
- Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
- Dune::FieldVector < Scalar, dimWorld > velocityNW(unitOuterNormal);
+ Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
+ Dune::FieldVector<Scalar, dimWorld> velocityNW(unitOuterNormal);
//calculate unit distVec
distVec /= dist;
@@ -650,8 +667,8 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
}
case pglobal:
{
- velocityW *= (lambdaW + lambdaNW) * scalarPerm * (cellData.globalPressure() - pressBound) / dist +
- scalarPerm * (lambdaW * gravityTermW + lambdaNW * gravityTermNW);
+ velocityW *= (lambdaW + lambdaNW) * scalarPerm * (cellData.globalPressure() - pressBound) / dist
+ + scalarPerm * (lambdaW * gravityTermW + lambdaNW * gravityTermNW);
velocityNW = 0;
break;
}
@@ -668,8 +685,8 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
{
problem_.neumann(boundValues, intersection);
- Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
- Dune::FieldVector < Scalar, dimWorld > velocityNW(unitOuterNormal);
+ Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
+ Dune::FieldVector<Scalar, dimWorld> velocityNW(unitOuterNormal);
velocityW *= boundValues[wPhaseIdx];
velocityNW *= boundValues[nPhaseIdx];
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh
index 72678fc842..670b563f18 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh
@@ -154,6 +154,19 @@ public:
void initialize(bool solveTwice = true)
{
ParentType::initialize();
+
+ const Element& 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);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+
updateMaterialLaws();
assemble();
@@ -248,18 +261,6 @@ public:
{
DUNE_THROW(Dune::NotImplemented, "MPFA method only implemented for 2-d!");
}
-
- const Element& 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);
- viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
- viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
}
private:
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh
index fce17a7732..7d5bc2a0dc 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh
@@ -125,18 +125,6 @@ public:
{
DUNE_THROW(Dune::NotImplemented, "MPFA method only implemented for 2-d!");
}
-
- const Element& 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);
- viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
- viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
}
//Calculates the velocities at all cell-cell interfaces.
@@ -150,6 +138,18 @@ public:
{
ParentType::initialize(solveTwice);
+ const Element& 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);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+
calculateVelocity();
return;
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
index c637f77ce8..bd9b23c946 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
@@ -1,11 +1,11 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
-* Copyright (C) 2008-2009 by Bernd Flemisch *
-* Institute for Modelling Hydraulic and Environmental Systems *
-* University of Stuttgart, Germany *
-* email: <givenname>.<name>@iws.uni-stuttgart.de *
-* *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * 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 *
@@ -18,7 +18,7 @@
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
-*****************************************************************************/
+ *****************************************************************************/
/*!
* \file
*
@@ -44,12 +44,12 @@
namespace Dumux
{
/*!
-* \ingroup MimeticPressure2p
-*/
+ * \ingroup MimeticPressure2p
+ */
/**
-* @brief compute local stiffness matrix for conforming finite elements for diffusion equation
-*
-*/
+ * @brief compute local stiffness matrix for conforming finite elements for diffusion equation
+ *
+ */
//! A class for computing local stiffness matrices
/*! A class for computing local stiffness matrix for the diffusion equation
@@ -62,11 +62,9 @@ namespace Dumux
* \tparam TypeTag The problem TypeTag
*/
template<class TypeTag>
-class MimeticGroundwaterEquationLocalStiffness
-:
-public LocalStiffness<TypeTag, 1>
+class MimeticGroundwaterEquationLocalStiffness: public LocalStiffness<TypeTag, 1>
{
- 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;
@@ -113,6 +111,10 @@ public:
bool levelBoundaryAsDirichlet, const GridView& gridView,
bool procBoundaryAsDirichlet=true)
: problem_(problem), gridView_(gridView)
+ {}
+
+ //! For initialization
+ void initialize()
{
ElementIterator element = problem_.gridView().template begin<0> ();
FluidState fluidState;
@@ -127,14 +129,14 @@ public:
//! assemble local stiffness matrix for given element and order
/*! On exit the following things have been done:
- - The stiffness matrix for the given entity and polynomial degree has been assembled and is
- accessible with the mat() method.
- - The boundary conditions have been evaluated and are accessible with the bc() method
- - The right hand side has been assembled. It contains either the value of the essential boundary
- condition or the assembled source term and neumann boundary condition. It is accessible via the rhs() method.
- @param[in] element a codim 0 entity reference
- @param[in] k order of CR basis (only k = 1 is implemented)
- */
+ - The stiffness matrix for the given entity and polynomial degree has been assembled and is
+ accessible with the mat() method.
+ - The boundary conditions have been evaluated and are accessible with the bc() method
+ - The right hand side has been assembled. It contains either the value of the essential boundary
+ condition or the assembled source term and neumann boundary condition. It is accessible via the rhs() method.
+ @param[in] element a codim 0 entity reference
+ @param[in] k order of CR basis (only k = 1 is implemented)
+ */
void assemble (const Element& element, int k=1)
{
unsigned int numFaces = element.template count<1>();
@@ -146,7 +148,7 @@ public:
this->b[i] = 0;
this->bctype[i].reset();
for (unsigned int j=0; j<numFaces; j++)
- this->A[i][j] = 0;
+ this->A[i][j] = 0;
}
assembleV(element,k);
@@ -165,12 +167,12 @@ public:
//! assemble only boundary conditions for given element
/*! On exit the following things have been done:
- - The boundary conditions have been evaluated and are accessible with the bc() method
- - The right hand side contains either the value of the essential boundary
- condition or the assembled neumann boundary condition. It is accessible via the rhs() method.
- @param[in] element a codim 0 entity reference
- @param[in] k order of CR basis
- */
+ - The boundary conditions have been evaluated and are accessible with the bc() method
+ - The right hand side contains either the value of the essential boundary
+ condition or the assembled neumann boundary condition. It is accessible via the rhs() method.
+ @param[in] element a codim 0 entity reference
+ @param[in] k order of CR basis
+ */
void assembleBoundaryCondition (const Element& element, int k=1)
{
unsigned int numFaces = element.template count<1>();
@@ -232,8 +234,8 @@ public:
N[i] = unitOuterNormal;
for (int k = 0; k < dim; k++)
- // move origin to the center of gravity
- R[i][k] = faceVol[i]*(faceGlobal[k] - centerGlobal[k]);
+ // move origin to the center of gravity
+ R[i][k] = faceVol[i]*(faceGlobal[k] - centerGlobal[k]);
}
// std::cout << "N =\dim" << N;
@@ -244,21 +246,21 @@ public:
// (1) orthonormalize columns of the matrix R
Scalar norm = R[0][0]*R[0][0];
for (unsigned int i = 1; i < numFaces; i++)
- norm += R[i][0]*R[i][0];
+ norm += R[i][0]*R[i][0];
norm = sqrt(norm);
for (unsigned int i = 0; i < numFaces; i++)
- R[i][0] /= norm;
+ R[i][0] /= norm;
Scalar weight = R[0][1]*R[0][0];
for (unsigned int i = 1; i < numFaces; i++)
- weight += R[i][1]*R[i][0];
+ weight += R[i][1]*R[i][0];
for (unsigned int i = 0; i < numFaces; i++)
- R[i][1] -= weight*R[i][0];
+ R[i][1] -= weight*R[i][0];
norm = R[0][1]*R[0][1];
for (unsigned int i = 1; i < numFaces; i++)
- norm += R[i][1]*R[i][1];
+ norm += R[i][1]*R[i][1];
norm = sqrt(norm);
for (unsigned int i = 0; i < numFaces; i++)
- R[i][1] /= norm;
+ R[i][1] /= norm;
if (dim == 3)
{
Scalar weight1 = R[0][2]*R[0][0];
@@ -269,13 +271,13 @@ public:
weight2 += R[i][2]*R[i][1];
}
for (unsigned int i = 0; i < numFaces; i++)
- R[i][1] -= weight1*R[i][0] + weight2*R[i][1];
+ R[i][1] -= weight1*R[i][0] + weight2*R[i][1];
norm = R[0][2]*R[0][2];
for (unsigned int i = 1; i < numFaces; i++)
- norm += R[i][2]*R[i][2];
+ norm += R[i][2]*R[i][2];
norm = sqrt(norm);
for (unsigned int i = 0; i < numFaces; i++)
- R[i][2] /= norm;
+ R[i][2] /= norm;
}
// std::cout << "~R =\dim" << R;
@@ -297,7 +299,7 @@ public:
Scalar traceK = K[0][0];
for (int i = 1; i < dim; i++)
- traceK += K[i][i];
+ traceK += K[i][i];
D *= 2.0*traceK/volume;
// std::cout << "u~D =\dim" << D;
@@ -310,7 +312,7 @@ public:
{
W[i][j] = NK[i][0]*N[j][0];
for (int k = 1; k < dim; k++)
- W[i][j] += NK[i][k]*N[j][k];
+ W[i][j] += NK[i][k]*N[j][k];
}
}
@@ -322,7 +324,6 @@ public:
// std::cout << D[3][2] << ", " << D[3][0] << ", " << D[3][3] << ", " << D[3][1] << std::endl;
// std::cout << D[1][2] << ", " << D[1][0] << ", " << D[1][3] << ", " << D[1][1] << std::endl;
-
// Now the notation is borrowed from Aarnes/Krogstadt/Lie 2006, Section 3.4.
// The matrix W developed so far corresponds to one element-associated
// block of the matrix B^{-1} there.
@@ -332,12 +333,12 @@ public:
// This is just a row vector of size numFaces.
// scale with volume
for (unsigned int i = 0; i < numFaces; i++)
- c[i] = faceVol[i];
+ c[i] = faceVol[i];
// Set up the element part of the matrix \Pi coupling velocities
// and pressure-traces. This is a diagonal matrix with entries given by faceVol.
for (unsigned int i = 0; i < numFaces; i++)
- Pi[i][i] = faceVol[i];
+ Pi[i][i] = faceVol[i];
// Calculate the element part of the matrix D^{-1} = (c W c^T)^{-1} which is just a scalar value.
Dune::FieldVector<Scalar,2*dim> Wc(0);
@@ -380,19 +381,19 @@ private:
// Calculate the element part of the matrix F D^{-1} F^T.
Dune::FieldMatrix<Scalar,2*dim,2*dim> FDinvFT(0);
for (unsigned int i = 0; i < numFaces; i++)
- for (unsigned int j = 0; j < numFaces; j++)
- FDinvFT[i][j] = dinv*F[i]*F[j];
+ for (unsigned int j = 0; j < numFaces; j++)
+ FDinvFT[i][j] = dinv*F[i]*F[j];
// Calculate the element part of the matrix S = Pi W Pi^T - F D^{-1} F^T.
for (unsigned int i = 0; i < numFaces; i++)
- for (unsigned int j = 0; j < numFaces; j++)
- this->A[i][j] = PiWPiT[i][j] - FDinvFT[i][j];
+ for (unsigned int j = 0; j < numFaces; j++)
+ this->A[i][j] = PiWPiT[i][j] - FDinvFT[i][j];
// Calculate the source term F D^{-1} f
// NOT WORKING AT THE MOMENT
Scalar factor = dinv*qmean;
for (unsigned int i = 0; i < numFaces; i++)
- this->b[i] = F[i]*factor;
+ this->b[i] = F[i]*factor;
// std::cout << "faceVol = " << faceVol << std::endl << "W = " << std::endl << W << std::endl
// << "c = " << c << std::endl << "Pi = " << std::endl << Pi << std::endl
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
index 984ffae108..a9a55bef7c 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh
@@ -128,6 +128,7 @@ template<class TypeTag> class MimeticPressure2P
void assemble(bool first)
{
LocalStiffness lstiff(problem_, false, problem_.gridView());
+ lstiff.initialize();
A_.assemble(lstiff, pressTrace_, f_);
return;
}
@@ -138,6 +139,7 @@ template<class TypeTag> class MimeticPressure2P
void postprocess()
{
LocalStiffness lstiff(problem_, false, problem_.gridView());
+ lstiff.initialize();
A_.calculatePressure(lstiff, pressTrace_, normalVelocity_, pressure_);
return;
}
@@ -168,6 +170,18 @@ public:
*/
void initialize(bool solveTwice = true)
{
+ 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);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+
updateMaterialLaws();
A_.initializeMatrix();
f_.resize(problem_.gridView().size(1));//resize to make sure the final grid size (after the problem was completely built) is used!
@@ -289,18 +303,6 @@ public:
{
DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
}
-
- 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);
- viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
- viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
}
private:
diff --git a/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh b/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
index f2d863f5da..dc215c1515 100644
--- a/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
+++ b/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
@@ -81,6 +81,7 @@ private:
typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix;
public:
+
//! Returns capillary diffusion term
/*! Returns capillary diffusion term for current element face
* \param flux Flux vector (gets the flux from the function)
diff --git a/dumux/decoupled/2p/transport/fv/convectivepart.hh b/dumux/decoupled/2p/transport/fv/convectivepart.hh
index 2572f07ff8..869fb2b8f7 100644
--- a/dumux/decoupled/2p/transport/fv/convectivepart.hh
+++ b/dumux/decoupled/2p/transport/fv/convectivepart.hh
@@ -52,6 +52,9 @@ private:
typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
public:
+ //! For initialization
+ void initialize(){};
+
/*! \brief Returns convective term for current element face
* \param intersection Intersection of two grid elements/global boundary
* \param sat Saturation of current element
diff --git a/dumux/decoupled/2p/transport/fv/diffusivepart.hh b/dumux/decoupled/2p/transport/fv/diffusivepart.hh
index 7dee5ba53c..6d7e7f8ee8 100644
--- a/dumux/decoupled/2p/transport/fv/diffusivepart.hh
+++ b/dumux/decoupled/2p/transport/fv/diffusivepart.hh
@@ -51,6 +51,10 @@ private:
typedef Dune::FieldVector<Scalar, dim> DimVector;
public:
+
+ //! For initialization
+ void initialize(){};
+
/*! \brief Returns diffusive term for current element face
*
* \param flux Flux vector (gets the flux from the function)
diff --git a/dumux/decoupled/2p/transport/fv/evalcflflux.hh b/dumux/decoupled/2p/transport/fv/evalcflflux.hh
index 2ac478bc72..74ffcb011a 100644
--- a/dumux/decoupled/2p/transport/fv/evalcflflux.hh
+++ b/dumux/decoupled/2p/transport/fv/evalcflflux.hh
@@ -69,6 +69,9 @@ private:
public:
+ //! For initialization
+ void initialize(){};
+
/*! \brief adds a flux to the cfl-criterion evaluation
*
* \param lambdaW wetting phase mobility
diff --git a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
index 985caf365c..a32297ecd7 100644
--- a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
+++ b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
@@ -371,21 +371,6 @@ public:
capillaryFlux_ = new CapillaryFlux(problem);
gravityFlux_ = new GravityFlux(problem);
velocity_ = new Velocity(problem);
-
- if (!compressibility_)
- {
- 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);
- viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
- viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
- }
}
//! Destructor
@@ -984,6 +969,23 @@ void FVSaturation2P<TypeTag>::getSource(Scalar& update, const Element& element,
template<class TypeTag>
void FVSaturation2P<TypeTag>::initialize()
{
+ ParentType::initialize();
+
+ if (!compressibility_)
+ {
+ 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);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+ }
+
// 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)
@@ -1013,6 +1015,10 @@ void FVSaturation2P<TypeTag>::initialize()
}
}
+ velocity_->initialize();
+ capillaryFlux_->initialize();
+ gravityFlux_->initialize();
+
return;
}
diff --git a/dumux/decoupled/2p/transport/fv/gravitypart.hh b/dumux/decoupled/2p/transport/fv/gravitypart.hh
index 8de8f7011f..cc0097e970 100644
--- a/dumux/decoupled/2p/transport/fv/gravitypart.hh
+++ b/dumux/decoupled/2p/transport/fv/gravitypart.hh
@@ -51,18 +51,18 @@ template<class TypeTag>
class GravityPart: public ConvectivePart<TypeTag>
{
private:
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView)GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
- typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
- typedef typename SpatialParams::MaterialLaw MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+ typedef typename SpatialParams::MaterialLaw MaterialLaw;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
- typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
enum
{
@@ -82,6 +82,7 @@ private:
typedef Dune::FieldMatrix<Scalar,dim,dim> DimMatrix;
public:
+
/*! \brief Returns convective term for current element face
*
* \param flux Flux vector (gets the flux from the function)
@@ -180,6 +181,10 @@ public:
*/
GravityPart (Problem& problem)
: ConvectivePart<TypeTag>(problem), problem_(problem), preComput_(GET_PROP_VALUE(TypeTag, PrecomputedConstRels))
+ {}
+
+ //! For initialization
+ void initialize()
{
ElementIterator element = problem_.gridView().template begin<0> ();
FluidState fluidState;
@@ -195,11 +200,10 @@ public:
}
private:
- Problem& problem_;//problem data
+ Problem& problem_; //problem data
const bool preComput_;//if preCompute = true the mobilities are taken from the variable object, if preCompute = false new mobilities will be taken (for implicit Scheme)
Scalar density_[numPhases];
Scalar viscosity_[numPhases];
-};
-}
+};}
#endif
diff --git a/dumux/decoupled/common/fv/fvpressure.hh b/dumux/decoupled/common/fv/fvpressure.hh
index 667c540856..f4812c7049 100644
--- a/dumux/decoupled/common/fv/fvpressure.hh
+++ b/dumux/decoupled/common/fv/fvpressure.hh
@@ -173,6 +173,16 @@ public:
const Scalar pressure(const int globalIdx) const
{ return pressure_[globalIdx];}
+ const Matrix& globalMatrix()
+ {
+ return A_;
+ }
+
+ const RHSVector& rightHandSide()
+ {
+ return f_;
+ }
+
/*! \brief Initialize pressure model
*
* Function initializes the sparse matrix to solve the global system of equations and sets/calculates the initial pressure
diff --git a/dumux/decoupled/common/fv/fvtransport.hh b/dumux/decoupled/common/fv/fvtransport.hh
index 80e5fffb09..55eda38727 100644
--- a/dumux/decoupled/common/fv/fvtransport.hh
+++ b/dumux/decoupled/common/fv/fvtransport.hh
@@ -123,8 +123,10 @@ public:
void getSource(Scalar& update, const Element& element, CellData& cellDataI);
//! Sets the initial solution \f$ S_0 \f$.
- void initialize();
-
+ void initialize()
+ {
+ evalCflFluxFunction_->initialize();
+ }
/*! \brief Updates constitutive relations and stores them in the variable class*/
void updateMaterialLaws();
diff --git a/dumux/decoupled/common/fv/fvvelocity.hh b/dumux/decoupled/common/fv/fvvelocity.hh
index 2fda6715bf..5651918dc4 100644
--- a/dumux/decoupled/common/fv/fvvelocity.hh
+++ b/dumux/decoupled/common/fv/fvvelocity.hh
@@ -69,6 +69,11 @@ template<class TypeTag, class Velocity> class FVVelocity
public:
+ void initialize()
+ {
+ velocity_.initialize();
+ }
+
//function which iterates through the grid and calculates the global velocity field
void calculateVelocity();
diff --git a/dumux/decoupled/common/fv/fvvelocitydefault.hh b/dumux/decoupled/common/fv/fvvelocitydefault.hh
index 041a6fc43f..76bee6510e 100644
--- a/dumux/decoupled/common/fv/fvvelocitydefault.hh
+++ b/dumux/decoupled/common/fv/fvvelocitydefault.hh
@@ -59,6 +59,9 @@ public:
FVVelocityDefault(Problem& problem)
{}
+ //! For initialization
+ void initialize(){};
+
//! Local velocity calculation
void calculateVelocity(const Intersection& intersection, CellData& cellData)
{}
--
GitLab