Commit 8b4dedd7 authored by Katherina Baber's avatar Katherina Baber
Browse files

included velocity output

added calculateVelocity() to fluxvariables (Kmvp, KmvpNormal) and
densityAtIp
localresidual uses now massfractions and the velocities calculated in
the fluxvariables


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@5345 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent f1188b2c
......@@ -183,10 +183,6 @@ protected:
potentialGrad_ -= tmp;
// Scalar k = problem.spatialParameters().intrinsicPermeability(element, fvElemGeom_, scvIdx_);
// VectorGradient K(0);
// K[0][0] = K[1][1] = k;
//TODO use meanK
VectorGradient K = problem.spatialParameters().intrinsicPermeability(element, fvElemGeom_, scvIdx_);
ScalarGradient Kmvp;
K.mv(potentialGrad_, Kmvp);
......
......@@ -97,14 +97,41 @@ public:
viscosityAtIP_ = Scalar(0);
molarDensityAtIP_ = Scalar(0);
densityAtIP_ = Scalar(0);
potentialGrad_ = Scalar(0);
concentrationGrad_ = Scalar(0);
molarConcGrad_ = Scalar(0);
calculateGradients_(problem, element, elemDat);
calculateK_(problem, element, elemDat);
calculateVelocities_(problem, element, elemDat);
calculateDiffCoeffPM_(problem, element, elemDat);
calculateDispersionTensor_(problem, element, elemDat);
};
public:
/*!
* \brief Return the pressure potential multiplied with the
* intrinsic permeability which goes from vertex i to
* vertex j.
*
* Note that the length of the face's normal is the area of the
* phase, so this is not the actual velocity by the integral of
* the velocity over the face's area. Also note that the phase
* mobility is not yet included here since this would require a
* decision on the upwinding approach (which is done in the
* actual model).
*/
Scalar KmvpNormal() const
{ return KmvpNormal_; }
/*!
* \brief Return the pressure potential multiplied with the
* intrinsic permeability as vector (for velocity output)
*/
Vector Kmvp() const
{ return Kmvp_; }
const SCVFace &face() const
{ return fvElemGeom_.subContVolFace[scvfIdx_]; }
......@@ -166,6 +193,8 @@ public:
Scalar molarDensityAtIP() const
{ return molarDensityAtIP_; }
Scalar densityAtIP() const
{ return densityAtIP_; }
/*!
* \brief Given the intrinisc permeability times the pressure
......@@ -189,6 +218,20 @@ public:
int downstreamIdx(Scalar normalFlux) const
{ return (normalFlux > 0)?face().j:face().i; }
/*!
* \brief Return the local index of the upstream control volume
* for a given phase.
*/
int upstreamIdx() const
{ return upstreamIdx_; }
/*!
* \brief Return the local index of the downstream control volume
* for a given phase.
*/
int downstreamIdx() const
{ return downstreamIdx_; }
protected:
/*!
......@@ -205,10 +248,6 @@ protected:
const VolumeVariables &vVars_i = elemDat[face().i];
const VolumeVariables &vVars_j = elemDat[face().j];
potentialGrad_ = 0.0;
concentrationGrad_ = 0.0;
molarConcGrad_ = 0.0;
Vector tmp;
//The decision of the if-statement depends on the function useTwoPointGradient(const Element &elem,
//int vertexI,int vertexJ) defined in test/tissue_tumor_spatialparameters.hh
......@@ -240,6 +279,9 @@ protected:
//phase moledensity
molarDensityAtIP_ += elemDat[idx].molarDensity()*face().shapeValue[idx];
//phase density
densityAtIP_ += elemDat[idx].density()*face().shapeValue[idx];
}
}
else {
......@@ -289,8 +331,26 @@ protected:
sp.intrinsicPermeability(element,
fvElemGeom_,
face().j));
}
void calculateVelocities_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
{
K_.mv(potentialGrad_, Kmvp_);
KmvpNormal_ = - (Kmvp_ * face().normal);
// set the upstream and downstream vertices
upstreamIdx_ = face().i;
downstreamIdx_ = face().j;
if (KmvpNormal_ < 0)
{
std::swap(upstreamIdx_,
downstreamIdx_);
}
}
/*!
* \brief Calculation of the effective diffusion coefficient
*
......@@ -307,8 +367,10 @@ protected:
// Diffusion coefficient in the porous medium
diffCoeffPM_
= 1./2*(vDat_i.porosity() * vDat_i.tortuosity() * vDat_i.diffCoeff() +
vDat_j.porosity() * vDat_j.tortuosity() * vDat_j.diffCoeff());
= harmonicMean(vDat_i.porosity() * vDat_i.tortuosity() * vDat_i.diffCoeff(),
vDat_j.porosity() * vDat_j.tortuosity() * vDat_j.diffCoeff());
// = 1./2*(vDat_i.porosity() * vDat_i.tortuosity() * vDat_i.diffCoeff() +
// vDat_j.porosity() * vDat_j.tortuosity() * vDat_j.diffCoeff());
}
/*!
......@@ -376,12 +438,21 @@ protected:
//! the intrinsic permeability tensor
Tensor K_;
// intrinsic permeability times pressure potential gradient
Vector Kmvp_;
// projected on the face normal
Scalar KmvpNormal_;
// local index of the upwind vertex for each phase
int upstreamIdx_;
// local index of the downwind vertex for each phase
int downstreamIdx_;
//! viscosity of the fluid at the integration point
Scalar viscosityAtIP_;
//! molar densities of the fluid at the integration point
Scalar molarDensityAtIP_;
Scalar molarDensityAtIP_, densityAtIP_;
};
} // end namepace
......
......@@ -29,6 +29,7 @@
#ifndef DUMUX_ONEP_TWOC_LOCAL_RESIDUAL_HH
#define DUMUX_ONEP_TWOC_LOCAL_RESIDUAL_HH
#define VELOCITY_OUTPUT 1 //1 turns velocity output on, 0 turns it off
#include <dumux/boxmodels/common/boxmodel.hh>
......@@ -62,15 +63,22 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
enum
{
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
// indices of the primary variables
pressureIdx = Indices::pressureIdx,
......@@ -118,10 +126,15 @@ public:
result[contiEqIdx] =
volVars.density()*volVars.porosity();
// storage term of the transport equation
// storage term of the transport equation - molefractions
// result[transEqIdx] =
// volVars.concentration(comp1Idx) *
// volVars.porosity();
//storage term of the transport equation - massfractions
result[transEqIdx] =
volVars.concentration(comp1Idx) *
volVars.porosity();
volVars.density() * volVars.massFrac(comp1Idx) * volVars.porosity();
}
/*!
......@@ -156,28 +169,35 @@ public:
////////
// advective fluxes of all components in all phases
////////
Vector tmpVec(0);
fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(), tmpVec);
// "intrinsic" flux from cell i to cell j
Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up =
this->curVolVars_(fluxVars.upstreamIdx());
const VolumeVariables &dn =
this->curVolVars_(fluxVars.downstreamIdx());
// total mass flux
//KmvpNormal is the Darcy velocity multiplied with the normal vector, calculated in 1p2cfluxvariables.hh
flux[contiEqIdx] =
normalFlux *
fluxVars.KmvpNormal() *
(( upwindAlpha)*up.density()/up.viscosity()
+
((1 - upwindAlpha)*dn.density()/dn.viscosity()));
// advective flux of the second component
// advective flux of the second component -molefraction
// flux[transEqIdx] +=
// fluxVars.KmvpNormal() *
// (( upwindAlpha)*up.concentration(comp1Idx)/up.viscosity()
// +
// (1 - upwindAlpha)*dn.concentration(comp1Idx)/dn.viscosity());
// advective flux of the second component - massfraction
flux[transEqIdx] +=
normalFlux *
(( upwindAlpha)*up.concentration(comp1Idx)/up.viscosity()
fluxVars.KmvpNormal() *
(( upwindAlpha)*up.density() * up.massFrac(comp1Idx)/up.viscosity()
+
(1 - upwindAlpha)*dn.concentration(comp1Idx)/dn.viscosity());
(1 - upwindAlpha)*dn.density()*dn.massFrac(comp1Idx)/dn.viscosity());
}
......@@ -190,16 +210,22 @@ public:
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
Scalar tmp(0);
// diffusive flux of second component
flux[transEqIdx] -=
tmp -=
fluxVars.porousDiffCoeff() *
(fluxVars.concentrationGrad(1) * fluxVars.face().normal);
(fluxVars.concentrationGrad(comp1Idx) * fluxVars.face().normal);
// dispersive flux of second component
Vector normalDisp;
fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp);
flux[transEqIdx] -=
tmp -=
(normalDisp * fluxVars.concentrationGrad(comp1Idx));
//molar
// flux[transEqIdx] += tmp;
//transform to mass fractions
flux[transEqIdx] += tmp * FluidSystem::molarMass(comp1Idx);
}
/*!
* \brief Calculate the source term of the equation
......@@ -214,7 +240,7 @@ public:
this->fvElemGeom_(),
localVertexIdx);
}
Implementation *asImp_()
{ return static_cast<Implementation *> (this); }
const Implementation *asImp_() const
......
......@@ -34,6 +34,7 @@
#include "1p2cproperties.hh"
#include "1p2cproblem.hh"
#include "1p2clocalresidual.hh"
#include <dumux/boxmodels/common/boxmodel.hh>
......@@ -92,6 +93,7 @@ class OnePTwoCBoxModel : public BoxModel<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
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(VertexMapper)) VertexMapper;
......@@ -100,7 +102,12 @@ class OnePTwoCBoxModel : public BoxModel<TypeTag>
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
static const Scalar upwindAlpha = GET_PROP_VALUE(TypeTag, PTAG(UpwindAlpha));
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
public:
/*!
......@@ -120,11 +127,36 @@ public:
// create the required scalar fields
unsigned numVertices = this->problem_().gridView().size(dim);
ScalarField *pressure = writer.template createField<Scalar, 1>(numVertices);
ScalarField *delp = writer.template createField<Scalar, 1>(numVertices);
ScalarField *moleFrac0 = writer.template createField<Scalar, 1>(numVertices);
ScalarField *moleFrac1 = writer.template createField<Scalar, 1>(numVertices);
ScalarField *massFrac0 = writer.template createField<Scalar, 1>(numVertices);
ScalarField *massFrac1 = writer.template createField<Scalar, 1>(numVertices);
ScalarField *rho = writer.template createField<Scalar, 1>(numVertices);
ScalarField *mu = writer.template createField<Scalar, 1>(numVertices);
#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
ScalarField *velocityX = writer.template createField<Scalar, 1>(numVertices);
ScalarField *velocityY = writer.template createField<Scalar, 1>(numVertices);
ScalarField *velocityZ = writer.template createField<Scalar, 1>(numVertices);
//use vertiacl faces for vx and horizontal faces for vy calculation
GlobalPosition boxSurface[numVertices];
// initialize velocity fields
for (int i = 0; i < numVertices; ++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
}
#endif
unsigned numElements = this->gridView_().size(0);
ScalarField *rank =
writer.template createField<Scalar, 1> (numElements);
......@@ -154,20 +186,138 @@ public:
i,
false);
(*pressure)[globalIdx] = volVars.pressure();
(*delp)[globalIdx] = volVars.pressure() - 1e5;
(*moleFrac0)[globalIdx] = volVars.moleFrac(0);
(*moleFrac1)[globalIdx] = volVars.moleFrac(1);
(*massFrac0)[globalIdx] = volVars.massFrac(0);
(*massFrac1)[globalIdx] = volVars.massFrac(1);
(*rho)[globalIdx] = volVars.density();
(*mu)[globalIdx] = volVars.viscosity();
};
#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
// 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;
ElementVolumeVariables elemVolVars;
elemVolVars.update(this->problem_(),
*elemIt,
fvElemGeom,
false /* isOldSol? */);
// loop over the phases
for (int faceIdx = 0; faceIdx < fvElemGeom.numEdges; faceIdx++)
{
velocity = 0.0;
//prepare the flux calculations (set up and prepare geometry, FE gradients)
FluxVariables fluxVars(this->problem_(),
*elemIt,
fvElemGeom,
faceIdx,
elemVolVars);
//use vertiacl faces for vx and horizontal faces for vy calculation
GlobalPosition xVector(0), yVector(0);
xVector[0] = 1; yVector[1] = 1;
Dune::SeqScalarProduct<GlobalPosition> 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 *= (upwindAlpha / up.viscosity() +
(1 - upwindAlpha)/ 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 *= (upwindAlpha / up.viscosity() +
(1 - upwindAlpha)/ dn.viscosity());
// add surface area for weighting purposes
boxSurface[vertIIdx][1] += scvfArea;
boxSurface[vertJIdx][1] += scvfArea;
(*velocityY)[vertJIdx] += velocity[1];
(*velocityY)[vertIIdx] += velocity[1];
}
}
#endif
}
#ifdef VELOCITY_OUTPUT
// normalize the velocities at the vertices
// calculate the bounding box of the grid view
int index = 0;
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 vertiacl 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];
}
}
#endif
writer.addVertexData(pressure, "p");
writer.addVertexData(moleFrac0, "x_0");
writer.addVertexData(moleFrac1, "x_1");
writer.addVertexData(massFrac0, "X_0");
writer.addVertexData(massFrac1, "X_1");
writer.addVertexData(delp, "delp");
#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
writer.addVertexData(velocityX, "Vx");
writer.addVertexData(velocityY, "Vy");
if (dim > 2)
writer.addVertexData(velocityZ, "Vz");
#endif
writer.addVertexData(moleFrac0, "x_if");
writer.addVertexData(moleFrac1, "x_TRAIL");
writer.addVertexData(massFrac0, "X_if");
writer.addVertexData(massFrac1, "X_TRAIL");
writer.addVertexData(rho, "rho");
writer.addVertexData(mu, "mu");
writer.addCellData(rank, "process rank");
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment