diff --git a/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh b/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh
index 7b4a16145a0e8e33041902ad4108f87734f34d35..77281bb0c05466566d6007361420c4df69c649a3 100644
--- a/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh
+++ b/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh
@@ -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_;
diff --git a/dumux/boxmodels/1p2c/1p2clocalresidual.hh b/dumux/boxmodels/1p2c/1p2clocalresidual.hh
index a5789a8ced0fe7a0c37afd676ff06ac13dff6bed..7e01b354f811d4bea17c75fdd0c89c7fc91bbf07 100644
--- a/dumux/boxmodels/1p2c/1p2clocalresidual.hh
+++ b/dumux/boxmodels/1p2c/1p2clocalresidual.hh
@@ -1,4 +1,5 @@
/*****************************************************************************
+ * 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); }
diff --git a/dumux/boxmodels/1p2c/1p2cmodel.hh b/dumux/boxmodels/1p2c/1p2cmodel.hh
index cbbd4263949f223485c7781f68207533c5193fcc..bda4e336c491bc12420c6551ba05508cf26d185b 100644
--- a/dumux/boxmodels/1p2c/1p2cmodel.hh
+++ b/dumux/boxmodels/1p2c/1p2cmodel.hh
@@ -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