From bb9ed0c3696d3753e136cd504c03cc6b315f0415 Mon Sep 17 00:00:00 2001
From: Klaus Mosthaf <klmos@env.dtu.dk>
Date: Sun, 23 Jan 2011 18:05:49 +0000
Subject: [PATCH] reanimated velocity output for the 2p2cmodel - velocity
output can be set in the 2p2clocalresidual - added Kmvp function (without
multiplication by the normal) in 2p2cfluxvariables - direction of velocity
looks good, values not checked yet
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5093 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/2p2c/2p2cfluxvariables.hh | 13 +++++++---
dumux/boxmodels/2p2c/2p2cmodel.hh | 31 +++++++++++++++--------
2 files changed, 31 insertions(+), 13 deletions(-)
diff --git a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
index c48f64207a..7e174493a2 100644
--- a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
@@ -213,7 +213,6 @@ private:
// multiply the pressure potential with the intrinsic
// permeability
Tensor K;
- Vector Kmvp;
for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
{
spatialParams.meanK(K,
@@ -223,8 +222,8 @@ private:
spatialParams.intrinsicPermeability(element,
fvElemGeom_,
face().j));
- K.mv(potentialGrad_[phaseIdx], Kmvp);
- KmvpNormal_[phaseIdx] = - (Kmvp * face().normal);
+ K.mv(potentialGrad_[phaseIdx], Kmvp_[phaseIdx]);
+ KmvpNormal_[phaseIdx] = - (Kmvp_[phaseIdx] * face().normal);
}
// set the upstream and downstream vertices
@@ -292,6 +291,13 @@ public:
Scalar KmvpNormal(int phaseIdx) const
{ return KmvpNormal_[phaseIdx]; }
+ /*!
+ * \brief Return the pressure potential multiplied with the
+ * intrinsic permeability as vector (for velocity output)
+ */
+ Vector Kmvp(int phaseIdx) const
+ { return Kmvp_[phaseIdx]; }
+
/*!
* \brief Return the local index of the upstream control volume
* for a given phase.
@@ -354,6 +360,7 @@ protected:
Scalar densityAtIP_[numPhases], molarDensityAtIP_[numPhases];
// intrinsic permeability times pressure potential gradient
+ Vector Kmvp_[numPhases];
// projected on the face normal
Scalar KmvpNormal_[numPhases];
diff --git a/dumux/boxmodels/2p2c/2p2cmodel.hh b/dumux/boxmodels/2p2c/2p2cmodel.hh
index 6ca7452c30..dae4031fe5 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -108,6 +108,7 @@ class TwoPTwoCModel: public BoxModel<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
@@ -152,6 +153,9 @@ class TwoPTwoCModel: public BoxModel<TypeTag>
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ static const Scalar mobilityUpwindAlpha =
+ GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
+
public:
/*!
* \brief Initialize the static data with the initial solution.
@@ -316,7 +320,7 @@ public:
ScalarField *velocityX = writer.template createField<Scalar, 1>(numVertices);
ScalarField *velocityY = writer.template createField<Scalar, 1>(numVertices);
ScalarField *velocityZ = writer.template createField<Scalar, 1>(numVertices);
- Scalar maxV=0.; // variable to store the maximum face velocity
+// Scalar maxV=0.; // variable to store the maximum face velocity
// initialize velocity fields
Scalar boxSurface[numVertices];
@@ -385,32 +389,39 @@ public:
// In the box method, the velocity is evaluated on the FE-Grid. However, to get an
// average apparent velocity at the vertex, all contributing velocities have to be interpolated.
GlobalPosition velocity(0.);
+ ElementVolumeVariables elemVolVars;
+
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ false /* isOldSol? */);
+
// loop over the phases
- for (int faceIdx = 0; faceIdx< this->fvElemGeom_().numEdges; faceIdx++)
+ for (int faceIdx = 0; faceIdx< fvElemGeom.numEdges; faceIdx++)
{
//prepare the flux calculations (set up and prepare geometry, FE gradients)
FluxVariables fluxDat(this->problem_(),
- this->elem_(),
- this->fvElemGeom_(),
+ *elemIt,
+ fvElemGeom,
faceIdx,
- elemDat);
+ elemVolVars);
// choose phase of interest. Alternatively, a loop over all phases would be possible.
int phaseIdx = gPhaseIdx;
// get darcy velocity
- velocity = fluxDat.KmvpNormal(phaseIdx); // mind the sign: vDarcy = kf grad p
+ velocity = fluxDat.Kmvp(phaseIdx); // mind the sign: vDarcy = kf grad p
// up+downstream mobility
- const VolumeVariables &up = this->curVolVars_(fluxDat.upstreamIdx(phaseIdx));
- const VolumeVariables &down = this->curVolVars_(fluxDat.downstreamIdx(phaseIdx));
+ const VolumeVariables &up = elemVolVars[fluxDat.upstreamIdx(phaseIdx)];
+ const VolumeVariables &down = elemVolVars[fluxDat.downstreamIdx(phaseIdx)];
Scalar scvfArea = fluxDat.face().normal.two_norm(); //get surface area to weight velocity at the IP with the surface area
velocity *= (mobilityUpwindAlpha*up.mobility(phaseIdx) + (1-mobilityUpwindAlpha)*down.mobility(phaseIdx))* scvfArea;
- int vertIIdx = this->problem().vertexMapper().map(this->elem_(),
+ int vertIIdx = this->problem_().vertexMapper().map(*elemIt,
fluxDat.face().i,
dim);
- int vertJIdx = this->problem().vertexMapper().map(this->elem_(),
+ int vertJIdx = this->problem_().vertexMapper().map(*elemIt,
fluxDat.face().j,
dim);
// add surface area for weighting purposes
--
GitLab