From 475e9e31edd4b73d4c9fcadf321194f8099e7f39 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Tue, 21 May 2013 13:51:55 +0000
Subject: [PATCH] implicit 1p: add velocity calculation and output. Reviewed by
Christoph.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10721 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/1p/1pmodel.hh | 43 ++++++++++++++++++-------
dumux/implicit/1p/1pproperties.hh | 1 +
dumux/implicit/1p/1ppropertydefaults.hh | 3 ++
3 files changed, 35 insertions(+), 12 deletions(-)
diff --git a/dumux/implicit/1p/1pmodel.hh b/dumux/implicit/1p/1pmodel.hh
index edda3fa5c1..132ab64234 100644
--- a/dumux/implicit/1p/1pmodel.hh
+++ b/dumux/implicit/1p/1pmodel.hh
@@ -27,6 +27,7 @@
#ifndef DUMUX_1P_MODEL_HH
#define DUMUX_1P_MODEL_HH
+#include <dumux/implicit/common/implicitvelocityoutput.hh>
#include "1pproperties.hh"
namespace Dumux
@@ -58,6 +59,7 @@ class OnePBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -79,19 +81,27 @@ public:
MultiWriter &writer)
{
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+ typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
// create the required scalar fields
unsigned numDofs = this->numDofs();
ScalarField *p = writer.allocateManagedBuffer(numDofs);
ScalarField *K = writer.allocateManagedBuffer(numDofs);
+ VectorField *velocity = writer.template allocateManagedBuffer<double, dim>(numDofs);
+ ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
+
+ if (velocityOutput.enableOutput())
+ {
+ // initialize velocity field
+ for (unsigned int i = 0; i < numDofs; ++i)
+ {
+ (*velocity)[i] = double(0);
+ }
+ }
unsigned numElements = this->gridView_().size(0);
ScalarField *rank = writer.allocateManagedBuffer(numElements);
- FVElementGeometry fvGeometry;
- VolumeVariables volVars;
- ElementBoundaryTypes elemBcTypes;
-
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt)
@@ -99,30 +109,39 @@ public:
int idx = this->problem_().model().elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
+ FVElementGeometry fvGeometry;
fvGeometry.update(this->gridView_(), *elemIt);
- elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
+
+ ElementVolumeVariables elemVolVars;
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvGeometry,
+ false /* oldSol? */);
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
- volVars.update(sol[globalIdx],
- this->problem_(),
- *elemIt,
- fvGeometry,
- scvIdx,
- false);
const SpatialParams &spatialParams = this->problem_().spatialParams();
- (*p)[globalIdx] = volVars.pressure();
+ (*p)[globalIdx] = elemVolVars[scvIdx].pressure();
(*K)[globalIdx] = spatialParams.intrinsicPermeability(*elemIt,
fvGeometry,
scvIdx);
}
+
+ // velocity output
+ velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *elemIt, /*phaseIdx=*/0);
}
+ velocityOutput.completeVelocityCalculation(*velocity);
+
writer.attachDofData(*p, "p", isBox);
writer.attachDofData(*K, "K", isBox);
+ if (velocityOutput.enableOutput())
+ {
+ writer.attachDofData(*velocity, "velocity", isBox, dim);
+ }
writer.attachCellData(*rank, "process rank");
}
};
diff --git a/dumux/implicit/1p/1pproperties.hh b/dumux/implicit/1p/1pproperties.hh
index 9762091c06..4eeb4475e7 100644
--- a/dumux/implicit/1p/1pproperties.hh
+++ b/dumux/implicit/1p/1pproperties.hh
@@ -60,6 +60,7 @@ NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered i
NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes
NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
// \}
}
diff --git a/dumux/implicit/1p/1ppropertydefaults.hh b/dumux/implicit/1p/1ppropertydefaults.hh
index 1a7ae7aced..a7a4b2aedf 100644
--- a/dumux/implicit/1p/1ppropertydefaults.hh
+++ b/dumux/implicit/1p/1ppropertydefaults.hh
@@ -91,6 +91,9 @@ public:
typedef Dumux::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
};
+// disable velocity output by default
+SET_BOOL_PROP(OneP, VtkAddVelocity, false);
+
// enable gravity by default
SET_BOOL_PROP(OneP, ProblemEnableGravity, true);
--
GitLab