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

extend/correct a few comments, remove artifacts of renamed concepts

(phaseState -> fluidState, material context -> material parameters, etc)

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4834 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent f2370a73
......@@ -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
......
......@@ -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_;
......
......@@ -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_;
}
......
......@@ -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
......
......@@ -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
......
......@@ -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) {
......
......@@ -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
......
......@@ -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:
......
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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_;
......
......@@ -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_;
......
......@@ -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_;
......
......@@ -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_;
......
......@@ -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_;
}
......
......@@ -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_;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment