From 5620f46d428bd1f0eaf18077f50b3e9717d0acdf Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Tue, 21 May 2013 12:45:41 +0000
Subject: [PATCH] implicit 1p2c: replace the specific velocity calculation by
the general one. Reviewed by Christoph.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10717 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/1p2c/1p2cfluxvariables.hh | 47 +++--
dumux/implicit/1p2c/1p2cmodel.hh | 202 ++++----------------
dumux/implicit/1p2c/1p2cproperties.hh | 1 +
dumux/implicit/1p2c/1p2cpropertydefaults.hh | 3 +
4 files changed, 74 insertions(+), 179 deletions(-)
diff --git a/dumux/implicit/1p2c/1p2cfluxvariables.hh b/dumux/implicit/1p2c/1p2cfluxvariables.hh
index 2c4368bbb8..491a651b84 100644
--- a/dumux/implicit/1p2c/1p2cfluxvariables.hh
+++ b/dumux/implicit/1p2c/1p2cfluxvariables.hh
@@ -95,6 +95,8 @@ public:
const bool onBoundary = false)
: fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
{
+ mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
+
viscosity_ = Scalar(0);
molarDensity_ = Scalar(0);
density_ = Scalar(0);
@@ -233,15 +235,31 @@ public:
* \brief Return the local index of the upstream control volume
* for a given phase.
*/
- int upstreamIdx() const
- { return upstreamIdx_; }
+ int upstreamIdx() const
+ { return upstreamIdx_; }
- /*!
- * \brief Return the local index of the downstream control volume
- * for a given phase.
- */
- int downstreamIdx() const
- { return downstreamIdx_; }
+ /*!
+ * \brief Return the local index of the downstream control volume
+ * for a given phase.
+ */
+ int downstreamIdx() const
+ { return downstreamIdx_; }
+
+ /*!
+ * \brief Return the volumetric flux over a face of a given phase.
+ *
+ * This is the calculated velocity multiplied by the unit normal
+ * and the area of the face.
+ * face().normal
+ * has already the magnitude of the area.
+ *
+ * \param phaseIdx index of the phase
+ */
+ Scalar volumeFlux(const unsigned int phaseIdx) const
+ {
+ assert (phaseIdx == Indices::phaseIdx);
+ return volumeFlux_;
+ }
protected:
@@ -401,6 +419,10 @@ protected:
std::swap(upstreamIdx_,
downstreamIdx_);
}
+
+ volumeFlux_ = KmvpNormal_;
+ volumeFlux_ *= mobilityUpwindWeight_/elemVolVars[upstreamIdx_].viscosity()
+ + (1.0 - mobilityUpwindWeight_)/elemVolVars[downstreamIdx_].viscosity();
}
/*!
* \brief Calculation of the effective diffusion coefficient.
@@ -493,15 +515,18 @@ protected:
Scalar KmvpNormal_;
// local index of the upwind vertex for each phase
- int upstreamIdx_;
- // local index of the downwind vertex for each phase
- int downstreamIdx_;
+ int upstreamIdx_;
+ // local index of the downwind vertex for each phase
+ int downstreamIdx_;
//! viscosity of the fluid at the integration point
Scalar viscosity_;
//! molar densities of the fluid at the integration point
Scalar molarDensity_, density_;
+
+ Scalar volumeFlux_; //!< Velocity multiplied with normal (magnitude=area)
+ Scalar mobilityUpwindWeight_; //!< Upwind weight for mobility. Set to one for full upstream weighting
};
} // end namepace
diff --git a/dumux/implicit/1p2c/1p2cmodel.hh b/dumux/implicit/1p2c/1p2cmodel.hh
index 174400ec0c..110070e540 100644
--- a/dumux/implicit/1p2c/1p2cmodel.hh
+++ b/dumux/implicit/1p2c/1p2cmodel.hh
@@ -27,6 +27,7 @@
#ifndef DUMUX_ONEP_TWOC_MODEL_HH
#define DUMUX_ONEP_TWOC_MODEL_HH
+#include <dumux/implicit/common/implicitvelocityoutput.hh>
#include "1p2cproperties.hh"
namespace Dumux
@@ -86,6 +87,9 @@ class OnePTwoCBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum { phaseIdx = Indices::phaseIdx };
+
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dim> DimVector;
@@ -94,15 +98,10 @@ class OnePTwoCBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
public:
/*!
- * \brief Constructor. Sets the upwind weight.
+ * \brief Constructor
*/
OnePTwoCBoxModel()
- {
- // retrieve the upwind weight for the mass conservation equations. Use the value
- // specified via the property system as default, and overwrite
- // it by the run-time parameter from the Dune::ParameterTree
- upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
- }
+ {}
/*!
* \brief \copybrief ImplicitModel::addOutputVtkFields
@@ -114,8 +113,8 @@ public:
void addOutputVtkFields(const SolutionVector &sol,
MultiWriter &writer)
{
- bool velocityOutput = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
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();
@@ -127,41 +126,20 @@ public:
ScalarField &massFraction1 = *writer.allocateManagedBuffer(numDofs);
ScalarField &rho = *writer.allocateManagedBuffer(numDofs);
ScalarField &mu = *writer.allocateManagedBuffer(numDofs);
- ScalarField &velocityX = *writer.allocateManagedBuffer(numDofs);
- ScalarField &velocityY = *writer.allocateManagedBuffer(numDofs);
- ScalarField &velocityZ = *writer.allocateManagedBuffer(numDofs);
- //use vertical faces for vx and horizontal faces for vy calculation
- std::vector<DimVector> boxSurface(numDofs);
-
- // velocity output currently only works for the box discretization
- if (!isBox)
- velocityOutput = false;
-
- if (velocityOutput)
+ VectorField *velocity = writer.template allocateManagedBuffer<double, dim>(numDofs);
+ ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
+
+ if (velocityOutput.enableOutput())
{
- // initialize velocity fields
+ // initialize velocity field
for (unsigned int i = 0; i < numDofs; ++i)
{
- velocityX[i] = 0;
- if (dim > 1)
- {
- velocityY[i] = 0;
- }
- if (dim > 2)
- {
- velocityZ[i] = 0;
- }
- boxSurface[i] = Scalar(0.0); // initialize the boundary surface of the fv-boxes
+ (*velocity)[i] = Scalar(0);
}
}
unsigned numElements = this->gridView_().size(0);
- ScalarField &rank =
- *writer.allocateManagedBuffer(numElements);
-
- FVElementGeometry fvGeometry;
- VolumeVariables volVars;
- ElementBoundaryTypes elemBcTypes;
+ ScalarField &rank = *writer.allocateManagedBuffer(numElements);
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
@@ -170,149 +148,40 @@ 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);
-
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
- volVars.update(sol[globalIdx],
- this->problem_(),
+ ElementVolumeVariables elemVolVars;
+ elemVolVars.update(this->problem_(),
*elemIt,
fvGeometry,
- scvIdx,
- false);
-
- pressure[globalIdx] = volVars.pressure();
- delp[globalIdx] = volVars.pressure() - 1e5;
- moleFraction0[globalIdx] = volVars.moleFraction(0);
- moleFraction1[globalIdx] = volVars.moleFraction(1);
- massFraction0[globalIdx] = volVars.massFraction(0);
- massFraction1[globalIdx] = volVars.massFraction(1);
- rho[globalIdx] = volVars.density();
- mu[globalIdx] = volVars.viscosity();
- }
+ false /* oldSol? */);
- if (velocityOutput)
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
- // 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.
- DimVector velocity;
-
- ElementVolumeVariables elemVolVars;
- elemVolVars.update(this->problem_(),
- *elemIt,
- fvGeometry,
- false /* isOldSol? */);
- // loop over the phases
- for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; faceIdx++)
- {
- velocity = 0.0;
- //prepare the flux calculations (set up and prepare geometry, FE gradients)
- FluxVariables fluxVars(this->problem_(),
- *elemIt,
- fvGeometry,
- faceIdx,
- elemVolVars);
-
- //use vertical faces for vx and horizontal faces for vy calculation
- DimVector xVector(0), yVector(0);
- xVector[0] = 1; yVector[1] = 1;
-
- Dune::SeqScalarProduct<DimVector> sp;
-
- Scalar xDir = std::abs(sp.dot(fluxVars.face().normal, xVector));
- Scalar yDir = std::abs(sp.dot(fluxVars.face().normal, yVector));
-
- // up+downstream node
- const VolumeVariables &up =
- elemVolVars[fluxVars.upstreamIdx()];
- const VolumeVariables &dn =
- elemVolVars[fluxVars.downstreamIdx()];
-
- //get surface area to weight velocity at the IP with the surface area
- Scalar scvfArea = fluxVars.face().normal.two_norm();
-
- int vertIIdx = this->problem_().vertexMapper().map(
- *elemIt, fluxVars.face().i, dim);
- int vertJIdx = this->problem_().vertexMapper().map(
- *elemIt, fluxVars.face().j, dim);
-
- //use vertical faces (horizontal noraml vector) to calculate vx
- //in case of heterogeneities it seams to be better to define intrinisc permeability elementwise
- if (xDir > yDir)//(fluxVars.face().normal[0] > 1e-10 || fluxVars.face().normal[0] < -1e-10)// (xDir > yDir)
- {
- // get darcy velocity
- //calculate (v n) n/A
- Scalar tmp = fluxVars.KmvpNormal();
- velocity = fluxVars.face().normal;
- velocity *= tmp;
- velocity /= scvfArea;
- velocity *= (upwindWeight_ / up.viscosity() +
- (1 - upwindWeight_)/ dn.viscosity());
-
- // add surface area for weighting purposes
- boxSurface[vertIIdx][0] += scvfArea;
- boxSurface[vertJIdx][0] += scvfArea;
-
- velocityX[vertJIdx] += velocity[0];
- velocityX[vertIIdx] += velocity[0];
-
- }
- if (yDir > xDir)//(fluxVars.face().normal[1] > 1e-10 || fluxVars.face().normal[1] < -1e-10)// (yDir > xDir)
- {
- // get darcy velocity
- //calculate (v n) n/A
- Scalar tmp = fluxVars.KmvpNormal();
- velocity = fluxVars.face().normal;
- velocity *= tmp;
- velocity /= scvfArea;
- velocity *= (upwindWeight_ / up.viscosity() +
- (1 - upwindWeight_)/ dn.viscosity());
-
- // add surface area for weighting purposes
- boxSurface[vertIIdx][1] += scvfArea;
- boxSurface[vertJIdx][1] += scvfArea;
+ int globalIdx = this->dofMapper().map(*elemIt, scvIdx, dofCodim);
- velocityY[vertJIdx] += velocity[1];
- velocityY[vertIIdx] += velocity[1];
- }
- }
+ pressure[globalIdx] = elemVolVars[scvIdx].pressure();
+ delp[globalIdx] = elemVolVars[scvIdx].pressure() - 1e5;
+ moleFraction0[globalIdx] = elemVolVars[scvIdx].moleFraction(0);
+ moleFraction1[globalIdx] = elemVolVars[scvIdx].moleFraction(1);
+ massFraction0[globalIdx] = elemVolVars[scvIdx].massFraction(0);
+ massFraction1[globalIdx] = elemVolVars[scvIdx].massFraction(1);
+ rho[globalIdx] = elemVolVars[scvIdx].density();
+ mu[globalIdx] = elemVolVars[scvIdx].viscosity();
}
- }
-
- if (velocityOutput)
- {
- // normalize the velocities at the vertices
- // calculate the bounding box of the grid view
- VertexIterator vIt = this->gridView_().template begin<dim>();
- const VertexIterator vEndIt = this->gridView_().template end<dim>();
- for (; vIt!=vEndIt; ++vIt)
- {
- int i = this->problem_().vertexMapper().map(*vIt);
- //use vertical faces for vx and horizontal faces for vy calculation
- velocityX[i] /= boxSurface[i][0];
- if (dim >= 2)
- {
- velocityY[i] /= boxSurface[i][1];
- }
- if (dim == 3)
- {
- velocityZ[i] /= boxSurface[i][2];
- }
- }
+ // velocity output
+ velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *elemIt, phaseIdx);
}
+ velocityOutput.completeVelocityCalculation(*velocity);
+
writer.attachDofData(pressure, "P", isBox);
writer.attachDofData(delp, "delp", isBox);
- if (velocityOutput)
+ if (velocityOutput.enableOutput())
{
- writer.attachDofData(velocityX, "Vx", isBox);
- writer.attachDofData(velocityY, "Vy", isBox);
- if (dim > 2)
- writer.attachDofData(velocityZ, "Vz", isBox);
+ writer.attachDofData(*velocity, "velocity", isBox, dim);
}
char nameMoleFraction0[42], nameMoleFraction1[42];
snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0));
@@ -329,9 +198,6 @@ public:
writer.attachDofData(mu, "mu", isBox);
writer.attachCellData(rank, "process rank");
}
-
-private:
- Scalar upwindWeight_;
};
}
diff --git a/dumux/implicit/1p2c/1p2cproperties.hh b/dumux/implicit/1p2c/1p2cproperties.hh
index f81859414e..5f28543b1a 100644
--- a/dumux/implicit/1p2c/1p2cproperties.hh
+++ b/dumux/implicit/1p2c/1p2cproperties.hh
@@ -59,6 +59,7 @@ NEW_PROP_TAG(Indices); //!< Enumerations for the model
NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The default value of the upwind weight
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(UseMoles); //!Defines whether mole (true) or mass (false) fractions are used
NEW_PROP_TAG(Scaling); //!Defines Scaling of the model
diff --git a/dumux/implicit/1p2c/1p2cpropertydefaults.hh b/dumux/implicit/1p2c/1p2cpropertydefaults.hh
index b61b67ca92..a84fd9c12d 100644
--- a/dumux/implicit/1p2c/1p2cpropertydefaults.hh
+++ b/dumux/implicit/1p2c/1p2cpropertydefaults.hh
@@ -69,6 +69,9 @@ SET_TYPE_PROP(OnePTwoC, FluxVariables, OnePTwoCFluxVariables<TypeTag>);
//! set default upwind weight to 1.0, i.e. fully upwind
SET_SCALAR_PROP(OnePTwoC, ImplicitMassUpwindWeight, 1.0);
+//! weight for the upwind mobility in the velocity calculation
+SET_SCALAR_PROP(OnePTwoC, ImplicitMobilityUpwindWeight, 1.0);
+
//! Set the indices used by the 1p2c model
SET_TYPE_PROP(OnePTwoC, Indices, Dumux::OnePTwoCIndices<TypeTag, 0>);
--
GitLab