From f7adc2d1d73d028bc52f5492c04a642810bebf68 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 15 Oct 2010 13:37:28 +0000
Subject: [PATCH] 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
---
dumux/boxmodels/1p2c/1p2cindices.hh | 9 ++-
dumux/boxmodels/1p2c/1p2cpropertydefaults.hh | 2 +-
.../boxmodels/2p2cni/2p2cnivolumevariables.hh | 4 +-
.../common/boxelementboundarytypes.hh | 76 +++++++++---------
dumux/boxmodels/common/boxlocalresidual.hh | 77 +++++++++++--------
dumux/boxmodels/common/boxmodel.hh | 63 +++++++++++++++
dumux/boxmodels/common/boxvolumevariables.hh | 3 -
dumux/boxmodels/common/pdelabboxassembler.hh | 3 +
dumux/material/components/h2o.hh | 14 +++-
dumux/material/components/iapws/common.hh | 6 ++
.../fluidmatrixinteractions/2p/efftoabslaw.hh | 8 +-
.../2p/regularizedvangenuchten.hh | 12 +--
dumux/material/fluidsystems/h2o_n2_system.hh | 12 +++
dumux/nonlinear/newtoncontroller.hh | 2 +-
test/boxmodels/1p/1ptestproblem.hh | 21 ++---
test/boxmodels/1p2c/tissue_tumor_problem.hh | 17 +---
test/boxmodels/2p/lensproblem.hh | 48 +++++-------
test/boxmodels/2p2c/injectionproblem.hh | 30 ++------
test/boxmodels/2p2cni/waterairproblem.hh | 35 ++-------
test/boxmodels/2pni/injectionproblem2pni.hh | 30 ++------
.../boxmodels/richards/richardslensproblem.hh | 48 ++++++------
.../richards/richardslensspatialparameters.hh | 7 +-
tutorial/tutorialproblem_coupled.hh | 16 +---
23 files changed, 279 insertions(+), 264 deletions(-)
diff --git a/dumux/boxmodels/1p2c/1p2cindices.hh b/dumux/boxmodels/1p2c/1p2cindices.hh
index 8c47e7d844..0ca7145536 100644
--- a/dumux/boxmodels/1p2c/1p2cindices.hh
+++ b/dumux/boxmodels/1p2c/1p2cindices.hh
@@ -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
};
diff --git a/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh b/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
index 393bf0dce4..a423eaf5a7 100644
--- a/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
+++ b/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
@@ -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);
diff --git a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
index 07ed08533b..0cb8b887db 100644
--- a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
@@ -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,
diff --git a/dumux/boxmodels/common/boxelementboundarytypes.hh b/dumux/boxmodels/common/boxelementboundarytypes.hh
index ecf24523b4..a24938b342 100644
--- a/dumux/boxmodels/common/boxelementboundarytypes.hh
+++ b/dumux/boxmodels/common/boxelementboundarytypes.hh
@@ -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
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index 08fc4ac55f..d2e72dba57 100644
--- a/dumux/boxmodels/common/boxlocalresidual.hh
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -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_()
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
index 86a6da17c2..7ee2d0fa82 100644
--- a/dumux/boxmodels/common/boxmodel.hh
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -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
diff --git a/dumux/boxmodels/common/boxvolumevariables.hh b/dumux/boxmodels/common/boxvolumevariables.hh
index 36c7924dd3..1dd6c8405d 100644
--- a/dumux/boxmodels/common/boxvolumevariables.hh
+++ b/dumux/boxmodels/common/boxvolumevariables.hh
@@ -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;
};
diff --git a/dumux/boxmodels/common/pdelabboxassembler.hh b/dumux/boxmodels/common/pdelabboxassembler.hh
index f3f87d80eb..bb5d5aa2b2 100644
--- a/dumux/boxmodels/common/pdelabboxassembler.hh
+++ b/dumux/boxmodels/common/pdelabboxassembler.hh
@@ -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_()
{
diff --git a/dumux/material/components/h2o.hh b/dumux/material/components/h2o.hh
index 9b47c29469..ab73c4e8b2 100644
--- a/dumux/material/components/h2o.hh
+++ b/dumux/material/components/h2o.hh
@@ -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.
*/
diff --git a/dumux/material/components/iapws/common.hh b/dumux/material/components/iapws/common.hh
index 67836c1c62..bcf5d6bfd6 100644
--- a/dumux/material/components/iapws/common.hh
+++ b/dumux/material/components/iapws/common.hh
@@ -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;
diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
index 20b43f3ef3..2fb4b9e3b7 100644
--- a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
@@ -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 ¶ms, 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 ¶ms, Scalar Sn)
{
return (Sn - params.Snr())/(1 - params.Swr() - params.Snr());
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
index b63f48ecb0..cc1905d816 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
@@ -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 ¶ms, Scalar Sw)
+ static Scalar pC(const Params ¶ms, 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);
}
/*!
diff --git a/dumux/material/fluidsystems/h2o_n2_system.hh b/dumux/material/fluidsystems/h2o_n2_system.hh
index 8cd8fd3763..5c9143ec7f 100644
--- a/dumux/material/fluidsystems/h2o_n2_system.hh
+++ b/dumux/material/fluidsystems/h2o_n2_system.hh
@@ -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
*/
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index cee80557a4..ecc43f909d 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -484,7 +484,7 @@ public:
if (error_ < 10*tolerance_)
reassembleTol = tolerance_/5;
this->model_().jacobianAssembler().computeColors(reassembleTol);
- }
+ }
}
/*!
diff --git a/test/boxmodels/1p/1ptestproblem.hh b/test/boxmodels/1p/1ptestproblem.hh
index 29bbf82157..f77bacbcb6 100644
--- a/test/boxmodels/1p/1ptestproblem.hh
+++ b/test/boxmodels/1p/1ptestproblem.hh
@@ -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;
diff --git a/test/boxmodels/1p2c/tissue_tumor_problem.hh b/test/boxmodels/1p2c/tissue_tumor_problem.hh
index c0910a13fe..b80077bf86 100644
--- a/test/boxmodels/1p2c/tissue_tumor_problem.hh
+++ b/test/boxmodels/1p2c/tissue_tumor_problem.hh
@@ -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
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
{
- const GlobalPosition &globalPos
- = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
initial_(values, globalPos);
}
diff --git a/test/boxmodels/2p/lensproblem.hh b/test/boxmodels/2p/lensproblem.hh
index dfb92bc49d..a8995efd46 100644
--- a/test/boxmodels/2p/lensproblem.hh
+++ b/test/boxmodels/2p/lensproblem.hh
@@ -109,14 +109,17 @@ SET_PROP(LensProblem, SpatialParameters)
};
// Enable the time step ramp up inside the newton method?
-SET_BOOL_PROP(LensProblem, EnableTimeStepRampUp, true);
+SET_BOOL_PROP(LensProblem, EnableTimeStepRampUp, false);
// Enable partial reassembly of the jacobian matrix?
-SET_BOOL_PROP(LensProblem, EnablePartialReassemble, true);
+SET_BOOL_PROP(LensProblem, EnablePartialReassemble, false);
// Enable reuse of jacobian matrices?
SET_BOOL_PROP(LensProblem, EnableJacobianRecycling, false);
+// Write the solutions of individual newton iterations?
+SET_BOOL_PROP(LensProblem, NewtonWriteConvergence, false);
+
// Enable gravity
SET_BOOL_PROP(LensProblem, EnableGravity, true);
}
@@ -217,6 +220,7 @@ public:
const GlobalPosition &lensUpperRight)
: ParentType(timeManager, gridView)
{
+ temperature_ = 273.15 + 20; // -> 20°C
this->spatialParameters().setLensCoords(lensLowerLeft, lensUpperRight);
}
@@ -261,7 +265,7 @@ public:
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
- return 273.15 + 10; // -> 10°C
+ return temperature_;
};
// \}
@@ -273,24 +277,16 @@ public:
/*!
* \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
+ * used for which equation on a given boundary control volume.
*
* \param values The boundary types for the conservation equations
- * \param element The finite element
- * \param fvElemGeom The finite-volume geometry in the box scheme
- * \param is The intersection between element and boundary
- * \param scvIdx The local vertex index
- * \param boundaryFaceIdx The index of the boundary face
+ * \param vertex The vertex on the boundary for which the
+ * conditions needs to be specified
*/
void boundaryTypes(BoundaryTypes &values,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
- const Intersection &is,
- int scvIdx,
- int boundaryFaceIdx) const
+ const Vertex &vertex) const
{
- const GlobalPosition &globalPos
- = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) {
values.setAllDirichlet();
@@ -304,30 +300,21 @@ public:
/*!
* \brief Evaluate the boundary conditions for a dirichlet
- * boundary segment.
+ * control volume.
*
* \param values The dirichlet values for the primary variables
- * \param element The finite element
- * \param fvElemGeom The finite-volume geometry in the box scheme
- * \param is The intersection between element and boundary
- * \param scvIdx The local vertex index
- * \param boundaryFaceIdx The index of the boundary face
+ * \param vertex The vertex representing the "half volume on the boundary"
*
* 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
+ const Vertex &vertex) const
{
- const GlobalPosition &globalPos
- = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
Scalar densityW = FluidSystem::componentDensity(wPhaseIdx,
wPhaseIdx,
- temperature(element, fvElemGeom, scvIdx),
+ temperature_,
1e5);
if (onLeftBoundary_(globalPos))
@@ -459,6 +446,7 @@ private:
return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
}
+ Scalar temperature_;
static const Scalar eps_ = 3e-6;
};
} //end namespace
diff --git a/test/boxmodels/2p2c/injectionproblem.hh b/test/boxmodels/2p2c/injectionproblem.hh
index 1b346e12f5..78ffbcc712 100644
--- a/test/boxmodels/2p2c/injectionproblem.hh
+++ b/test/boxmodels/2p2c/injectionproblem.hh
@@ -223,20 +223,11 @@ public:
* used for which equation on a given boundary segment.
*
* \param values The boundary types for the conservation equations
- * \param element The finite element
- * \param fvElemGeom The finite-volume geometry in the box scheme
- * \param is The intersection between element and boundary
- * \param scvIdx The local vertex index
- * \param boundaryFaceIdx The index of the boundary face
+ * \param vertex The vertex for which the boundary type is set
*/
- 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
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
if (globalPos[0] < eps_)
values.setAllDirichlet();
@@ -249,22 +240,13 @@ public:
* boundary segment.
*
* \param values The dirichlet values for the primary variables
- * \param element The finite element
- * \param fvElemGeom The finite-volume geometry in the box scheme
- * \param is The intersection between element and boundary
- * \param scvIdx The local vertex index
- * \param boundaryFaceIdx The index of the boundary face
+ * \param vertex The vertex for which the boundary type is set
*
* 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
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
initial_(values, globalPos);
}
diff --git a/test/boxmodels/2p2cni/waterairproblem.hh b/test/boxmodels/2p2cni/waterairproblem.hh
index 9f43629d4d..f1db731c83 100644
--- a/test/boxmodels/2p2cni/waterairproblem.hh
+++ b/test/boxmodels/2p2cni/waterairproblem.hh
@@ -226,20 +226,11 @@ public:
* used for which equation on a given boundary segment.
*
* \param values The boundary types for the conservation equations
- * \param element The finite element
- * \param fvElemGeom The finite-volume geometry in the box scheme
- * \param is The intersection between element and boundary
- * \param scvIdx The local vertex index
- * \param boundaryFaceIdx The index of the boundary face
+ * \param vertex The vertex for which the boundary type is set
*/
- 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
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
if(globalPos[0] > 40 - eps_ || globalPos[0] < eps_)
values.setAllDirichlet();
@@ -256,26 +247,14 @@ public:
* boundary segment.
*
* \param values The dirichlet values for the primary variables
- * \param element The finite element
- * \param fvElemGeom The finite-volume geometry in the box scheme
- * \param is The intersection between element and boundary
- * \param scvIdx The local vertex index
- * \param boundaryFaceIdx The index of the boundary face
+ * \param vertex The vertex for which the boundary type is set
*
* 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
{
- const GlobalPosition &globalPos
- = element.geometry().corner(scvIdx);
- // const LocalPosition &localPos
- // = DomainTraits::referenceElement(element.geometry().type()).position(dim,scvIdx);
-
+ const GlobalPosition globalPos = vertex.geometry().center();
+
initial_(values, globalPos);
}
diff --git a/test/boxmodels/2pni/injectionproblem2pni.hh b/test/boxmodels/2pni/injectionproblem2pni.hh
index 181b5609c3..f3a395a5bc 100644
--- a/test/boxmodels/2pni/injectionproblem2pni.hh
+++ b/test/boxmodels/2pni/injectionproblem2pni.hh
@@ -256,20 +256,11 @@ public:
* used for which equation on a given boundary segment.
*
* \param values The boundary types for the conservation equations
- * \param element The finite element
- * \param fvElemGeom The finite-volume geometry in the box scheme
- * \param is The intersection between element and boundary
- * \param scvIdx The local vertex index
- * \param boundaryFaceIdx The index of the boundary face
+ * \param vertex The vertex for which the boundary type is set
*/
- 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
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
if (globalPos[0] < eps_)
values.setAllDirichlet();
@@ -288,22 +279,13 @@ public:
* boundary segment.
*
* \param values The dirichlet values for the primary variables
- * \param element The finite element
- * \param fvElemGeom The finite-volume geometry in the box scheme
- * \param is The intersection between element and boundary
- * \param scvIdx The local vertex index
- * \param boundaryFaceIdx The index of the boundary face
+ * \param vertex The vertex for which the boundary type is set
*
* 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
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
Scalar densityW = 1000.0;
values[pressureIdx] = 1e5 + (depthBOR_ - globalPos[1])*densityW*9.81;
diff --git a/test/boxmodels/richards/richardslensproblem.hh b/test/boxmodels/richards/richardslensproblem.hh
index 80763d20db..cff1468081 100644
--- a/test/boxmodels/richards/richardslensproblem.hh
+++ b/test/boxmodels/richards/richardslensproblem.hh
@@ -121,12 +121,13 @@ class RichardsLensProblem : public RichardsBoxProblem<TypeTag>
contiEqIdx = Indices::contiEqIdx,
// Grid and world dimension
- dimWorld = GridView::dimensionworld,
+ dim = GridView::dimensionworld,
};
typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::Intersection Intersection;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
public:
RichardsLensProblem(TimeManager &timeManager,
@@ -185,16 +186,13 @@ public:
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
+ *
+ * \param values The boundary types for the conservation equations
+ * \param vertex The vertex for which the boundary type is set
*/
- 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
{
- const GlobalPosition &globalPos
- = element.geometry().corner(scvIdx);
+ const GlobalPosition globalPos = vertex.geometry().center();
if (onLeftBoundary_(globalPos) ||
onRightBoundary_(globalPos))
@@ -209,17 +207,16 @@ public:
* \brief Evaluate the boundary conditions for a dirichlet
* boundary segment.
*
+ * \param values The dirichlet values for the primary variables
+ * \param vertex The vertex for which the boundary type is set
+ *
* 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
{
+ const GlobalPosition globalPos = vertex.geometry().center();
// use initial values as boundary conditions
- initial(values, element, fvElemGeom, scvIdx);
+ initial_(values, globalPos);
}
/*!
@@ -280,20 +277,23 @@ public:
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
- //const GlobalPosition &pos = element.geometry().corner(scvIdx);
+ const GlobalPosition pos = element.geometry().corner(scvIdx);
+ initial_(values, pos);
+ };
+
+ // \}
+
+private:
+ void initial_(PrimaryVariables &values, const GlobalPosition &pos) const
+ {
Scalar Sw = 0.0;
Scalar pc =
- MaterialLaw::pC(this->spatialParameters().materialLawParams(element,
- fvElemGeom,
- scvIdx),
+ MaterialLaw::pC(this->spatialParameters().materialLawParams(pos),
Sw);
values[pwIdx] = pNreference() - pc;
}
- // \}
-
-private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[0] < this->bboxMin()[0] + eps_;
diff --git a/test/boxmodels/richards/richardslensspatialparameters.hh b/test/boxmodels/richards/richardslensspatialparameters.hh
index 9a4c7cbf5c..d8feb7c014 100644
--- a/test/boxmodels/richards/richardslensspatialparameters.hh
+++ b/test/boxmodels/richards/richardslensspatialparameters.hh
@@ -119,7 +119,12 @@ public:
const FVElementGeometry &fvElemGeom,
int scvIdx) const
{
- const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ return materialLawParams(fvElemGeom.subContVol[scvIdx].global);
+ }
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const GlobalPosition &globalPos) const
+ {
if (isInLens_(globalPos))
return lensMaterialParams_;
diff --git a/tutorial/tutorialproblem_coupled.hh b/tutorial/tutorialproblem_coupled.hh
index 4cec37ffe4..f7b21113eb 100644
--- a/tutorial/tutorialproblem_coupled.hh
+++ b/tutorial/tutorialproblem_coupled.hh
@@ -127,14 +127,9 @@ public:
// Specifies which kind of boundary condition should be used for
// which equation on a given boundary segment.
- void boundaryTypes(BoundaryTypes &BCtype,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
- const Intersection &isIt,
- int scvIdx,
- int boundaryFaceIdx) const
+ void boundaryTypes(BoundaryTypes &BCtype, const Vertex &vertex) const
{
- const GlobalPosition &pos = element.geometry().corner(scvIdx);
+ const GlobalPosition &pos = vertex.geometry().center();
if (pos[0] < eps_) // dirichlet conditions on left boundary
BCtype.setAllDirichlet();
else // neuman for the remaining boundaries
@@ -145,12 +140,7 @@ public:
// Evaluate the boundary conditions for a dirichlet boundary
// segment. For this method, the 'values' parameter stores
// primary variables.
- void dirichlet(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvElemGeom,
- const Intersection &isIt,
- int scvIdx,
- int boundaryFaceIdx) const
+ void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
{
values[Indices::pwIdx] = 200.0e3; // 200 kPa = 2 bar
values[Indices::SnIdx] = 0.0; // 0 % oil saturation on left boundary
--
GitLab