Commit ace814e6 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

incorporated naming conventions in boxmodels/common

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8116 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent fdb398a5
......@@ -81,11 +81,9 @@ public:
* \param problem The problem object which needs to be simulated
* \param element The DUNE Codim<0> entity for which the boundary
* types should be collected
* \param fvElemGeom The element's finite volume geometry
*/
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvElemGeom)
const Element &element)
{
int numVerts = element.template count<dim>();
this->resize(numVerts);
......@@ -109,6 +107,20 @@ public:
}
}
/*!
* \brief Update the boundary types for all vertices of an element.
*
* \param problem The problem object which needs to be simulated
* \param element The DUNE Codim<0> entity for which the boundary
* types should be collected
* \param fvGeometry The element's finite volume geometry
*/
DUMUX_DEPRECATED_MSG("use update(problem, element) instead")
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry)
{ update(problem, element); }
/*!
* \brief Returns whether the element has a vertex which contains
* a Dirichlet value.
......
......@@ -66,12 +66,12 @@ public:
*
* \param problem The problem which needs to be simulated.
* \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
* \param fvElemGeom The finite volume geometry of the element
* \param fvGeometry The finite volume geometry of the element
* \param oldSol Tells whether the model's previous or current solution should be used.
*/
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvElemGeom,
const FVElementGeometry &fvGeometry,
bool oldSol)
{
const SolutionVector &globalSol =
......@@ -81,16 +81,16 @@ public:
const VertexMapper &vertexMapper = problem.vertexMapper();
// we assert that the i-th shape function is
// associated to the i-th vert of the element.
int n = element.template count<dim>();
this->resize(n);
for (int i = 0; i < n; i++) {
int numVertices = element.template count<dim>();
this->resize(numVertices);
for (int scvIdx = 0; scvIdx < numVertices; scvIdx++) {
const PrimaryVariables &solI
= globalSol[vertexMapper.map(element, i, dim)];
(*this)[i].update(solI,
= globalSol[vertexMapper.map(element, scvIdx, dim)];
(*this)[scvIdx].update(solI,
problem,
element,
fvElemGeom,
i,
fvGeometry,
scvIdx,
oldSol);
}
}
......@@ -99,33 +99,33 @@ public:
* \brief Construct the volume variables for all of vertices of an
* element given a solution vector computed by PDELab.
*
* \tparam ElemSolVectorType The container type which stores the
* \tparam ElementSolutionVector The container type which stores the
* primary variables of the element
* using _local_ indices
*
* \param problem The problem which needs to be simulated.
* \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
* \param fvElemGeom The finite volume geometry of the element
* \param fvGeometry The finite volume geometry of the element
* \param elementSolVector The local solution for the element using PDELab ordering
*/
template<typename ElemSolVectorType>
template<typename ElementSolutionVector>
void updatePDELab(const Problem &problem,
const Element &element,
const FVElementGeometry &fvElemGeom,
const ElemSolVectorType& elementSolVector)
const FVElementGeometry &fvGeometry,
const ElementSolutionVector& elementSolVector)
{
int n = element.template count<dim>();
this->resize(n);
for (int vertexIdx = 0; vertexIdx < n; vertexIdx++)
int numVertices = element.template count<dim>();
this->resize(numVertices);
for (int scvIdx = 0; scvIdx < numVertices; scvIdx++)
{
PrimaryVariables solI(0);
for (int eqnIdx=0; eqnIdx<numEq; eqnIdx++)
solI[eqnIdx] = elementSolVector[vertexIdx + eqnIdx*n];
(*this)[vertexIdx].update(solI,
for (int eqnIdx = 0; eqnIdx < numEq; eqnIdx++)
solI[eqnIdx] = elementSolVector[scvIdx + eqnIdx*numVertices];
(*this)[scvIdx].update(solI,
problem,
element,
fvElemGeom,
vertexIdx,
fvGeometry,
scvIdx,
false);
}
......
......@@ -141,14 +141,14 @@ public:
fvElemGeom_.update(gridView_(), element);
reset_();
bcTypes_.update(problem_(), elem_(), fvElemGeom_);
bcTypes_.update(problem_(), element_(), fvElemGeom_);
// this is pretty much a HACK because the internal state of
// the problem is not supposed to be changed during the
// evaluation of the residual. (Reasons: It is a violation of
// abstraction, makes everything more prone to errors and is
// not thread save.) The real solution are context objects!
problem_().updateCouplingParams(elem_());
problem_().updateCouplingParams(element_());
// set the hints for the volume variables
model_().setHints(element, prevVolVars_, curVolVars_);
......@@ -156,12 +156,12 @@ public:
// update the secondary variables for the element at the last
// and the current time levels
prevVolVars_.update(problem_(),
elem_(),
element_(),
fvElemGeom_,
true /* isOldSol? */);
curVolVars_.update(problem_(),
elem_(),
element_(),
fvElemGeom_,
false /* isOldSol? */);
......@@ -169,7 +169,7 @@ public:
model_().updateCurHints(element, curVolVars_);
// calculate the local residual
localResidual().eval(elem_(),
localResidual().eval(element_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
......@@ -177,7 +177,7 @@ public:
residual_ = localResidual().residual();
storageTerm_ = localResidual().storageTerm();
model_().updatePVWeights(elem_(), curVolVars_);
model_().updatePVWeights(element_(), curVolVars_);
// calculate the local jacobian matrix
int numVertices = fvElemGeom_.numVertices;
......@@ -303,6 +303,16 @@ protected:
/*!
* \brief Returns a reference to the element.
*/
const Element &element_() const
{
Valgrind::CheckDefined(elemPtr_);
return *elemPtr_;
};
/*!
* \brief Returns a reference to the element.
*/
DUMUX_DEPRECATED_MSG("use element_() instead")
const Element &elem_() const
{
Valgrind::CheckDefined(elemPtr_);
......@@ -332,7 +342,7 @@ protected:
*/
void reset_()
{
int n = elem_().template count<dim>();
int n = element_().template count<dim>();
for (int i = 0; i < n; ++ i) {
storageJacobian_[i] = 0.0;
for (int j = 0; j < n; ++ j) {
......@@ -375,9 +385,9 @@ protected:
* is the value of a sub-control volume's primary variable at the
* evaluation point and \f$\epsilon\f$ is a small value larger than 0.
*
* \param dest The vector storing the partial derivatives of all
* \param partialDeriv The vector storing the partial derivatives of all
* equations
* \param destStorage the mass matrix contributions
* \param storageDeriv the mass matrix contributions
* \param scvIdx The sub-control volume index of the current
* finite element for which the partial derivative
* ought to be calculated
......@@ -386,12 +396,12 @@ protected:
* for which the partial derivative ought to be
* calculated
*/
void evalPartialDerivative_(ElementSolutionVector &dest,
PrimaryVariables &destStorage,
void evalPartialDerivative_(ElementSolutionVector &partialDeriv,
PrimaryVariables &storageDeriv,
int scvIdx,
int pvIdx)
{
int globalIdx = vertexMapper_().map(elem_(), scvIdx, dim);
int globalIdx = vertexMapper_().map(element_(), scvIdx, dim);
PrimaryVariables priVars(model_().curSol()[globalIdx]);
VolumeVariables origVolVars(curVolVars_[scvIdx]);
......@@ -411,26 +421,26 @@ protected:
// calculate the residual
curVolVars_[scvIdx].update(priVars,
problem_(),
elem_(),
element_(),
fvElemGeom_,
scvIdx,
false);
localResidual().eval(elem_(),
localResidual().eval(element_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
// store the residual and the storage term
dest = localResidual().residual();
destStorage = localResidual().storageTerm()[scvIdx];
partialDeriv = localResidual().residual();
storageDeriv = localResidual().storageTerm()[scvIdx];
}
else {
// we are using backward differences, i.e. we don't need
// to calculate f(x + \epsilon) and we can recycle the
// (already calculated) residual f(x)
dest = residual_;
destStorage = storageTerm_[scvIdx];
partialDeriv = residual_;
storageDeriv = storageTerm_[scvIdx];
}
......@@ -445,37 +455,37 @@ protected:
// calculate residual again
curVolVars_[scvIdx].update(priVars,
problem_(),
elem_(),
element_(),
fvElemGeom_,
scvIdx,
false);
localResidual().eval(elem_(),
localResidual().eval(element_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
dest -= localResidual().residual();
destStorage -= localResidual().storageTerm()[scvIdx];
partialDeriv -= localResidual().residual();
storageDeriv -= localResidual().storageTerm()[scvIdx];
}
else {
// we are using forward differences, i.e. we don't need to
// calculate f(x - \epsilon) and we can recycle the
// (already calculated) residual f(x)
dest -= residual_;
destStorage -= storageTerm_[scvIdx];
partialDeriv -= residual_;
storageDeriv -= storageTerm_[scvIdx];
}
// divide difference in residuals by the magnitude of the
// deflections between the two function evaluation
dest /= delta;
destStorage /= delta;
partialDeriv /= delta;
storageDeriv /= delta;
// restore the original state of the element's volume variables
curVolVars_[scvIdx] = origVolVars;
#if HAVE_VALGRIND
for (unsigned i = 0; i < dest.size(); ++i)
Valgrind::CheckDefined(dest[i]);
for (unsigned i = 0; i < partialDeriv.size(); ++i)
Valgrind::CheckDefined(partialDeriv[i]);
#endif
}
......@@ -486,7 +496,7 @@ protected:
*/
void updateLocalJacobian_(int scvIdx,
int pvIdx,
const ElementSolutionVector &deriv,
const ElementSolutionVector &partialDeriv,
const PrimaryVariables &storageDeriv)
{
// store the derivative of the storage term
......@@ -496,7 +506,7 @@ protected:
for (int i = 0; i < fvElemGeom_.numVertices; i++)
{
if (jacAsm_().vertexColor(elem_(), i) == Green) {
if (jacAsm_().vertexColor(element_(), i) == Green) {
// Green vertices are not to be changed!
continue;
}
......@@ -506,7 +516,7 @@ protected:
// the residual of equation 'eqIdx' at vertex 'i'
// depending on the primary variable 'pvIdx' at vertex
// 'scvIdx'.
this->A_[i][scvIdx][eqIdx][pvIdx] = deriv[i][eqIdx];
this->A_[i][scvIdx][eqIdx][pvIdx] = partialDeriv[i][eqIdx];
Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
}
}
......
......@@ -92,11 +92,11 @@ public:
* This assumes that all objects of the simulation have been fully
* allocated but not necessarily initialized completely.
*
* \param prob The representation of the physical problem to be
* \param problem The representation of the physical problem to be
* solved.
*/
void init(Problem &prob)
{ problemPtr_ = &prob; }
void init(Problem &problem)
{ problemPtr_ = &problem; }
/*!
* \brief Compute the local residual, i.e. the deviation of the
......@@ -107,10 +107,10 @@ public:
*/
void eval(const Element &element)
{
FVElementGeometry fvElemGeom;
FVElementGeometry fvGeometry;
fvElemGeom.update(gridView_(), element);
fvElemGeomPtr_ = &fvElemGeom;
fvGeometry.update(gridView_(), element);
fvElemGeomPtr_ = &fvGeometry;
ElementVolumeVariables volVarsPrev, volVarsCur;
// update the hints
......@@ -151,9 +151,9 @@ public:
{
elemPtr_ = &element;
FVElementGeometry fvElemGeom;
fvElemGeom.update(gridView_(), element);
fvElemGeomPtr_ = &fvElemGeom;
FVElementGeometry fvGeometry;
fvGeometry.update(gridView_(), element);
fvElemGeomPtr_ = &fvGeometry;
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
......@@ -187,9 +187,9 @@ public:
{
elemPtr_ = &element;
FVElementGeometry fvElemGeom;
fvElemGeom.update(gridView_(), element);
fvElemGeomPtr_ = &fvElemGeom;
FVElementGeometry fvGeometry;
fvGeometry.update(gridView_(), element);
fvElemGeomPtr_ = &fvGeometry;
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
......@@ -209,7 +209,7 @@ public:
*
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
* \param fvElemGeom The finite-volume geometry of the element
* \param fvGeometry The finite-volume geometry of the element
* \param prevVolVars The volume averaged variables for all
* sub-control volumes of the element at the previous
* time level
......@@ -220,7 +220,7 @@ public:
* vertices of the element
*/
void eval(const Element &element,
const FVElementGeometry &fvElemGeom,
const FVElementGeometry &fvGeometry,
const ElementVolumeVariables &prevVolVars,
const ElementVolumeVariables &curVolVars,
const ElementBoundaryTypes &bcTypes)
......@@ -229,14 +229,14 @@ public:
Valgrind::CheckDefined(curVolVars);
#if !defined NDEBUG && HAVE_VALGRIND
for (int i=0; i < fvElemGeom.numVertices; i++) {
for (int i=0; i < fvGeometry.numVertices; i++) {
prevVolVars[i].checkDefined();
curVolVars[i].checkDefined();
}
#endif // HAVE_VALGRIND
elemPtr_ = &element;
fvElemGeomPtr_ = &fvElemGeom;
fvElemGeomPtr_ = &fvGeometry;
bcTypesPtr_ = &bcTypes;
prevVolVarsPtr_ = &prevVolVars;
curVolVarsPtr_ = &curVolVars;
......@@ -353,16 +353,16 @@ protected:
*/
void evalDirichlet_()
{
PrimaryVariables tmp(0);
for (int i = 0; i < fvGeometry_().numVertices; ++i) {
const BoundaryTypes &bcTypes = bcTypes_(i);
PrimaryVariables dirichletValues(0);
for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; ++scvIdx) {
const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
if (! bcTypes.hasDirichlet())
continue;
// ask the problem for the dirichlet values
const VertexPointer vPtr = element_().template subEntity<dim>(i);
Valgrind::SetUndefined(tmp);
asImp_().problem_().dirichlet(tmp, *vPtr);
const VertexPointer vPtr = element_().template subEntity<dim>(scvIdx);
Valgrind::SetUndefined(dirichletValues);
asImp_().problem_().dirichlet(dirichletValues, *vPtr);
// set the dirichlet conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
......@@ -370,12 +370,12 @@ protected:
continue;
int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
assert(0 <= pvIdx && pvIdx < numEq);
Valgrind::CheckDefined(tmp[pvIdx]);
Valgrind::CheckDefined(dirichletValues[pvIdx]);
residual_[i][eqIdx] =
curPrimaryVar_(i, pvIdx) - tmp[pvIdx];
residual_[scvIdx][eqIdx] =
curPrimaryVar_(scvIdx, pvIdx) - dirichletValues[pvIdx];
storageTerm_[i][eqIdx] = 0.0;
storageTerm_[scvIdx][eqIdx] = 0.0;
};
};
}
......@@ -405,7 +405,7 @@ protected:
faceVertIdx < numFaceVerts;
++faceVertIdx)
{
int elemVertIdx = refElem.subEntity(faceIdx,
int scvIdx = refElem.subEntity(faceIdx,
1,
faceVertIdx,
dim);
......@@ -416,13 +416,13 @@ protected:
// add the residual of all vertices of the boundary
// segment
asImp_().evalNeumannSegment_(isIt,
elemVertIdx,
scvIdx,
boundaryFaceIdx);
// evaluate the outflow conditions at the boundary face
// ATTENTION: This is so far a beta version that is only for the 2p2c and 2p2cni model
// available and not thoroughly tested.
asImp_().evalOutflowSegment(isIt,
elemVertIdx,
scvIdx,
boundaryFaceIdx);
}
}
......@@ -437,29 +437,29 @@ protected:
int boundaryFaceIdx)
{
// temporary vector to store the neumann boundary fluxes
PrimaryVariables values(0.0);
PrimaryVariables neumannFlux(0.0);
const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
// deal with neumann boundaries
if (bcTypes.hasNeumann()) {
Valgrind::SetUndefined(values);
problem_().boxSDNeumann(values,
Valgrind::SetUndefined(neumannFlux);
problem_().boxSDNeumann(neumannFlux,
element_(),
fvGeometry_(),
*isIt,
scvIdx,
boundaryFaceIdx,
curVolVars_());
values *=
neumannFlux *=
fvGeometry_().boundaryFace[boundaryFaceIdx].area
* curVolVars_(scvIdx).extrusionFactor();
Valgrind::CheckDefined(values);
Valgrind::CheckDefined(neumannFlux);
// set the neumann conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!bcTypes.isNeumann(eqIdx))
continue;
residual_[scvIdx][eqIdx] += values[eqIdx];
residual_[scvIdx][eqIdx] += neumannFlux[eqIdx];
}
}
}
......@@ -472,15 +472,15 @@ protected:
{
// calculate the mass flux over the faces and subtract
// it from the local rates
for (int k = 0; k < fvGeometry_().numEdges; k++)
for (int scvfIdx = 0; scvfIdx < fvGeometry_().numEdges; scvfIdx++)
{
int i = fvGeometry_().subContVolFace[k].i;
int j = fvGeometry_().subContVolFace[k].j;
int i = fvGeometry_().subContVolFace[scvfIdx].i;
int j = fvGeometry_().subContVolFace[scvfIdx].j;
PrimaryVariables flux;
Valgrind::SetUndefined(flux);
asImp_().computeFlux(flux, k);
asImp_().computeFlux(flux, scvfIdx);
Valgrind::CheckDefined(flux);
Scalar extrusionFactor =
......@@ -519,13 +519,13 @@ protected:
// calculate the amount of conservation each quantity inside
// all sub control volumes
for (int i=0; i < fvGeometry_().numVertices; i++) {
Valgrind::SetUndefined(storageTerm_[i]);
asImp_().computeStorage(storageTerm_[i], i, /*isOldSol=*/false);
storageTerm_[i] *=
fvGeometry_().subContVol[i].volume
* curVolVars_(i).extrusionFactor();
Valgrind::CheckDefined(storageTerm_[i]);
for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++) {
Valgrind::SetUndefined(storageTerm_[scvIdx]);
asImp_().computeStorage(storageTerm_[scvIdx], scvIdx, /*isOldSol=*/false);
storageTerm_[scvIdx] *=
fvGeometry_().subContVol[scvIdx].volume
* curVolVars_(scvIdx).extrusionFactor();
Valgrind::CheckDefined(storageTerm_[scvIdx]);
}
}
......@@ -537,12 +537,12 @@ protected:
void evalVolumeTerms_()
{
// evaluate the volume terms (storage + source terms)
for (int i=0; i < fvGeometry_().numVertices; i++)
for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++)
{
Scalar extrusionFactor =
curVolVars_(i).extrusionFactor();
curVolVars_(scvIdx).extrusionFactor();
PrimaryVariables tmp(0.);
PrimaryVariables values(0.0);
// mass balance within the element. this is the
// \f$\frac{m}{\partial t}\f$ term if using implicit
......@@ -550,30 +550,30 @@ protected:
//
// TODO (?): we might need a more explicit way for
// doing the time discretization...
Valgrind::SetUndefined(storageTerm_[i]);
Valgrind::SetUndefined(tmp);
asImp_().computeStorage(storageTerm_[i], i, false);
asImp_().computeStorage(tmp, i, true);
Valgrind::CheckDefined(storageTerm_[i]);
Valgrind::CheckDefined(tmp);
storageTerm_[i] -= tmp;
storageTerm_[i] *=
fvGeometry_().subContVol[i].volume
Valgrind::SetUndefined(storageTerm_[scvIdx]);
Valgrind::SetUndefined(values);
asImp_().computeStorage(storageTerm_[scvIdx], scvIdx, false);
asImp_().computeStorage(values, scvIdx, true);
Valgrind::CheckDefined(storageTerm_[scvIdx]);
Valgrind::CheckDefined(values);
storageTerm_[scvIdx] -= values;
storageTerm_[scvIdx] *=
fvGeometry_().subContVol[scvIdx].volume
/ problem_().timeManager().timeStepSize()
* extrusionFactor;
residual_[i] += storageTerm_[i];
residual_[scvIdx] += storageTerm_[scvIdx];
// subtract the source term from the local rate
Valgrind::SetUndefined(tmp);
asImp_().computeSource(tmp, i);