diff --git a/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh b/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
index 50bb1459f8e52cd411918fdadc27ed25b28192fc..c0cab07b07314876cfd2d32665b2199d3804e955 100644
--- a/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
+++ b/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
@@ -144,7 +144,7 @@ public:
return outerPorosity_;
}
- // return the brooks-corey context depending on the position
+ // return the parameter object for the Brooks-Corey material law which depends on the position
const MaterialLawParams& materialLawParams(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh b/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
index 150d22ece9be2fcdfa13cef2a49e898ee8ebe499..7a0d85f01c6e21011def2e5ca3135e1d93c9ff40 100644
--- a/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
@@ -65,7 +65,7 @@ public:
}
- // return the brooks-corey context depending on the position
+ // return the parameter object for the Brooks-Corey material law which depends on the position
const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
{
return materialLawParams_;
diff --git a/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh b/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
index 868caa590d7666e83e9de80344b9bd188f6e6bad..1e6be79da4bde0a22de66226dfac6dd8a824ea0c 100644
--- a/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
+++ b/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
@@ -65,10 +65,10 @@ public:
}
- // return the brooks-corey context depending on the position
+ // return the brooks-corey material parameter object which depends on the position
const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
{
- return materialLawParams_;
+ return materialLawParams_;
}
diff --git a/dumux/boxmodels/2p/2plocalresidual.hh b/dumux/boxmodels/2p/2plocalresidual.hh
index ac135ff374ec2b437383b977375d71c651ee9bfc..953a657512c340580340491f5d58da7f0261ebce 100644
--- a/dumux/boxmodels/2p/2plocalresidual.hh
+++ b/dumux/boxmodels/2p/2plocalresidual.hh
@@ -152,9 +152,13 @@ public:
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
// calculate the flux in the normal direction of the
- // current sub control volume face
- fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(phaseIdx),
- tmpVec);
+ // current sub control volume face:
+ //
+ // v = - (K grad p) * n
+ //
+ // (the minus comes from the Darcy law which states that
+ // the flux is from high to low pressure potentials.)
+ fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(phaseIdx), tmpVec);
Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
// data attached to upstream and the downstream vertices
diff --git a/dumux/boxmodels/2p2c/2p2cproperties.hh b/dumux/boxmodels/2p2c/2p2cproperties.hh
index dfd1c9e18316cc3b87aaa56f0ea7711b0dbb93ed..9ebcfa6fed96f85e42f67780ab4cc81fa8513a0c 100644
--- a/dumux/boxmodels/2p2c/2p2cproperties.hh
+++ b/dumux/boxmodels/2p2c/2p2cproperties.hh
@@ -53,7 +53,7 @@ NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters
NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
-NEW_PROP_TAG(MaterialLawParams); //!< The context material law (extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters)
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
NEW_PROP_TAG(MobilityUpwindAlpha); //!< The value of the upwind parameter for the mobility
diff --git a/dumux/boxmodels/common/pdelabboxassembler.hh b/dumux/boxmodels/common/pdelabboxassembler.hh
index 2f3cce1dbf30966e5b9f089762a96021ab723fc2..fbac1a198535cf0852463a4f31317611ae6f4178 100644
--- a/dumux/boxmodels/common/pdelabboxassembler.hh
+++ b/dumux/boxmodels/common/pdelabboxassembler.hh
@@ -337,9 +337,9 @@ public:
*
* \param relTol The relative error below which a vertex won't be
* reassembled. Note that this specifies the
- * relative error between the current evaluation
- * point and the current solution and _not_ the
- * delta vector of the Newton iteration!
+ * worst-case relative error between the last
+ * linearization point and the current solution and
+ * _not_ the delta vector of the Newton iteration!
*/
void computeColors(Scalar relTol)
{
@@ -349,10 +349,9 @@ public:
ElementIterator elemIt = gridView_().template begin<0>();
ElementIterator elemEndIt = gridView_().template end<0>();
- greenElems_ = 0;
+ // mark the red vertices and update the tolerance of the
+ // linearization which actually will get achieved
nextReassembleTolerance_ = 0;
-
- // mark the red vertices
for (int i = 0; i < vertexColor_.size(); ++i) {
vertexColor_[i] = Green;
if (vertexDelta_[i] > relTol) {
@@ -366,26 +365,33 @@ public:
// Mark all red elements
for (; elemIt != elemEndIt; ++elemIt) {
- bool needReassemble = false;
+ // find out whether the current element features a red
+ // vertex
+ bool isRed = false;
int numVerts = elemIt->template count<dim>();
for (int i=0; i < numVerts; ++i) {
int globalI = vertexMapper_().map(*elemIt, i, dim);
if (vertexColor_[globalI] == Red) {
- needReassemble = true;
+ isRed = true;
break;
}
};
+ // if yes, the element color is also red, else it is not
+ // red, i.e. green for the mean time
int globalElemIdx = elementMapper_().map(*elemIt);
- elementColor_[globalElemIdx] = needReassemble?Red:Green;
+ if (isRed)
+ elementColor_[globalElemIdx] = Red;
+ else
+ elementColor_[globalElemIdx] = Green;
}
- // Mark all yellow vertices
+ // Mark yellow vertices (as orange for the mean time)
elemIt = gridView_().template begin<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = this->elementMapper_().map(*elemIt);
- if (elementColor_[elemIdx] == Green)
- continue; // green elements do not tint vertices
+ if (elementColor_[elemIdx] != Red)
+ continue; // non-red elements do not tint vertices
// yellow!
int numVerts = elemIt->template count<dim>();
@@ -398,8 +404,7 @@ public:
};
}
- // Mark all yellow elements
- int numGreen = 0;
+ // Mark yellow elements
elemIt = gridView_().template begin<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = this->elementMapper_().map(*elemIt);
@@ -407,6 +412,8 @@ public:
continue; // element is red already!
}
+ // check whether the element features a yellow
+ // (resp. orange at this point) vertex
bool isYellow = false;
int numVerts = elemIt->template count<dim>();
for (int i=0; i < numVerts; ++i) {
@@ -417,15 +424,14 @@ public:
}
};
- if (isYellow) {
+ if (isYellow)
elementColor_[elemIdx] = Yellow;
- }
- else // elementColor_[elemIdx] == Green
- ++ numGreen;
}
// Demote orange vertices to yellow ones if it has at least
- // one green element as a neighbor
+ // one green element as a neighbor. also count the total
+ // number of green elements for statistical purposes.
+ greenElems_ = 0;
elemIt = gridView_().template begin<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = this->elementMapper_().map(*elemIt);
@@ -444,17 +450,19 @@ public:
};
}
- // promote the remaining orange vertices to red ones
+ // promote the remaining orange vertices to red
for (int i=0; i < vertexColor_.size(); ++i) {
- // if a vertex is green or yellow don't recolor it!
+ // if a vertex is green or yellow don't do anything!
if (vertexColor_[i] == Green || vertexColor_[i] == Yellow)
continue;
+
+ // make sure the vertex is red (this is a no-op vertices
+ // which are already red!)
+ vertexColor_[i] = Red;
- // set discrepancy for this vertex to 0 because the
- // system will be consistently linearized for this
- // vertex
+ // set the error of this vertex to 0 because the system
+ // will be consistently linearized at this vertex
vertexDelta_[i] = 0.0;
- vertexColor_[i] = Red;
};
greenElems_ = gridView_().comm().sum(greenElems_);
@@ -555,9 +563,10 @@ private:
// allocate raw matrix
matrix_ = new Matrix(nVerts, nVerts, Matrix::random);
- // find out how many neighbors each vertex has
+ // find out the global indices of the neighboring vertices of
+ // each vertex
typedef std::set<int> NeighborSet;
- std::vector<NeighborSet > neighbors(nVerts);
+ std::vector<NeighborSet> neighbors(nVerts);
ElementIterator eIt = gridView_().template begin<0>();
const ElementIterator eEndIt = gridView_().template end<0>();
for (; eIt != eEndIt; ++eIt) {
@@ -565,25 +574,32 @@ private:
// loop over all element vertices
int n = elem.template count<dim>();
- for (int i = 0; i < n; ++i) {
+ for (int i = 0; i < n - 1; ++i) {
int globalI = vertexMapper_().map(*eIt, i, dim);
- for (int j = 0; j < n; ++j) {
+ for (int j = i + 1; j < n; ++j) {
int globalJ = vertexMapper_().map(*eIt, j, dim);
- // insert into BCRS matrix
+ // make sure that vertex j is in the neighbor set
+ // of vertex i and vice-versa
neighbors[globalI].insert(globalJ);
+ neighbors[globalJ].insert(globalI);
}
}
};
- // allocate space for the rows
+ // make vertices neighbors to themselfs
+ for (int i = 0; i < nVerts; ++i)
+ neighbors[i].insert(i);
+
+ // allocate space for the rows of the matrix
for (int i = 0; i < nVerts; ++i) {
matrix_->setrowsize(i, neighbors[i].size());
}
matrix_->endrowsizes();
- // allocate space for the rows
+ // fill the rows with indices. each vertex talks to all of its
+ // neighbors. (it also talks to itself since vertices are
+ // sometimes quite egocentric.)
for (int i = 0; i < nVerts; ++i) {
- // off-diagonal entries
typename NeighborSet::iterator nIt = neighbors[i].begin();
typename NeighborSet::iterator nEndIt = neighbors[i].end();
for (; nIt != nEndIt; ++nIt) {
diff --git a/dumux/boxmodels/richards/richardsproperties.hh b/dumux/boxmodels/richards/richardsproperties.hh
index 52f64170aa0d8e9cb103b5ac1ccc123d0c30009e..b76b4a48a3a3a03f8ce6f781d22b72127c8a0cf7 100644
--- a/dumux/boxmodels/richards/richardsproperties.hh
+++ b/dumux/boxmodels/richards/richardsproperties.hh
@@ -50,7 +50,7 @@ NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
NEW_PROP_TAG(RichardsIndices); //!< Enumerations used by the Richards models
NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters object
NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (by default extracted from the spatial parameters)
-NEW_PROP_TAG(MaterialLawParams); //!< The context material law (by default extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The type of the parameter object for the material law (by default extracted from the spatial parameters)
NEW_PROP_TAG(FluidSystem); //!< The fluid system to be used for the Richards model
NEW_PROP_TAG(FluidState); //!< The fluid state to be used for the Richards model
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
diff --git a/dumux/material/fluidsystems/2p_system.hh b/dumux/material/fluidsystems/2p_system.hh
index ad65252308b3d26514125e5a522e9515ade8905f..583f121575e735b9c9b2d962b3747747bc9ff834 100644
--- a/dumux/material/fluidsystems/2p_system.hh
+++ b/dumux/material/fluidsystems/2p_system.hh
@@ -135,7 +135,7 @@ public:
* \param phaseIdx index of the phase
* \param temperature phase temperature in [K]
* \param pressure phase pressure in [Pa]
- * \param phaseState The fluid state of the two-phase model
+ * \param fluidState The fluid state of the two-phase model
* \tparam FluidState the fluid state class of the two-phase model
* \return returns the density of the phase
*/
@@ -143,7 +143,7 @@ public:
static Scalar phaseDensity(int phaseIdx,
Scalar temperature,
Scalar pressure,
- const FluidState &phaseState)
+ const FluidState &fluidState)
{
switch (phaseIdx) {
case wPhaseIdx:
@@ -160,7 +160,7 @@ public:
* \param phaseIdx index of the phase
* \param temperature phase temperature in [K]
* \param pressure phase pressure in [Pa]
- * \param phaseState The fluid state of the two-phase model
+ * \param fluidState The fluid state of the two-phase model
* \tparam FluidState the fluid state class of the two-phase model
* \return returns the viscosity of the phase [Pa*s]
*/
@@ -168,7 +168,7 @@ public:
static Scalar phaseViscosity(int phaseIdx,
Scalar temperature,
Scalar pressure,
- const FluidState &phaseState)
+ const FluidState &fluidState)
{
switch (phaseIdx) {
case wPhaseIdx:
@@ -185,15 +185,15 @@ public:
* \param phaseIdx index of the phase
* \param temperature phase temperature in [K]
* \param pressure phase pressure in [Pa]
- * \param phaseState The fluid state of the two-phase model
+ * \param fluidState The fluid state of the two-phase model
* \tparam FluidState the fluid state class of the two-phase model
* \return returns the specific enthalpy of the phase [J/kg]
*/
template <class FluidState>
static Scalar phaseEnthalpy(int phaseIdx,
- Scalar temperature,
- Scalar pressure,
- const FluidState &phaseState)
+ Scalar temperature,
+ Scalar pressure,
+ const FluidState &fluidState)
{
switch (phaseIdx) {
case wPhaseIdx:
@@ -210,15 +210,15 @@ public:
* \param phaseIdx index of the phase
* \param temperature phase temperature in [K]
* \param pressure phase pressure in [Pa]
- * \param phaseState The fluid state of the two-phase model
+ * \param fluidState The fluid state of the two-phase model
* \tparam FluidState the fluid state class of the two-phase model
* \return returns the specific internal energy of the phase [J/kg]
*/
template <class FluidState>
static Scalar phaseInternalEnergy(int phaseIdx,
- Scalar temperature,
- Scalar pressure,
- const FluidState &phaseState)
+ Scalar temperature,
+ Scalar pressure,
+ const FluidState &fluidState)
{
switch (phaseIdx) {
case wPhaseIdx:
diff --git a/test/boxmodels/2p/lensspatialparameters.hh b/test/boxmodels/2p/lensspatialparameters.hh
index e74ad8493b6af23fa476eca3dbf9c58376941a96..7812b97ea451afa0773c4115848fee9bc0e2c821 100644
--- a/test/boxmodels/2p/lensspatialparameters.hh
+++ b/test/boxmodels/2p/lensspatialparameters.hh
@@ -125,7 +125,7 @@ public:
int scvIdx) const
{ return 0.4; }
- // return the brooks-corey context depending on the position
+ // return the parameter object for the Brooks-Corey material law which depends on the position
const MaterialLawParams& materialLawParams(const Element &element,
const FVElementGeometry &fvElemGeom,
int scvIdx) const
diff --git a/test/boxmodels/2p2c/injectionspatialparameters.hh b/test/boxmodels/2p2c/injectionspatialparameters.hh
index 392df67e5a97f3853ba0ab560593c6ea5d289c35..38ecf667f55d1efda2599370b0b6c61da850b70b 100644
--- a/test/boxmodels/2p2c/injectionspatialparameters.hh
+++ b/test/boxmodels/2p2c/injectionspatialparameters.hh
@@ -160,7 +160,7 @@ public:
/*!
- * \brief return the brooks-corey context depending on the position
+ * \brief return the parameter object for the Brooks-Corey material law which depends on the position
*
* \param element The current finite element
* \param fvElemGeom The current finite volume geometry of the element
diff --git a/test/boxmodels/2p2cni/waterairspatialparameters.hh b/test/boxmodels/2p2cni/waterairspatialparameters.hh
index 2e9dcde0a4cf1b4b4f63012e916500de7a6fd9ff..a3d41703132544ad1d11b7e2a1d859be4b9eeb65 100644
--- a/test/boxmodels/2p2cni/waterairspatialparameters.hh
+++ b/test/boxmodels/2p2cni/waterairspatialparameters.hh
@@ -160,7 +160,7 @@ public:
/*!
- * \brief return the brooks-corey context depending on the position
+ * \brief return the parameter object for the Brooks-Corey material law which depends on the position
*
* \param element The current finite element
* \param fvElemGeom The current finite volume geometry of the element
diff --git a/test/decoupled/1p/test_diffusion_spatialparams.hh b/test/decoupled/1p/test_diffusion_spatialparams.hh
index fd392174413e82e9c8c827786866678720d14a13..b18507c83792e7350f6deb9c11fa20d09daded43 100644
--- a/test/decoupled/1p/test_diffusion_spatialparams.hh
+++ b/test/decoupled/1p/test_diffusion_spatialparams.hh
@@ -73,7 +73,7 @@ public:
}
- // return the brooks-corey context depending on the position
+ // return the parameter object for the Brooks-Corey material law which depends on the position
const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
{
return materialLawParams_;
diff --git a/test/decoupled/2p/test_impes_spatialparams.hh b/test/decoupled/2p/test_impes_spatialparams.hh
index cab3c4125bec45439ee67f3e70a6dfe9a8113e7d..9049f81a582613014e34b7e127eb558e350ef63a 100644
--- a/test/decoupled/2p/test_impes_spatialparams.hh
+++ b/test/decoupled/2p/test_impes_spatialparams.hh
@@ -70,7 +70,7 @@ public:
}
- // return the brooks-corey context depending on the position
+ // return the parameter object for the Brooks-Corey material law which depends on the position
const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
{
return materialLawParams_;
diff --git a/test/decoupled/2p/test_transport_spatialparams.hh b/test/decoupled/2p/test_transport_spatialparams.hh
index 14e3fd6d0e278b8ec66b99f4b9b54bba317fcc28..8c17642e688a2b9a35442523ea21a00cdd5813f2 100644
--- a/test/decoupled/2p/test_transport_spatialparams.hh
+++ b/test/decoupled/2p/test_transport_spatialparams.hh
@@ -68,7 +68,7 @@ public:
}
- // return the brooks-corey context depending on the position
+ // return the parameter object for the Brooks-Corey material law which depends on the position
const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
{
return materialLawParams_;
diff --git a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
index 4c78fccda37e041f64a8fbf8ef454643b9bf3a1d..088e31c295b4f36a548d04f0bfaee82463f0f9da 100644
--- a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
+++ b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
@@ -70,7 +70,7 @@ public:
}
- // return the brooks-corey context depending on the position
+ // return the parameter object for the Brooks-Corey material law which depends on the position
const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
{
return materialLawParams_;
diff --git a/tutorial/tutorialspatialparameters_coupled.hh b/tutorial/tutorialspatialparameters_coupled.hh
index b6f62f29aa7bbeec48d6a0d2c0eb4decb89ac945..f595b6c5bd0b6b747d605b75f470b27b6b2de199 100644
--- a/tutorial/tutorialspatialparameters_coupled.hh
+++ b/tutorial/tutorialspatialparameters_coupled.hh
@@ -56,11 +56,11 @@ class TutorialSpatialParametersCoupled: public BoxSpatialParameters<TypeTag> /*@
typedef typename Grid::Traits::template Codim<0>::Entity Element;
// select material law to be used
- typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw; /*@\label{tutorial-coupled:rawlaw}@*/
+ typedef RegularizedBrooksCorey<Scalar> EffectiveMaterialLaw; /*@\label{tutorial-coupled:rawlaw}@*/
public:
// adapter for absolute law
- typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw; /*@\label{tutorial-coupled:eff2abs}@*/
+ typedef EffToAbsLaw<EffectiveMaterialLaw> MaterialLaw; /*@\label{tutorial-coupled:eff2abs}@*/
// determine appropriate parameters depening on selected materialLaw
typedef typename MaterialLaw::Params MaterialLawParams; /*@\label{tutorial-coupled:matLawObjectType}@*/
@@ -83,10 +83,11 @@ public:
return 0.2;
}
- // return the materialLaw context (i.e. BC, regularizedVG, etc) depending on the position
+ // return the parameter object for the material law (i.e. Brooks-Corey)
+ // which may vary with the spatial position
const MaterialLawParams& materialLawParams(const Element &element, /*@\label{tutorial-coupled:matLawParams}@*/
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
{
return materialParams_;
}
diff --git a/tutorial/tutorialspatialparameters_decoupled.hh b/tutorial/tutorialspatialparameters_decoupled.hh
index 22479f852bd639e300d846b2ae0679e2f5930872..bc021ae22ed5ab6cb290f12f15170e2bdf6051cc 100644
--- a/tutorial/tutorialspatialparameters_decoupled.hh
+++ b/tutorial/tutorialspatialparameters_decoupled.hh
@@ -47,10 +47,10 @@ class TutorialSpatialParametersDecoupled
typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
// material law typedefs
- typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw;
-// typedef LinearMaterial<Scalar> RawMaterialLaw;
+ typedef RegularizedBrooksCorey<Scalar> EffectiveMaterialLaw;
+// typedef LinearMaterial<Scalar> EffectiveMaterialLaw;
public:
- typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw;
+ typedef EffToAbsLaw<EffectiveMaterialLaw> MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
//! Update the spatial parameters with the flow solution after a timestep.
@@ -73,7 +73,8 @@ public:
return 0.2;
}
- //! return the material law context (i.e. BC, regularizedVG, etc) depending on the position
+ //! return the parameter object for the material law (i.e. Brooks-Corey)
+ //! which may vary with the spatial position
const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
{
return materialLawParams_;