From 54b2641014bc17bd8f567326d13621c6745d0f58 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 22 May 2013 14:16:51 +0000
Subject: [PATCH] fix implicit velocity output, which counted too many
surrounding cells in case of more than one phase. This also made the function
completeVelocityCalculation obsolete.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10733 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/1p/1pmodel.hh | 2 -
dumux/implicit/1p2c/1p2cmodel.hh | 2 -
dumux/implicit/2p/2pmodel.hh | 3 -
dumux/implicit/2p2c/2p2cmodel.hh | 3 -
dumux/implicit/3p3c/3p3cmodel.hh | 4 --
.../implicit/common/implicitvelocityoutput.hh | 64 ++++++++-----------
dumux/implicit/mpnc/mpncvtkwritercommon.hh | 7 +-
dumux/implicit/richards/richardsmodel.hh | 2 -
8 files changed, 30 insertions(+), 57 deletions(-)
diff --git a/dumux/implicit/1p/1pmodel.hh b/dumux/implicit/1p/1pmodel.hh
index 132ab64234..8e5ad1e622 100644
--- a/dumux/implicit/1p/1pmodel.hh
+++ b/dumux/implicit/1p/1pmodel.hh
@@ -134,8 +134,6 @@ public:
velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *elemIt, /*phaseIdx=*/0);
}
- velocityOutput.completeVelocityCalculation(*velocity);
-
writer.attachDofData(*p, "p", isBox);
writer.attachDofData(*K, "K", isBox);
if (velocityOutput.enableOutput())
diff --git a/dumux/implicit/1p2c/1p2cmodel.hh b/dumux/implicit/1p2c/1p2cmodel.hh
index 110070e540..be916cee7c 100644
--- a/dumux/implicit/1p2c/1p2cmodel.hh
+++ b/dumux/implicit/1p2c/1p2cmodel.hh
@@ -175,8 +175,6 @@ public:
velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *elemIt, phaseIdx);
}
- velocityOutput.completeVelocityCalculation(*velocity);
-
writer.attachDofData(pressure, "P", isBox);
writer.attachDofData(delp, "delp", isBox);
if (velocityOutput.enableOutput())
diff --git a/dumux/implicit/2p/2pmodel.hh b/dumux/implicit/2p/2pmodel.hh
index 3f39e85688..3d64da3563 100644
--- a/dumux/implicit/2p/2pmodel.hh
+++ b/dumux/implicit/2p/2pmodel.hh
@@ -187,9 +187,6 @@ public:
velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *elemIt, nPhaseIdx);
}
- velocityOutput.completeVelocityCalculation(*velocityW);
- velocityOutput.completeVelocityCalculation(*velocityN);
-
writer.attachDofData(*Sn, "Sn", isBox);
writer.attachDofData(*Sw, "Sw", isBox);
writer.attachDofData(*pN, "pn", isBox);
diff --git a/dumux/implicit/2p2c/2p2cmodel.hh b/dumux/implicit/2p2c/2p2cmodel.hh
index 2af0d88b85..97ab736e73 100644
--- a/dumux/implicit/2p2c/2p2cmodel.hh
+++ b/dumux/implicit/2p2c/2p2cmodel.hh
@@ -375,9 +375,6 @@ public:
} // loop over elements
- velocityOutput.completeVelocityCalculation(*velocityW);
- velocityOutput.completeVelocityCalculation(*velocityN);
-
writer.attachDofData(*sN, "Sn", isBox);
writer.attachDofData(*sW, "Sw", isBox);
writer.attachDofData(*pN, "pN", isBox);
diff --git a/dumux/implicit/3p3c/3p3cmodel.hh b/dumux/implicit/3p3c/3p3cmodel.hh
index 2add58906c..b171e3283c 100644
--- a/dumux/implicit/3p3c/3p3cmodel.hh
+++ b/dumux/implicit/3p3c/3p3cmodel.hh
@@ -379,10 +379,6 @@ public:
}
- velocityOutput.completeVelocityCalculation(*velocityW);
- velocityOutput.completeVelocityCalculation(*velocityN);
- velocityOutput.completeVelocityCalculation(*velocityG);
-
writer.attachDofData(*saturation[wPhaseIdx], "Sw", isBox);
writer.attachDofData(*saturation[nPhaseIdx], "Sn", isBox);
writer.attachDofData(*saturation[gPhaseIdx], "Sg", isBox);
diff --git a/dumux/implicit/common/implicitvelocityoutput.hh b/dumux/implicit/common/implicitvelocityoutput.hh
index a8b752bdb8..455aadcebf 100644
--- a/dumux/implicit/common/implicitvelocityoutput.hh
+++ b/dumux/implicit/common/implicitvelocityoutput.hh
@@ -70,26 +70,41 @@ public:
* \param problem The problem to be solved
*/
ImplicitVelocityOutput(const Problem& problem)
- : problem_(problem)
+ : problem_(problem)
{
// check, if velocity output can be used (works only for cubes so far)
velocityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
+ if (velocityOutput_ && isBox)
+ {
+ cellNum_.assign(problem_.gridView().size(dofCodim), 0);
+ }
+
ElementIterator elemIt = problem_.gridView().template begin<0>();
ElementIterator elemEndIt = problem_.gridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt)
{
- if (elemIt->geometry().type().isCube() == false){
+ if (elemIt->geometry().type().isCube() == false)
+ {
velocityOutput_ = false;
}
+
+ if (isBox)
+ {
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(problem_.gridView(), *elemIt);
+
+ // transform vertex velocities from local to global coordinates
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ int globalIdx = problem_.vertexMapper().map(*elemIt, scvIdx, dofCodim);
+
+ cellNum_[globalIdx] += 1;
+ }
+ }
}
if (velocityOutput_ != GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity))
std::cout << "ATTENTION: Velocity output only works for cubes and is set to false for simplices\n";
-
- if (velocityOutput_ && isBox)
- {
- cellNum_.resize(problem_.gridView().size(dofCodim), 0);
- }
}
bool enableOutput()
@@ -137,7 +152,7 @@ public:
fvGeometry,
faceIdx,
elemVolVars);
-
+
const GlobalPosition globalNormal = fluxVars.face().normal;
GlobalPosition localNormal(0);
@@ -149,8 +164,7 @@ public:
// Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolume
// face in the reference element.
- Scalar flux;
- flux = fluxVars.volumeFlux(phaseIdx) / localArea;
+ Scalar flux = fluxVars.volumeFlux(phaseIdx) / localArea;
// transform the normal Darcy velocity into a vector
tmpVelocity = localNormal;
@@ -168,11 +182,9 @@ public:
Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
jacobianT2.mtv(scvVelocities[scvIdx], scvVelocity);
- scvVelocity /= element.geometry().integrationElement(localPos);
+ scvVelocity /= element.geometry().integrationElement(localPos)*cellNum_[globalIdx];
// add up the wetting phase subcontrolvolume velocities for each vertex
velocity[globalIdx] += scvVelocity;
-
- cellNum_[globalIdx] +=1;
}
}
else
@@ -206,12 +218,12 @@ public:
fvGeometry,
faceIdx,
elemVolVars,true);
-
+
Scalar flux = fluxVars.volumeFlux(phaseIdx);
scvVelocities[faceIdx] = flux;
}
}
-
+
Dune::FieldVector < CoordScalar, dim > refVelocity(0);
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (scvVelocities[2*i + 1] - scvVelocities[2*i]);
@@ -228,32 +240,10 @@ public:
} // velocity output
}
- template<class VelocityVector>
- void completeVelocityCalculation(VelocityVector& velocity)
- {
- if (velocityOutput_ && isBox)
- {
- unsigned int size = velocity.size();
-
- assert(size != cellNum_.size());
-
-
- // divide the vertex velocities by the number of adjacent scvs i.e. cells
- for(unsigned int globalIdx = 0; globalIdx < size; ++globalIdx)
- {
- velocity[globalIdx] /= cellNum_[globalIdx];
- }
- fluxVariables_.clear();
- }
- }
protected:
const Problem& problem_;
-
- // parameters given in constructor
bool velocityOutput_;
std::vector<int> cellNum_;
-
- std::vector<FluxVariables> fluxVariables_;
};
}
diff --git a/dumux/implicit/mpnc/mpncvtkwritercommon.hh b/dumux/implicit/mpnc/mpncvtkwritercommon.hh
index c386545c1c..df346904bc 100644
--- a/dumux/implicit/mpnc/mpncvtkwritercommon.hh
+++ b/dumux/implicit/mpnc/mpncvtkwritercommon.hh
@@ -65,7 +65,7 @@ class MPNCVtkWriterCommon : public MPNCVtkWriterModule<TypeTag>
typedef Dune::array<DimField, numPhases> PhaseDimField;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
-
+
public:
MPNCVtkWriterCommon(const Problem &problem)
: ParentType(problem), velocityOutput_(problem)
@@ -99,7 +99,7 @@ public:
if (velocityOutput_.enableOutput()) {
Scalar nDofs = this->problem_.model().numDofs();
- for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
velocity_[phaseIdx].resize(nDofs);
velocity_[phaseIdx] = 0;
}
@@ -127,7 +127,7 @@ public:
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
int globalIdx = this->problem_.model().dofMapper().map(element, scvIdx, dofCodim);
-
+
const VolumeVariables &volVars = elemVolVars[scvIdx];
if (porosityOutput_) porosity_[globalIdx] = volVars.porosity();
@@ -203,7 +203,6 @@ public:
if (velocityOutput_.enableOutput()) {
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- velocityOutput_.completeVelocityCalculation(velocity_[phaseIdx]);
std::ostringstream oss;
oss << "velocity_" << FluidSystem::phaseName(phaseIdx);
writer.attachDofData(velocity_[phaseIdx],
diff --git a/dumux/implicit/richards/richardsmodel.hh b/dumux/implicit/richards/richardsmodel.hh
index 59ac1d62ef..bc75f3279c 100644
--- a/dumux/implicit/richards/richardsmodel.hh
+++ b/dumux/implicit/richards/richardsmodel.hh
@@ -206,8 +206,6 @@ public:
velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *elemIt, /*phaseIdx=*/0);
}
- velocityOutput.completeVelocityCalculation(*velocity);
-
writer.attachDofData(*Sn, "Sn", isBox);
writer.attachDofData(*Sw, "Sw", isBox);
writer.attachDofData(*pN, "pn", isBox);
--
GitLab