Commit 79ed65c9 authored by Klaus Mosthaf's avatar Klaus Mosthaf
Browse files

Further simplified the evaluation of the outflow boundary conditions in the...

Further simplified the evaluation of the outflow boundary conditions in the 2p2c and 2p2cni model. Added bool onBoundary as argument to computeFlux(). Then, the fluxVariables object that is generated in a correct way, depending whether we are on a boundary or not and computeOutflowValues() becomes superfluous (it did the same as computeFlux(), just with a BoundaryVariables object and evaluation at the boundary instead of the integration point of the sub-control-volume face).
Removed the <FaceType> template parameter in the fluxvariables, it is not required anymore. The face is also not needed any further as argument, it should to be called via the new face() method, which returns either a boundaryFace or a inner SCV face.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7403 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent cecb7bb4
......@@ -104,25 +104,21 @@ public:
molarConcGrad_[phaseIdx] = Scalar(0);
}
calculateValues_(problem, element, face(), elemVolVars);
calculateValues_(problem, element, elemVolVars);
}
protected:
template <class FaceType>
void calculateValues_(const Problem &problem,
const Element &element,
const FaceType &face,
const ElementVolumeVariables &elemVolVars)
{
calculateGradients_(problem, element, face, elemVolVars);
calculateVelocities_(problem, element, face, elemVolVars);
calculateDiffCoeffPM_(problem, element, face, elemVolVars);
calculateGradients_(problem, element, elemVolVars);
calculateVelocities_(problem, element, elemVolVars);
calculateDiffCoeffPM_(problem, element, elemVolVars);
}
template <class FaceType>
void calculateGradients_(const Problem &problem,
const Element &element,
const FaceType &face,
const ElementVolumeVariables &elemVolVars)
{
// calculate gradients
......@@ -132,7 +128,7 @@ protected:
idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
const Vector &feGrad = face.grad[idx];
const Vector &feGrad = face().grad[idx];
// compute sum of pressure gradients for each phase
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
......@@ -170,18 +166,18 @@ protected:
if (GET_PARAM(TypeTag, bool, EnableGravity)) {
// estimate the gravitational acceleration at a given SCV face
// using the arithmetic mean
Vector g(problem.boxGravity(element, fvGeom_, face.i));
g += problem.boxGravity(element, fvGeom_, face.j);
Vector g(problem.boxGravity(element, fvGeom_, face().i));
g += problem.boxGravity(element, fvGeom_, face().j);
g /= 2;
for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
{
// calculate the phase density at the integration point. we
// only do this if the wetting phase is present in both cells
Scalar SI = elemVolVars[face.i].saturation(phaseIdx);
Scalar SJ = elemVolVars[face.j].saturation(phaseIdx);
Scalar rhoI = elemVolVars[face.i].density(phaseIdx);
Scalar rhoJ = elemVolVars[face.j].density(phaseIdx);
Scalar SI = elemVolVars[face().i].saturation(phaseIdx);
Scalar SJ = elemVolVars[face().j].saturation(phaseIdx);
Scalar rhoI = elemVolVars[face().i].density(phaseIdx);
Scalar rhoJ = elemVolVars[face().j].density(phaseIdx);
Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
if (fI + fJ == 0)
......@@ -215,10 +211,8 @@ protected:
return sp.eval(sat);
}
template <class FaceType>
void calculateVelocities_(const Problem &problem,
const Element &element,
const FaceType &face,
const ElementVolumeVariables &elemVolVars)
{
const SpatialParameters &spatialParams = problem.spatialParameters();
......@@ -229,19 +223,19 @@ protected:
spatialParams.meanK(K,
spatialParams.intrinsicPermeability(element,
fvGeom_,
face.i),
face().i),
spatialParams.intrinsicPermeability(element,
fvGeom_,
face.j));
face().j));
K.mv(potentialGrad_[phaseIdx], Kmvp_[phaseIdx]);
KmvpNormal_[phaseIdx] = -(Kmvp_[phaseIdx]*face.normal);
KmvpNormal_[phaseIdx] = -(Kmvp_[phaseIdx]*face().normal);
}
// set the upstream and downstream vertices
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
upstreamIdx_[phaseIdx] = face.i;
downstreamIdx_[phaseIdx] = face.j;
upstreamIdx_[phaseIdx] = face().i;
downstreamIdx_[phaseIdx] = face().j;
if (KmvpNormal_[phaseIdx] < 0) {
std::swap(upstreamIdx_[phaseIdx],
......@@ -250,14 +244,12 @@ protected:
}
}
template <class FaceType>
void calculateDiffCoeffPM_(const Problem &problem,
const Element &element,
const FaceType &face,
const ElementVolumeVariables &elemVolVars)
{
const VolumeVariables &vDat_i = elemVolVars[face.i];
const VolumeVariables &vDat_j = elemVolVars[face.j];
const VolumeVariables &vDat_i = elemVolVars[face().i];
const VolumeVariables &vDat_j = elemVolVars[face().j];
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
......@@ -355,12 +347,16 @@ public:
const Vector &molarConcGrad(int phaseIdx) const
{ return molarConcGrad_[phaseIdx]; };
/*!
* \brief The face of the current sub-control volume. This may be either
* an inner sub-control-volume face or a face on the boundary.
*/
const SCVFace &face() const
{
if (onBoundary_)
return fvGeom_.boundaryFace[faceIdx_];
else
return fvGeom_.subContVolFace[faceIdx_];
if (onBoundary_)
return fvGeom_.boundaryFace[faceIdx_];
else
return fvGeom_.subContVolFace[faceIdx_];
}
protected:
......
......@@ -146,21 +146,13 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
int scvIdx,
int boundaryFaceIdx)
{
// temporary vector to store the outflow boundary fluxes
PrimaryVariables flux(0.0);
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
// deal with outflow boundaries
if (bcTypes.hasOutflow())
{
const BoundaryVariables boundaryVars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
boundaryFaceIdx,
this->curVolVars_());
PrimaryVariables values(0.0);
asImp_()->computeOutflowValues(values, boundaryVars, scvIdx, boundaryFaceIdx);
asImp_()->computeFlux(values, boundaryFaceIdx, true);
Valgrind::CheckDefined(values);
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
......@@ -223,13 +215,14 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
* \param flux The flux over the SCV (sub-control-volume) face for each component
* \param faceIdx The index of the SCV face
*/
void computeFlux(PrimaryVariables &flux, int faceIdx) const
void computeFlux(PrimaryVariables &flux, int faceIdx, bool onBoundary=false) const
{
FluxVariables vars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
faceIdx,
this->curVolVars_());
this->curVolVars_(),
onBoundary);
flux = 0;
asImp_()->computeAdvectiveFlux(flux, vars);
......@@ -368,78 +361,6 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
this->curVolVars_());
}
/*!
* \brief Compute the fluxes at outflow boundaries. This does essentially the same
* as computeFluxes, but the fluxes are evaluated at the integration point
* of the boundary face. However, some variables are evaluated
* at the vertex (usually the ones which are upwinded).
*
* \param flux A temporary vector, where the outflow boundary fluxes are stored
* \param boundaryVars The boundary variables object
* \param scvIdx The index of the SCV containing the outflow boundary face
* \param boundaryFaceIdx The index of the boundary face
*/
void computeOutflowValues(PrimaryVariables &flux,
const BoundaryVariables &boundaryVars,
const int scvIdx,
const int boundaryFaceIdx)
{
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
const VolumeVariables& vertVars = this->curVolVars_()[scvIdx];
// component transport
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (int compIdx = contiCompIdx1_(); compIdx <= contiCompIdx2_(); ++compIdx)
{
int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
if (bcTypes.isOutflow(eqIdx))
{
// advective flux
flux[eqIdx] += boundaryVars.KmvpNormal(phaseIdx)
* vertVars.density(phaseIdx)
* vertVars.mobility(phaseIdx)
* vertVars.fluidState().massFraction(phaseIdx, compIdx);
}
}
}
// add diffusive flux of gas component in liquid phase
Scalar tmp = boundaryVars.molarConcGrad(lPhaseIdx)*boundaryVars.face().normal;
tmp *= -1;
tmp *= boundaryVars.porousDiffCoeff(lPhaseIdx) *
boundaryVars.molarDensityAtIP(lPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiGEqIdx && bcTypes.isOutflow(contiGEqIdx))
flux[contiGEqIdx] += tmp * FluidSystem::molarMass(gCompIdx);
if (replaceCompEqIdx != contiLEqIdx && bcTypes.isOutflow(contiLEqIdx))
flux[contiLEqIdx] -= tmp * FluidSystem::molarMass(lCompIdx);
// add diffusive flux of liquid component in gas phase
tmp = boundaryVars.molarConcGrad(gPhaseIdx)*boundaryVars.face().normal;
tmp *= -1;
tmp *= boundaryVars.porousDiffCoeff(gPhaseIdx) *
boundaryVars.molarDensityAtIP(gPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiLEqIdx && bcTypes.isOutflow(contiLEqIdx))
flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
if (replaceCompEqIdx != contiGEqIdx && bcTypes.isOutflow(contiGEqIdx))
flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
// flux of the total mass balance;
// this is only processed, if one component mass-balance equation
// is replaced by a total mass-balance equation
if (replaceCompEqIdx < numComponents)
{
if (bcTypes.isOutflow(replaceCompEqIdx))
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
flux[replaceCompEqIdx] += boundaryVars.KmvpNormal(phaseIdx)
*vertVars.density(phaseIdx)
*vertVars.mobility(phaseIdx);
}
}
protected:
void evalPhaseStorage_(int phaseIdx)
{
......
......@@ -85,7 +85,7 @@ public:
{
scvfIdx_ = faceIdx;
calculateValues_(problem, element, this->face(), elemVolVars);
calculateValues_(problem, element, elemVolVars);
}
/*!
......@@ -100,10 +100,8 @@ public:
{ return temperatureGrad_; }
protected:
template<class FaceType>
void calculateValues_(const Problem &problem,
const Element &element,
const FaceType &face,
const ElementVolumeVariables &elemVolVars)
{
// calculate temperature gradient using finite element
......@@ -112,13 +110,13 @@ protected:
Vector tmp(0.0);
for (int vertIdx = 0; vertIdx < this->fvGeom_.numVertices; vertIdx++)
{
tmp = face.grad[vertIdx];
tmp = this->face().grad[vertIdx];
tmp *= elemVolVars[vertIdx].temperature();
temperatureGrad_ += tmp;
}
// The spatial parameters calculates the actual heat flux vector
if (face.i != face.j)
if (this->face().i != this->face().j)
problem.spatialParameters().matrixHeatFlux(tmp,
*this,
elemVolVars,
......@@ -130,12 +128,12 @@ protected:
problem.spatialParameters().boundaryMatrixHeatFlux(tmp,
*this,
elemVolVars,
face,
this->face(),
element,
this->fvGeom_);
// project the heat flux vector on the face's normal vector
normalMatrixHeatFlux_ = tmp*face.normal;
normalMatrixHeatFlux_ = tmp*this->face().normal;
}
private:
......
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