Commit f7adc2d1 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

box: fix for dirichlet boundaries

as described in my e-mail from wednesday. Everything in the stable
repository compiles.

also, some compile fixes for the Mp-Nc model where made to some fluid
systems and the 1p2c index structure got an offset template argument
like the other models.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4446 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 56b8036d
......@@ -31,15 +31,16 @@ namespace Dumux
/*!
* \brief The indices for the isothermal single-phase, two-component model.
*/
template <int PVOffset = 0>
struct OnePTwoCIndices
{
// Equation indices
static const int contiEqIdx = 0; //!< continuity equation index
static const int transEqIdx = 1; //!< transport equation index
static const int contiEqIdx = PVOffset + 0; //!< continuity equation index
static const int transEqIdx = PVOffset + 1; //!< transport equation index
// primary variable indices
static const int pressureIdx = 0; //!< pressure
static const int x1Idx = 1; //!< mole fraction of the second component
static const int pressureIdx = PVOffset + 0; //!< pressure
static const int x1Idx = PVOffset + 1; //!< mole fraction of the second component
//in my case the therapeutic agent
};
......
......@@ -61,7 +61,7 @@ SET_TYPE_PROP(BoxOnePTwoC, FluxVariables, OnePTwoCFluxVariables<TypeTag>);
SET_SCALAR_PROP(BoxOnePTwoC, UpwindAlpha, 1.0);
//! Set the indices used by the 1p2c model
SET_TYPE_PROP(BoxOnePTwoC, OnePTwoCIndices, Dumux::OnePTwoCIndices);
SET_TYPE_PROP(BoxOnePTwoC, OnePTwoCIndices, Dumux::OnePTwoCIndices<0>);
//! Set the default phase used by the fluid system to the first one
SET_INT_PROP(BoxOnePTwoC, PhaseIndex, 0);
......
......@@ -140,7 +140,9 @@ public:
{ return heatCapacity_; };
protected:
// this method gets called by the parent class
// this method gets called by the parent class. since this method
// is protected, we are friends with our parent..
friend class TwoPTwoCVolumeVariables<TypeTag>;
void updateTemperature_(const PrimaryVariables &sol,
const Element &element,
const FVElementGeometry &elemGeom,
......
......@@ -39,8 +39,6 @@ class BoxElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTa
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
typedef typename RefElemProp::Container ReferenceElements;
......@@ -50,6 +48,11 @@ class BoxElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTa
enum { dim = GridView::dimension };
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
public:
// copying a the boundary types of an element should be explicitly
// requested
......@@ -61,54 +64,53 @@ public:
* \brief The constructor.
*/
BoxElementBoundaryTypes()
{ }
{
hasDirichlet_ = false;
hasNeumann_ = false;
}
/*!
* \brief Update the boundary types for all vertices of an element.
*/
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvElemGeom)
{
Dune::GeometryType geoType = element.geometry().type();
const ReferenceElement &refElem = ReferenceElements::general(geoType);
int numVerts = element.template count<dim>();
this->resize(numVerts);
for (int i = 0; i < numVerts; ++i)
hasDirichlet_ = false;
hasNeumann_ = false;
for (int i = 0; i < numVerts; ++i) {
(*this)[i].reset();
// evaluate boundary conditions
IntersectionIterator isIt = problem.gridView().ibegin(element);
const IntersectionIterator &endIt = problem.gridView().iend(element);
for (; isIt != endIt; ++isIt) {
// Ignore non- boundary faces.
if (!isIt->boundary())
if (!problem.model().onBoundary(element, i))
continue;
// Set the boundary type for all vertices of the face
int faceIdx = isIt->indexInInside();
int numFaceVerts = refElem.size(faceIdx, 1, dim);
for (int faceVertIdx = 0;
faceVertIdx < numFaceVerts;
faceVertIdx++)
{
int elemVertIdx = refElem.subEntity(faceIdx,
1,
faceVertIdx,
dim);
int boundaryFaceIdx =
fvElemGeom.boundaryFaceIndex(faceIdx,
faceVertIdx);
// set the boundary types
problem.boundaryTypes((*this)[elemVertIdx],
element,
fvElemGeom,
*isIt,
elemVertIdx,
boundaryFaceIdx);
(*this)[elemVertIdx].checkWellPosed();
Valgrind::CheckDefined((*this)[elemVertIdx]);
}
const VertexPointer vptr = element.template subEntity<dim>(i);
problem.boundaryTypes((*this)[i], *vptr);
hasDirichlet_ = hasDirichlet_ or (*this)[i].hasDirichlet();
hasNeumann_ = hasNeumann_ or (*this)[i].hasNeumann();
}
};
/*!
* \brief Returns whether the element has a Dirichlet value.
*/
bool hasDirichlet() const
{ return hasDirichlet_; }
/*!
* \brief Returns whether the element potentially has a Neumann
* boundary segment.
*/
bool hasNeumann() const
{ return hasNeumann_; }
protected:
bool hasDirichlet_;
bool hasNeumann_;
};
} // namespace Dumux
......
......@@ -67,6 +67,8 @@ private:
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
typedef typename RefElemProp::Container ReferenceElements;
......@@ -215,7 +217,10 @@ public:
asImp_().evalVolumeTerms_();
// evaluate the boundary
asImp_().evalBoundary_();
if (bcTypes.hasNeumann())
asImp_().evalNeumann_();
if (bcTypes.hasDirichlet())
asImp_().evalDirichlet_();
#if HAVE_VALGRIND
for (int i=0; i < fvElemGeom_().numVertices; i++)
......@@ -234,7 +239,7 @@ public:
* \brief Returns the local residual for a given sub-control
* volume of the element.
*/
const PrimaryVariables residual(int scvIdx) const
const PrimaryVariables &residual(int scvIdx) const
{ return residual_[scvIdx]; }
/*!
......@@ -251,13 +256,36 @@ protected:
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
void evalBoundary_()
// set the values of the Dirichlet control volumes
void evalDirichlet_()
{
PrimaryVariables tmp;
for (int i = 0; i < fvElemGeom_().numVertices; ++i) {
const BoundaryTypes &bcTypes = bcTypes_(i);
if (! bcTypes.hasDirichlet())
continue;
// ask the problem for the dirichlet values
const VertexPointer vPtr = elem_().template subEntity<dim>(i);
asImp_().problem_().dirichlet(tmp, *vPtr);
// set the dirichlet conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!bcTypes.isDirichlet(eqIdx))
continue;
int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
this->residual_[i][eqIdx] =
curPrimaryVars_(i)[pvIdx] - tmp[pvIdx];
};
};
}
// evaluate the neumann boundary segments
void evalNeumann_()
{
Dune::GeometryType geoType = elem_().geometry().type();
const ReferenceElement &refElem = ReferenceElements::general(geoType);
// evaluate boundary conditions for all intersections of
// the current element
IntersectionIterator isIt = gridView_().ibegin(elem_());
const IntersectionIterator &endIt = gridView_().iend(elem_());
for (; isIt != endIt; ++isIt)
......@@ -278,22 +306,23 @@ protected:
1,
faceVertIdx,
dim);
int boundaryFaceIdx =
fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);
// add the neuman fluxes of a single boundary segment
evalBoundarySegment_(isIt,
elemVertIdx,
boundaryFaceIdx);
// add the residual of all vertices of the boundary
// segment
evalNeumannSegment_(isIt,
elemVertIdx,
boundaryFaceIdx);
}
}
}
// handle boundary conditions for a single sub-control volume face
void evalBoundarySegment_(const IntersectionIterator &isIt,
int scvIdx,
int boundaryFaceIdx)
// handle Neumann boundary conditions for a single sub-control volume face
void evalNeumannSegment_(const IntersectionIterator &isIt,
int scvIdx,
int boundaryFaceIdx)
{
// temporary vector to store the neumann boundary fluxes
PrimaryVariables values(0.0);
......@@ -311,26 +340,6 @@ protected:
Valgrind::CheckDefined(values);
residual_[scvIdx] += values;
}
// deal with dirichlet boundaries
if (bcTypes.hasDirichlet()) {
problem_().dirichlet(values,
elem_(),
fvElemGeom_(),
*isIt,
scvIdx,
boundaryFaceIdx);
Valgrind::CheckDefined(values);
for (int i = 0; i < numEq; ++i) {
if (!bcTypes.isDirichlet(i))
continue;
int pvIdx = bcTypes.eqToDirichletIndex(i);
residual_[scvIdx][i] =
curPrimaryVars_(scvIdx)[pvIdx] - values[pvIdx];
}
}
}
void evalFluxes_()
......
......@@ -109,6 +109,8 @@ public:
void init(Problem &prob)
{
problemPtr_ = &prob;
updateBoundaryIndices_();
int nDofs = asImp_().numDofs();
uCur_.resize(nDofs);
......@@ -569,6 +571,27 @@ public:
const GridView &gridView() const
{ return problem_().gridView(); }
/*!
* \brief Returns true if the vertex with 'globalVertIdx' is
* located on the grid's boundary.
*/
bool onBoundary(int globalVertIdx) const
{ return boundaryIndices_.count(globalVertIdx) > 0; }
/*!
* \brief Returns true if a vertex is located on the grid's
* boundary.
*/
bool onBoundary(const Vertex &vertex) const
{ return onBoundary(vertexMapper().map(vertex)); }
/*!
* \brief Returns true if a vertex is located on the grid's
* boundary.
*/
bool onBoundary(const Element &elem, int vIdx) const
{ return onBoundary(vertexMapper().map(elem, vIdx, dim)); }
protected:
/*!
* \brief A reference to the problem on which the model is applied.
......@@ -644,6 +667,43 @@ protected:
}
}
// find all indices of boundary vertices. for this we need to loop
// over all intersections. if the DUNE grid interface would
// provide a onBoundary() method for entities this could be done
// in a much nicer way (actually this would not be necessary)
void updateBoundaryIndices_()
{
boundaryIndices_.clear();
ElementIterator eIt = gridView_().template begin<0>();
ElementIterator eEndIt = gridView_().template end<0>();
for (; eIt != eEndIt; ++eIt) {
Dune::GeometryType geoType = eIt->geometry().type();
const ReferenceElement &refElem = ReferenceElements::general(geoType);
IntersectionIterator isIt = gridView_().ibegin(*eIt);
IntersectionIterator isEndIt = gridView_().iend(*eIt);
for (; isIt != isEndIt; ++isIt) {
if (!isIt->boundary())
continue;
// add all vertices on the intersection to the set of
// boundary vertices
int faceIdx = isIt->indexInInside();
int numFaceVerts = refElem.size(faceIdx, 1, dim);
for (int faceVertIdx = 0;
faceVertIdx < numFaceVerts;
++faceVertIdx)
{
int elemVertIdx = refElem.subEntity(faceIdx,
1,
faceVertIdx,
dim);
int globalVertIdx = vertexMapper().map(*eIt, elemVertIdx, dim);
boundaryIndices_.insert(globalVertIdx);
}
}
}
}
bool verbose_() const
{ return gridView_().comm().rank() == 0; };
......@@ -656,6 +716,9 @@ protected:
// relations, matxerial laws, etc.
Problem *problemPtr_;
// the set of all indices of vertices on the boundary
std::set<int> boundaryIndices_;
// calculates the local jacobian matrix for a given element
LocalJacobian localJacobian_;
// Linearizes the problem at the current time step using the
......
......@@ -55,8 +55,6 @@ public:
{
evalPoint_ = 0;
primaryVars_ = v.primaryVars_;
Valgrind::CheckDefined(*this);
};
// assignment operator
......@@ -64,7 +62,6 @@ public:
{
evalPoint_ = 0;
primaryVars_ = v.primaryVars_;
Valgrind::CheckDefined(*this);
return *this;
};
......
......@@ -45,10 +45,12 @@ class BoxAssembler
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
enum{dim = GridView::dimension};
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
......@@ -381,6 +383,7 @@ public:
const SolutionVector& residual() const
{ return residual_; }
private:
void resetMatrix_()
{
......
......@@ -77,7 +77,13 @@ public:
{ return Common::molarMass; }
/*!
* \brief Returns the critical temperature [K] of water.
* \brief The acentric factor [] of water.
*/
static Scalar acentricFactor()
{ return Common::acentricFactor; }
/*!
* \brief Returns the critical temperature [K] of water
*/
static Scalar criticalTemperature()
{ return Common::criticalTemperature; }
......@@ -88,6 +94,12 @@ public:
static Scalar criticalPressure()
{ return Common::criticalPressure; }
/*!
* \brief Returns the molar volume [m^3/mol] of water at the critical point
*/
static Scalar criticalMolarVolume()
{ return Common::criticalMolarVolume; }
/*!
* \brief Returns the temperature [K] at water's triple point.
*/
......
......@@ -73,6 +73,12 @@ public:
//! Critical pressure of water [Pa]
static const Scalar criticalPressure = 22.064e6;
//! Critical molar volume of water [m^3/mol]
static const Scalar criticalMolarVolume = molarMass/322.0;
//! The acentric factor of water []
static const Scalar acentricFactor = 0.344;
//! Density of water at the critical point [kg/m^3]
static const Scalar criticalDensity = 322;
......
......@@ -101,13 +101,17 @@ public:
return EffLaw::krn(params, SwToSwe(params, Sw));
}
// convert an absolute wetting saturation to an effective one
/*!
* \brief Convert an absolute wetting saturation to an effective one.
*/
static Scalar SwToSwe(const Params &params, Scalar Sw)
{
return (Sw - params.Swr())/(1 - params.Swr() - params.Snr());
}
// convert an absolute wetting saturation to an effective one
/*!
* \brief convert an absolute wetting saturation to an effective one
*/
static Scalar SnToSne(const Params &params, Scalar Sn)
{
return (Sn - params.Snr())/(1 - params.Swr() - params.Snr());
......
......@@ -60,7 +60,7 @@ public:
\f]
* \param Sw Effective saturation of of the wetting phase \f$\overline{S}_w\f$
*/
static Scalar pC(const Params &params, Scalar Sw)
static Scalar pC(const Params &params, Scalar Swe)
{
// retrieve the low and the high threshold saturations for the
// unregularized capillary pressure curve from the parameters
......@@ -72,16 +72,16 @@ public:
// newton solver (if the derivative is calculated numerically)
// in order to get the saturation moving to the right
// direction if it temporarily is in an 'illegal' range.
if (Sw < SwThLow) {
return VanGenuchten::pC(params, SwThLow) + mLow_(params)*(Sw - SwThLow);
if (Swe < SwThLow) {
return VanGenuchten::pC(params, SwThLow) + mLow_(params)*(Swe - SwThLow);
}
else if (Sw > SwThHigh) {
return VanGenuchten::pC(params, SwThHigh) + mHigh_(params)*(Sw - SwThHigh);
else if (Swe > SwThHigh) {
return VanGenuchten::pC(params, SwThHigh) + mHigh_(params)*(Swe - SwThHigh);
}
// if the effective saturation is in an 'reasonable'
// range, we use the real van genuchten law...
return VanGenuchten::pC(params, Sw);
return VanGenuchten::pC(params, Swe);
}
/*!
......
......@@ -73,6 +73,18 @@ public:
static void init()
{ Components::init(); }
/*!
* \brief Return the human readable name of a phase
*/
static const char *phaseName(int phaseIdx)
{
switch (phaseIdx) {
case lPhaseIdx: return "l";
case gPhaseIdx: return "g";
};
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
* \brief Return the human readable name of a component
*/
......
......@@ -484,7 +484,7 @@ public:
if (error_ < 10*tolerance_)
reassembleTol = tolerance_/5;
this->model_().jacobianAssembler().computeColors(reassembleTol);
}
}
}
/*!
......
......@@ -165,17 +165,11 @@ public:
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*/
void boundaryTypes(BoundaryTypes &values,
const Element &element,
const FVElementGeometry &fvElemGeom,
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const
void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
{
double eps = 1.0e-3;
const GlobalPosition &globalPos
= fvElemGeom.boundaryFace[boundaryFaceIdx].ipGlobal;
const GlobalPosition globalPos = vertex.geometry().center();
double eps = 1.0e-3;
if (globalPos[dim-1] < eps || globalPos[dim-1] > this->bboxMax()[dim-1] - eps)
values.setAllDirichlet();
else
......@@ -188,15 +182,10 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichlet(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvElemGeom,
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const
void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
{
double eps = 1.0e-3;
const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
const GlobalPosition globalPos = vertex.geometry().center();
if (globalPos[dim-1] < eps) {
values[pressureIdx] = 2.0e+5;
......
......@@ -205,12 +205,7 @@ public:
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*/
void boundaryTypes(BoundaryTypes &values,
const Element &element,
const FVElementGeometry &fvElemGeom,
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const
void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
{
values.setAllDirichlet();
}
......@@ -221,15 +216,9 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichlet(PrimaryVariables &values,
const Element &element,
const FVElementGeometry &fvElemGeom,
const Intersection &is,
int scvIdx,
int boundaryFaceIdx) const