Commit 74ebe0e1 authored by Katherina Baber's avatar Katherina Baber
Browse files

moved all functions relevant for coupling to 1p2cCoupling model


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6769 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 1b425287
......@@ -88,6 +88,8 @@ public:
{
boundaryFace_ = &fvElemGeom_.boundaryFace[boundaryFaceIdx];
pressureAtBIP_ = Scalar(0);
massFractionAtBIP_ = Scalar(0);
densityAtIP_ = Scalar(0);
viscosityAtIP_ = Scalar(0);
molarDensityAtIP_ = Scalar(0);
......@@ -121,6 +123,20 @@ public:
Scalar porousDiffCoeff() const
{ return porousDiffCoeff_; };
/*!
* \brief Return pressure \f$\mathrm{[Pa]}\f$ of a phase at the integration
* point.
*/
Scalar pressureAtBIP() const
{ return pressureAtBIP_; }
/*!
* \brief Return massFraction \f$\mathrm{[-]}\f$ of component 1 at the integration
* point.
*/
Scalar massFractionAtBIP() const
{ return massFractionAtBIP_; }
/*!
* \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
* point.
......@@ -231,12 +247,17 @@ protected:
// FE gradient at vertex idx
const ScalarGradient& feGrad = boundaryFace_->grad[idx];
pressureAtBIP_ += elemDat[idx].pressure() *
boundaryFace_->shapeValue[idx];
// compute sum of pressure gradients for each phase
// the pressure gradient
tmp = feGrad;
tmp *= elemDat[idx].pressure();
potentialGrad_ += tmp;
massFractionAtBIP_ += elemDat[idx].massFrac(comp1Idx) *
boundaryFace_->shapeValue[idx];
// the concentration gradient of the non-wetting
// component in the wetting phase
tmp = feGrad;
......@@ -251,6 +272,7 @@ protected:
tmp *= elemDat[idx].fluidState().massFrac(phaseIdx, comp1Idx);
massFracGrad_ += tmp;
densityAtIP_ += elemDat[idx].density()*boundaryFace_->shapeValue[idx];
viscosityAtIP_ += elemDat[idx].viscosity()*boundaryFace_->shapeValue[idx];
molarDensityAtIP_ += elemDat[idx].molarDensity()*boundaryFace_->shapeValue[idx];
......@@ -285,7 +307,9 @@ protected:
ScalarGradient moleFracGrad_;
ScalarGradient massFracGrad_;
// density of each face at the integration point
// quanitities at the integration point
Scalar pressureAtBIP_;
Scalar massFractionAtBIP_;
Scalar densityAtIP_;
Scalar viscosityAtIP_;
Scalar molarDensityAtIP_;
......
/*****************************************************************************
* Copyright (C) 2011 by Katherina Baber *
* Copyright (C) 2009 by Karin Erbertseder *
* Copyright (C) 2009 by Andreas Lauser *
* Copyright (C) 2008 by Bernd Flemisch *
......@@ -60,7 +61,6 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
typedef BoxLocalResidual<TypeTag> ParentType;
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;
......@@ -176,12 +176,12 @@ public:
}
/*!
* \brief Evaluates the advective mass flux of all components over
* a face of a subcontrol volume.
*
* \param flux The advective flux over the sub-control-volume face for each component
* \param vars The flux variables at the current SCV
*/
* \brief Evaluates the advective mass flux of all components over
* a face of a subcontrol volume.
*
* \param flux The advective flux over the sub-control-volume face for each component
* \param vars The flux variables at the current SCV
*/
void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
////////
......@@ -189,10 +189,10 @@ public:
////////
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up =
// of the current phase
const VolumeVariables &up =
this->curVolVars_(fluxVars.upstreamIdx());
const VolumeVariables &dn =
const VolumeVariables &dn =
this->curVolVars_(fluxVars.downstreamIdx());
if(!useMoles)
......@@ -233,12 +233,12 @@ public:
}
/*!
* \brief Adds the diffusive mass flux of all components over
* a face of a subcontrol volume.
*
* \param flux The diffusive flux over the sub-control-volume face for each component
* \param vars The flux variables at the current SCV
*/
* \brief Adds the diffusive mass flux of all components over
* a face of a subcontrol volume.
*
* \param flux The diffusive flux over the sub-control-volume face for each component
* \param vars The flux variables at the current SCV
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
Scalar tmp(0);
......@@ -273,6 +273,7 @@ public:
flux[transEqIdx] += tmp;
}
}
/*!
* \brief Calculate the source term of the equation
* \param q The source/sink in the SCV for each component
......@@ -287,239 +288,142 @@ public:
localVertexIdx);
}
/*!
* \brief Evaluate Neuman, Outflow and Dirichlet conditions.
*
*/
void evalBoundary_()
{
ParentType::evalBoundary_();
if (this->bcTypes_().hasNeumann())
this->evalNeumann_();
if (this->bcTypes_().hasOutflow())
{
evalOutflow_();
}
evalOutflow_();
typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
const ReferenceElement &refElem = ReferenceElements::general(this->elem_().geometry().type());
if (this->bcTypes_().hasDirichlet())
this->evalDirichlet_();
}
// loop over vertices of the element
for (int idx = 0; idx < this->fvElemGeom_().numVertices; idx++)
protected:
/*!
* \brief Add all Outflow boundary conditions to the local
* residual.
*/
void evalOutflow_()
{
Dune::GeometryType geoType = this->elem_().geometry().type();
typedef typename Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
typedef typename Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
const ReferenceElement &refElem = ReferenceElements::general(geoType);
IntersectionIterator isIt = this->gridView_().ibegin(this->elem_());
const IntersectionIterator &endIt = this->gridView_().iend(this->elem_());
for (; isIt != endIt; ++isIt)
{
// consider only SCVs on the boundary
if (this->fvElemGeom_().subContVol[idx].inner)
// handle only faces on the boundary
if (!isIt->boundary())
continue;
int numberOfOuterFaces = 0;
// evaluate boundary conditions for the intersections of
// the current element
IntersectionIterator isIt = this->gridView_().ibegin(this->elem_());
const IntersectionIterator &endIt = this->gridView_().iend(this->elem_());
for (; isIt != endIt; ++isIt)
// Assemble the boundary for all vertices of the current
// face
int faceIdx = isIt->indexInInside();
int numFaceVerts = refElem.size(faceIdx, 1, dim);
for (int faceVertIdx = 0;
faceVertIdx < numFaceVerts;
++faceVertIdx)
{
// handle only intersections on the boundary
if (!isIt->boundary())
continue;
// assemble the boundary for all vertices of the current face
const int faceIdx = isIt->indexInInside();
const int numFaceVertices = refElem.size(faceIdx, 1, dim);
// loop over the single vertices on the current face
for (int faceVertIdx = 0; faceVertIdx < numFaceVertices; ++faceVertIdx)
{
const int elemVertIdx = refElem.subEntity(faceIdx, 1, faceVertIdx, dim);
// only evaluate, if we consider the same face vertex as in the outer
// loop over the element vertices
if (elemVertIdx != idx)
continue;
if (boundaryHasCoupling_(this->bcTypes_(idx)))
{
const int boundaryFaceIdx = this->fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);
asImp_()->evalCouplingVertex_(isIt, elemVertIdx, boundaryFaceIdx);
}
// count the number of outer faces to determine, if we are on
// a corner point and if an interpolation should be done
numberOfOuterFaces++;
}
} // end loop over intersections
if (numberOfOuterFaces == 2 && this->fvElemGeom_().numVertices == 4)
// do interpolation at corners of the grid
interpolateCornerPoints_(idx);
} // end loop over vertices
int elemVertIdx = refElem.subEntity(faceIdx,
1,
faceVertIdx,
dim);
int boundaryFaceIdx =
this->fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);
// add the residual of all vertices of the boundary
// segment
evalOutflowSegment_(isIt,
elemVertIdx,
boundaryFaceIdx);
}
}
}
protected:
/*!
* \brief Add all Outflow boundary conditions to the local
* residual.
*/
void evalOutflow_()
/*!
* \brief Add Outflow boundary conditions for a single sub-control
* volume face to the local residual.
*/
void evalOutflowSegment_(const IntersectionIterator &isIt,
int scvIdx,
int boundaryFaceIdx)
{
// temporary vector to store the neumann boundary fluxes
PrimaryVariables flux(0.0);
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
// deal with neumann boundaries
if (bcTypes.hasOutflow())
{
Dune::GeometryType geoType = this->elem_().geometry().type();
const BoundaryVariables boundaryVars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
boundaryFaceIdx,
this->curVolVars_(),
scvIdx);
typedef typename Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
typedef typename Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
const ReferenceElement &refElem = ReferenceElements::general(geoType);
const VolumeVariables& vertVars = this->curVolVars_()[scvIdx];
IntersectionIterator isIt = this->gridView_().ibegin(this->elem_());
const IntersectionIterator &endIt = this->gridView_().iend(this->elem_());
for (; isIt != endIt; ++isIt)
// mass balance
if (bcTypes.isOutflow(contiEqIdx))
{
// handle only faces on the boundary
if (!isIt->boundary())
continue;
// Assemble the boundary for all vertices of the current
// face
int faceIdx = isIt->indexInInside();
int numFaceVerts = refElem.size(faceIdx, 1, dim);
for (int faceVertIdx = 0;
faceVertIdx < numFaceVerts;
++faceVertIdx)
if(!useMoles) //use massfractions
{
int elemVertIdx = refElem.subEntity(faceIdx,
1,
faceVertIdx,
dim);
int boundaryFaceIdx =
this->fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);
// add the residual of all vertices of the boundary
// segment
evalOutflowSegment_(isIt,
elemVertIdx,
boundaryFaceIdx);
flux[contiEqIdx] += boundaryVars.KmvpNormal()*vertVars.density()/vertVars.viscosity();
}
else //use molefractions
{
flux[contiEqIdx] += boundaryVars.KmvpNormal()*vertVars.molarDensity()/vertVars.viscosity();
}
}
}
/*!
* \brief Add Outflow boundary conditions for a single sub-control
* volume face to the local residual.
*/
void evalOutflowSegment_(const IntersectionIterator &isIt,
int scvIdx,
int boundaryFaceIdx)
{
// temporary vector to store the neumann boundary fluxes
PrimaryVariables flux(0.0);
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
// deal with neumann boundaries
if (bcTypes.hasOutflow())
{
const BoundaryVariables boundaryVars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
boundaryFaceIdx,
this->curVolVars_(),
scvIdx);
const VolumeVariables& vertVars = this->curVolVars_()[scvIdx];
// mass balance
if (bcTypes.isOutflow(contiEqIdx))
{
if(!useMoles) //use massfractions
{
flux[contiEqIdx] += boundaryVars.KmvpNormal()*vertVars.density()/vertVars.viscosity();
}
else //use molefractions
{
flux[contiEqIdx] += boundaryVars.KmvpNormal()*vertVars.molarDensity()/vertVars.viscosity();
}
}
// component transport
if (bcTypes.isOutflow(transEqIdx))
{
if(!useMoles)//use massfractions
{
// advective flux
flux[transEqIdx]+= boundaryVars.KmvpNormal()*vertVars.density()/vertVars.viscosity()
*vertVars.fluidState().massFrac(phaseIdx, comp1Idx);
// diffusive flux of comp1 component in phase0
Scalar tmp = 0;
for (int i = 0; i < Vector::size; ++ i)
tmp += boundaryVars.massFracGrad(comp1Idx)[i]*boundaryVars.boundaryFace().normal[i];
tmp *= -1;
tmp *= boundaryVars.porousDiffCoeff()*boundaryVars.densityAtIP();
flux[transEqIdx] += tmp;//* FluidSystem::molarMass(comp1Idx);
}
else //use molefractions
{
// advective flux
flux[transEqIdx]+= boundaryVars.KmvpNormal()*vertVars.molarDensity()/vertVars.viscosity()
*vertVars.fluidState().moleFrac(phaseIdx, comp1Idx);
// diffusive flux of comp1 component in phase0
Scalar tmp = 0;
for (int i = 0; i < Vector::size; ++ i)
tmp += boundaryVars.moleFracGrad(comp1Idx)[i]*boundaryVars.boundaryFace().normal[i];
tmp *= -1;
tmp *= boundaryVars.porousDiffCoeff()*boundaryVars.molarDensityAtIP();
flux[transEqIdx] += tmp;
}
}
this->residual_[scvIdx] += flux;
}
}
void evalCouplingVertex_(const IntersectionIterator &isIt,
const int scvIdx,
const int boundaryFaceIdx)
{
const VolumeVariables &volVars = this->curVolVars_()[scvIdx];
// set pressure as part of the momentum coupling
if (this->bcTypes_(scvIdx).isCouplingOutflow(contiEqIdx))
{
this->residual_[scvIdx][contiEqIdx] = volVars.pressure();
}
if (this->bcTypes_(scvIdx).isCouplingOutflow(transEqIdx))
{
if(!useMoles)
this->residual_[scvIdx][transEqIdx] = volVars.fluidState().massFrac(phaseIdx, comp1Idx);
else
this->residual_[scvIdx][transEqIdx] = volVars.fluidState().moleFrac(phaseIdx, comp1Idx);
}
}
/*!
* \brief Interpolate the corner points of the grid.
*/
void interpolateCornerPoints_(const int idx)
{
// for (int eqnIdx=0; eqnIdx<numEq; ++eqnIdx)
// {
// //do not interpolate in case of Beavers-Joseph or Dirichlet condition for coupling interface
// if (boundaryHasCoupling_(this->bcTypes_(idx)))
// continue;
//
// this->residual_[idx][eqnIdx] = this->curPrimaryVars_(0)[eqnIdx]+this->curPrimaryVars_(3)[eqnIdx]
// -this->curPrimaryVars_(1)[eqnIdx]-this->curPrimaryVars_(2)[eqnIdx];
// }
}
/*!
* \brief Check if one of the boundary conditions is coupling.
*/
bool boundaryHasCoupling_(const BoundaryTypes& bcTypes) const
{
bool hasCoupling = false;
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
if (bcTypes.isCouplingInflow(eqIdx) || bcTypes.isCouplingOutflow(eqIdx))
hasCoupling = true;
return hasCoupling;
}
// component transport
if (bcTypes.isOutflow(transEqIdx))
{
if(!useMoles)//use massfractions
{
// advective flux
flux[transEqIdx]+= boundaryVars.KmvpNormal()*vertVars.density()/vertVars.viscosity()
*vertVars.fluidState().massFrac(phaseIdx, comp1Idx);
// diffusive flux of comp1 component in phase0
Scalar tmp = 0;
for (int i = 0; i < Vector::size; ++ i)
tmp += boundaryVars.massFracGrad(comp1Idx)[i]*boundaryVars.boundaryFace().normal[i];
tmp *= -1;
tmp *= boundaryVars.porousDiffCoeff()*boundaryVars.densityAtIP();
flux[transEqIdx] += tmp;//* FluidSystem::molarMass(comp1Idx);
}
else //use molefractions
{
// advective flux
flux[transEqIdx]+= boundaryVars.KmvpNormal()*vertVars.molarDensity()/vertVars.viscosity()
*vertVars.fluidState().moleFrac(phaseIdx, comp1Idx);
// diffusive flux of comp1 component in phase0
Scalar tmp = 0;
for (int i = 0; i < Vector::size; ++ i)
tmp += boundaryVars.moleFracGrad(comp1Idx)[i]*boundaryVars.boundaryFace().normal[i];
tmp *= -1;
tmp *= boundaryVars.porousDiffCoeff()*boundaryVars.molarDensityAtIP();
flux[transEqIdx] += tmp;
}
}
this->residual_[scvIdx] += flux;
}
}
Implementation *asImp_()
{ return static_cast<Implementation *> (this); }
......
......@@ -84,10 +84,8 @@ class OnePTwoCBoxModel : public BoxModel<TypeTag>
typedef OnePTwoCBoxModel<TypeTag> ThisType;
typedef BoxModel<TypeTag> ParentType;
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::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
......@@ -95,8 +93,6 @@ class OnePTwoCBoxModel : public BoxModel<TypeTag>
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;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
enum {
......@@ -104,13 +100,10 @@ class OnePTwoCBoxModel : public BoxModel<TypeTag>
dimWorld = GridView::dimensionworld
};
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
const Scalar scale_;
public:
OnePTwoCBoxModel():
scale_(GET_PROP_VALUE(TypeTag, PTAG(Scaling)))
OnePTwoCBoxModel()
{
// retrieve the upwind weight for the mass conservation equations. Use the value
// specified via the property system as default, and overwrite
......@@ -194,14 +187,14 @@ public:
i,
false);
pressure[globalIdx] = volVars.pressure()*scale_;
delp[globalIdx] = volVars.pressure()*scale_ - 1e5;
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()*scale_*scale_*scale_;
mu[globalIdx] = volVars.viscosity()*scale_;
rho[globalIdx] = volVars.density();
mu[globalIdx] = volVars.viscosity();
delFrac[globalIdx] = volVars.massFrac(1)-volVars.moleFrac(1);
};
......@@ -302,16 +295,13 @@ public:
//use vertiacl faces for vx and horizontal faces for vy calculation
velocityX[i] /= boxSurface[i][0];
velocityX[i] /= scale_;
if (dim >= 2)
{
velocityY[i] /= boxSurface[i][1];
velocityY[i] /= scale_;
}
if (dim == 3)
{
velocityZ[i] /= boxSurface[i][2];
velocityZ[i] /= scale_;
}
}
#endif
......
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