diff --git a/dumux/boxmodels/common/porousmediaboxproblem.hh b/dumux/boxmodels/common/porousmediaboxproblem.hh
index 63cd95cb08977333ad9fdab9814b3b4b3f695f80..c86c92e6d4a1fcd3c0f5c00a362823a30dbabe59 100644
--- a/dumux/boxmodels/common/porousmediaboxproblem.hh
+++ b/dumux/boxmodels/common/porousmediaboxproblem.hh
@@ -1,7 +1,6 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
- * 2012 by Bernd Flemisch *
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
diff --git a/dumux/implicit/box/boxlocalresidual.hh b/dumux/implicit/box/boxlocalresidual.hh
index 8c02754133eec4a5cde1c8afdbff9bbc3d93ece1..3b8fb1a63a8578a11a886719f474c48ca09f53e2 100644
--- a/dumux/implicit/box/boxlocalresidual.hh
+++ b/dumux/implicit/box/boxlocalresidual.hh
@@ -27,6 +27,7 @@
#include <dune/grid/common/geometry.hh>
#include <dumux/common/valgrind.hh>
+#include <dumux/implicit/common/implicitlocalresidual.hh>
#include "boxproperties.hh"
@@ -35,18 +36,16 @@ namespace Dumux
/*!
* \ingroup BoxModel
* \ingroup BoxLocalResidual
- * \brief Element-wise calculation of the residual matrix for models
- * based on the box scheme.
+ * \brief Element-wise calculation of the residual for models
+ * based on the fully implicit box scheme.
*
* \todo Please doc me more!
*/
template<class TypeTag>
-class BoxLocalResidual
+class BoxLocalResidual : public ImplicitLocalResidual<TypeTag>
{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ typedef ImplicitLocalResidual<TypeTag> ParentType;
+ friend class ImplicitLocalResidual<TypeTag>;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -63,273 +62,18 @@ private:
typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
// copying the local residual class is not a good idea
BoxLocalResidual(const BoxLocalResidual &);
public:
- BoxLocalResidual()
+ BoxLocalResidual() : ParentType()
{ }
- ~BoxLocalResidual()
- { }
-
- /*!
- * \brief Initialize the local residual.
- *
- * This assumes that all objects of the simulation have been fully
- * allocated but not necessarily initialized completely.
- *
- * \param problem The representation of the physical problem to be
- * solved.
- */
- void init(Problem &problem)
- { problemPtr_ = &problem; }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- */
- void eval(const Element &element)
- {
- FVElementGeometry fvGeometry;
-
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
-
- ElementVolumeVariables volVarsPrev, volVarsCur;
- // update the hints
- model_().setHints(element, volVarsPrev, volVarsCur);
-
- volVarsPrev.update(problem_(),
- element,
- fvGeometry_(),
- true /* oldSol? */);
- volVarsCur.update(problem_(),
- element,
- fvGeometry_(),
- false /* oldSol? */);
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
-
- // 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(element);
-
- asImp_().eval(element, fvGeometry_(), volVarsPrev, volVarsCur, bcTypes);
- }
-
- /*!
- * \brief Compute the storage term for the current solution.
- *
- * This can be used to figure out how much of each conservation
- * quantity is inside the element.
- *
- * \param element The DUNE Codim<0> entity for which the storage
- * term ought to be calculated
- */
- void evalStorage(const Element &element)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvGeometry;
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
- bcTypesPtr_ = &bcTypes;
-
- // no previous volume variables!
- prevVolVarsPtr_ = 0;
-
- ElementVolumeVariables volVars;
-
- // update the hints
- model_().setHints(element, volVars);
-
- // calculate volume current variables
- volVars.update(problem_(), element, fvGeometry_(), false);
- curVolVarsPtr_ = &volVars;
-
- asImp_().evalStorage_();
- }
-
- /*!
- * \brief Compute the flux term for the current solution.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param curVolVars The volume averaged variables for all
- * sub-contol volumes of the element
- */
- void evalFluxes(const Element &element,
- const ElementVolumeVariables &curVolVars)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvGeometry;
- fvGeometry.update(gridView_(), element);
- fvElemGeomPtr_ = &fvGeometry;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
-
- residual_.resize(fvGeometry_().numVertices);
- residual_ = 0;
-
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = 0;
- curVolVarsPtr_ = &curVolVars;
- asImp_().evalFluxes_();
- }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \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
- * \param curVolVars The volume averaged variables for all
- * sub-control volumes of the element at the current
- * time level
- * \param bcTypes The types of the boundary conditions for all
- * vertices of the element
- */
- void eval(const Element &element,
- const FVElementGeometry &fvGeometry,
- const ElementVolumeVariables &prevVolVars,
- const ElementVolumeVariables &curVolVars,
- const ElementBoundaryTypes &bcTypes)
- {
- Valgrind::CheckDefined(prevVolVars);
- Valgrind::CheckDefined(curVolVars);
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry.numVertices; i++) {
- prevVolVars[i].checkDefined();
- curVolVars[i].checkDefined();
- }
-#endif // HAVE_VALGRIND
-
- elemPtr_ = &element;
- fvElemGeomPtr_ = &fvGeometry;
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = &prevVolVars;
- curVolVarsPtr_ = &curVolVars;
-
- // resize the vectors for all terms
- int numVerts = fvGeometry_().numVertices;
- residual_.resize(numVerts);
- storageTerm_.resize(numVerts);
-
- residual_ = 0.0;
- storageTerm_ = 0.0;
-
- asImp_().evalFluxes_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
-#endif // HAVE_VALGRIND
-
- asImp_().evalVolumeTerms_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++) {
- Valgrind::CheckDefined(residual_[i]);
- }
-#endif // HAVE_VALGRIND
-
- // evaluate the boundary conditions
- asImp_().evalBoundary_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
-#endif // HAVE_VALGRIND
- }
-
- /*!
- * \brief Returns the local residual for all sub-control
- * volumes of the element.
- */
- const ElementSolutionVector &residual() const
- { return residual_; }
-
- /*!
- * \brief Returns the local residual for a given sub-control
- * volume of the element.
- *
- * \param scvIdx The local index of the sub-control volume
- * (i.e. the element's local vertex index)
- */
- const PrimaryVariables &residual(const int scvIdx) const
- { return residual_[scvIdx]; }
-
- /*!
- * \brief Returns the storage term for all sub-control volumes of the
- * element.
- */
- const ElementSolutionVector &storageTerm() const
- { return storageTerm_; }
-
- /*!
- * \brief Returns the storage term for a given sub-control volumes
- * of the element.
- */
- const PrimaryVariables &storageTerm(const int scvIdx) const
- { return storageTerm_[scvIdx]; }
-
protected:
- Implementation &asImp_()
- {
- assert(static_cast<Implementation*>(this) != 0);
- return *static_cast<Implementation*>(this);
- }
-
- const Implementation &asImp_() const
- {
- assert(static_cast<const Implementation*>(this) != 0);
- return *static_cast<const Implementation*>(this);
- }
-
- /*!
- * \brief Evaluate the boundary conditions
- * of the current element.
- */
- void evalBoundary_()
- {
- if (bcTypes_().hasNeumann() || bcTypes_().hasOutflow())
- asImp_().evalBoundaryFluxes_();
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
- Valgrind::CheckDefined(residual_[i]);
-#endif // HAVE_VALGRIND
-
- if (bcTypes_().hasDirichlet())
- asImp_().evalDirichlet_();
- }
-
/*!
* \brief Set the values of the Dirichlet boundary control volumes
* of the current element.
@@ -337,14 +81,14 @@ protected:
void evalDirichlet_()
{
PrimaryVariables dirichletValues(0);
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; ++scvIdx) {
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
+ for (int scvIdx = 0; scvIdx < this->fvGeometry_().numVertices; ++scvIdx) {
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
if (bcTypes.hasDirichlet()) {
// ask the problem for the dirichlet values
- const VertexPointer vPtr = element_().template subEntity<dim>(scvIdx);
+ const VertexPointer vPtr = this->element_().template subEntity<dim>(scvIdx);
Valgrind::SetUndefined(dirichletValues);
- asImp_().problem_().dirichlet(dirichletValues, *vPtr);
+ this->asImp_().problem_().dirichlet(dirichletValues, *vPtr);
// set the dirichlet conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
@@ -353,10 +97,10 @@ protected:
assert(0 <= pvIdx && pvIdx < numEq);
Valgrind::CheckDefined(dirichletValues[pvIdx]);
- residual_[scvIdx][eqIdx] =
- curPriVar_(scvIdx, pvIdx) - dirichletValues[pvIdx];
+ this->residual_[scvIdx][eqIdx] =
+ this->curPriVar_(scvIdx, pvIdx) - dirichletValues[pvIdx];
- storageTerm_[scvIdx][eqIdx] = 0.0;
+ this->storageTerm_[scvIdx][eqIdx] = 0.0;
}
}
}
@@ -369,11 +113,11 @@ protected:
*/
void evalBoundaryFluxes_()
{
- Dune::GeometryType geoType = element_().geometry().type();
+ Dune::GeometryType geoType = this->element_().geometry().type();
const ReferenceElement &refElement = ReferenceElements::general(geoType);
- IntersectionIterator isIt = gridView_().ibegin(element_());
- const IntersectionIterator &endIt = gridView_().iend(element_());
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ const IntersectionIterator &endIt = this->gridView_().iend(this->element_());
for (; isIt != endIt; ++isIt)
{
// handle only faces on the boundary
@@ -392,19 +136,17 @@ protected:
dim);
int boundaryFaceIdx =
- fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
+ this->fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
// add the residual of all vertices of the boundary
// segment
- asImp_().evalNeumannSegment_(isIt,
- scvIdx,
- boundaryFaceIdx);
+ this->asImp_().evalNeumannSegment_(isIt,
+ 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,
- scvIdx,
- boundaryFaceIdx);
+ this->asImp_().evalOutflowSegment_(isIt,
+ scvIdx,
+ boundaryFaceIdx);
}
}
}
@@ -420,28 +162,28 @@ protected:
{
// temporary vector to store the neumann boundary fluxes
PrimaryVariables neumannFlux(0.0);
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
// deal with neumann boundaries
if (bcTypes.hasNeumann()) {
Valgrind::SetUndefined(neumannFlux);
- problem_().boxSDNeumann(neumannFlux,
- element_(),
- fvGeometry_(),
- *isIt,
- scvIdx,
- boundaryFaceIdx,
- curVolVars_());
+ this->problem_().boxSDNeumann(neumannFlux,
+ this->element_(),
+ this->fvGeometry_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx,
+ this->curVolVars_());
neumannFlux *=
- fvGeometry_().boundaryFace[boundaryFaceIdx].area
- * curVolVars_(scvIdx).extrusionFactor();
+ this->fvGeometry_().boundaryFace[boundaryFaceIdx].area
+ * this->curVolVars_(scvIdx).extrusionFactor();
Valgrind::CheckDefined(neumannFlux);
// set the neumann conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!bcTypes.isNeumann(eqIdx))
continue;
- residual_[scvIdx][eqIdx] += neumannFlux[eqIdx];
+ this->residual_[scvIdx][eqIdx] += neumannFlux[eqIdx];
}
}
}
@@ -464,7 +206,7 @@ protected:
{
//calculate outflow fluxes
PrimaryVariables values(0.0);
- asImp_().computeFlux(values, boundaryFaceIdx, true);
+ this->asImp_().computeFlux(values, boundaryFaceIdx, true);
Valgrind::CheckDefined(values);
for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
@@ -485,20 +227,20 @@ protected:
{
// calculate the mass flux over the faces and subtract
// it from the local rates
- for (int scvfIdx = 0; scvfIdx < fvGeometry_().numEdges; scvfIdx++)
+ for (int scvfIdx = 0; scvfIdx < this->fvGeometry_().numEdges; scvfIdx++)
{
- int i = fvGeometry_().subContVolFace[scvfIdx].i;
- int j = fvGeometry_().subContVolFace[scvfIdx].j;
+ int i = this->fvGeometry_().subContVolFace[scvfIdx].i;
+ int j = this->fvGeometry_().subContVolFace[scvfIdx].j;
PrimaryVariables flux;
Valgrind::SetUndefined(flux);
- asImp_().computeFlux(flux, scvfIdx);
+ this->asImp_().computeFlux(flux, scvfIdx);
Valgrind::CheckDefined(flux);
Scalar extrusionFactor =
- (curVolVars_(i).extrusionFactor()
- + curVolVars_(j).extrusionFactor())
+ (this->curVolVars_(i).extrusionFactor()
+ + this->curVolVars_(j).extrusionFactor())
/ 2;
flux *= extrusionFactor;
@@ -516,225 +258,13 @@ protected:
// of sub-control volume i and _INTO_ sub-control volume
// j, we need to add the flux to finite volume i and
// subtract it from finite volume j
- residual_[i] += flux;
- residual_[j] -= flux;
+ this->residual_[i] += flux;
+ this->residual_[j] -= flux;
}
}
-
- /*!
- * \brief Set the local residual to the storage terms of all
- * sub-control volumes of the current element.
- */
- void evalStorage_()
- {
- storageTerm_.resize(fvGeometry_().numVertices);
- storageTerm_ = 0;
-
- // calculate the amount of conservation each quantity inside
- // all sub control volumes
- 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]);
- }
- }
-
- /*!
- * \brief Add the change in the storage terms and the source term
- * to the local residual of all sub-control volumes of the
- * current element.
- */
- void evalVolumeTerms_()
- {
- // evaluate the volume terms (storage + source terms)
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++)
- {
- Scalar extrusionFactor =
- curVolVars_(scvIdx).extrusionFactor();
-
- PrimaryVariables values(0.0);
-
- // mass balance within the element. this is the
- // \f$\frac{m}{\partial t}\f$ term if using implicit
- // euler as time discretization.
- //
- // TODO (?): we might need a more explicit way for
- // doing the time discretization...
- 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_[scvIdx] += storageTerm_[scvIdx];
-
- // subtract the source term from the local rate
- Valgrind::SetUndefined(values);
- asImp_().computeSource(values, scvIdx);
- Valgrind::CheckDefined(values);
- values *= fvGeometry_().subContVol[scvIdx].volume * extrusionFactor;
- residual_[scvIdx] -= values;
-
- // make sure that only defined quantities were used
- // to calculate the residual.
- Valgrind::CheckDefined(residual_[scvIdx]);
- }
- }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- const Problem &problem_() const
- { return *problemPtr_; };
-
- /*!
- * \brief Returns a reference to the model.
- */
- const Model &model_() const
- { return problem_().model(); };
-
- /*!
- * \brief Returns a reference to the vertex mapper.
- */
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); };
-
- /*!
- * \brief Returns a reference to the grid view.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns a reference to the current element.
- */
- const Element &element_() const
- {
- Valgrind::CheckDefined(elemPtr_);
- return *elemPtr_;
- }
- /*!
- * \brief Returns a reference to the current element's finite
- * volume geometry.
- */
- const FVElementGeometry &fvGeometry_() const
- {
- Valgrind::CheckDefined(fvElemGeomPtr_);
- return *fvElemGeomPtr_;
- }
-
- /*!
- * \brief Returns a reference to the primary variables of
- * the last time step of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &prevPriVars_(const int i) const
- {
- return prevVolVars_(i).priVars();
- }
-
- /*!
- * \brief Returns a reference to the primary variables of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &curPriVars_(const int i) const
- {
- return curVolVars_(i).priVars();
- }
-
- /*!
- * \brief Returns the j'th primary of the i'th sub-control volume
- * of the current element.
- */
- Scalar curPriVar_(const int i, const int j) const
- {
- return curVolVars_(i).priVar(j);
- }
-
- /*!
- * \brief Returns a reference to the current volume variables of
- * all sub-control volumes of the current element.
- */
- const ElementVolumeVariables &curVolVars_() const
- {
- Valgrind::CheckDefined(curVolVarsPtr_);
- return *curVolVarsPtr_;
- }
-
- /*!
- * \brief Returns a reference to the volume variables of the i-th
- * sub-control volume of the current element.
- */
- const VolumeVariables &curVolVars_(const int i) const
- {
- return curVolVars_()[i];
- }
-
- /*!
- * \brief Returns a reference to the previous time step's volume
- * variables of all sub-control volumes of the current
- * element.
- */
- const ElementVolumeVariables &prevVolVars_() const
- {
- Valgrind::CheckDefined(prevVolVarsPtr_);
- return *prevVolVarsPtr_;
- }
-
- /*!
- * \brief Returns a reference to the previous time step's volume
- * variables of the i-th sub-control volume of the current
- * element.
- */
- const VolumeVariables &prevVolVars_(const int i) const
- {
- return prevVolVars_()[i];
- }
-
- /*!
- * \brief Returns a reference to the boundary types of all
- * sub-control volumes of the current element.
- */
- const ElementBoundaryTypes &bcTypes_() const
- {
- Valgrind::CheckDefined(bcTypesPtr_);
- return *bcTypesPtr_;
- }
-
- /*!
- * \brief Returns a reference to the boundary types of the i-th
- * sub-control volume of the current element.
- */
- const BoundaryTypes &bcTypes_(const int i) const
- {
- return bcTypes_()[i];
- }
-
-protected:
- ElementSolutionVector storageTerm_;
- ElementSolutionVector residual_;
-
- // The problem we would like to solve
- Problem *problemPtr_;
-
- const Element *elemPtr_;
- const FVElementGeometry *fvElemGeomPtr_;
-
- // current and previous secondary variables for the element
- const ElementVolumeVariables *prevVolVarsPtr_;
- const ElementVolumeVariables *curVolVarsPtr_;
-
- const ElementBoundaryTypes *bcTypesPtr_;
+ const VertexMapper &vertexMapper_() const
+ { return this->problem_().vertexMapper(); };
};
}
diff --git a/dumux/implicit/cellcentered/cclocalresidual.hh b/dumux/implicit/cellcentered/cclocalresidual.hh
index cdc914e8b2d0b9c446f6486eb2e8a92dcfca1881..abe2f9f0b53319fc00e1f2c0c5806eb0672132cc 100644
--- a/dumux/implicit/cellcentered/cclocalresidual.hh
+++ b/dumux/implicit/cellcentered/cclocalresidual.hh
@@ -31,6 +31,7 @@
#include <dune/grid/common/geometry.hh>
#include <dumux/common/valgrind.hh>
+#include <dumux/implicit/common/implicitlocalresidual.hh>
#include "ccproperties.hh"
@@ -39,18 +40,16 @@ namespace Dumux
/*!
* \ingroup CCModel
* \ingroup CCLocalResidual
- * \brief Element-wise calculation of the residual matrix for models
- * based on the box scheme.
+ * \brief Element-wise calculation of the residual for models
+ * based on the fully implicit cell-centered scheme.
*
* \todo Please doc me more!
*/
template<class TypeTag>
-class CCLocalResidual
+class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+ typedef ImplicitLocalResidual<TypeTag> ParentType;
+ friend class ImplicitLocalResidual<TypeTag>;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -62,302 +61,17 @@ private:
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
// copying the local residual class is not a good idea
CCLocalResidual(const CCLocalResidual &);
public:
- CCLocalResidual()
- { }
-
- ~CCLocalResidual()
+ CCLocalResidual() : ParentType()
{ }
- /*!
- * \brief Initialize the local residual.
- *
- * 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
- * solved.
- */
- void init(Problem &prob)
- { problemPtr_ = &prob; }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- */
- void eval(const Element &element)
- {
- FVElementGeometry fvElemGeom;
-
- fvElemGeom.update(gridView_(), element);
- fvElemGeomPtr_ = &fvElemGeom;
-
- ElementVolumeVariables volVarsPrev, volVarsCur;
- // update the hints
- model_().setHints(element, volVarsPrev, volVarsCur);
-
- volVarsPrev.update(problem_(),
- element,
- fvGeometry_(),
- true /* oldSol? */);
- volVarsCur.update(problem_(),
- element,
- fvGeometry_(),
- false /* oldSol? */);
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
-
- // 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(element);
-
- asImp_().eval(element, fvGeometry_(), volVarsPrev, volVarsCur, bcTypes);
- }
-
- /*!
- * \brief Compute the storage term for the current solution.
- *
- * This can be used to figure out how much of each conservation
- * quantity is inside the element.
- *
- * \param element The DUNE Codim<0> entity for which the storage
- * term ought to be calculated
- */
- void evalStorage(const Element &element)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvElemGeom;
- fvElemGeom.update(gridView_(), element);
- fvElemGeomPtr_ = &fvElemGeom;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
- bcTypesPtr_ = &bcTypes;
-
- // no previous volume variables!
- prevVolVarsPtr_ = 0;
-
- ElementVolumeVariables volVars;
-
- // update the hints
- model_().setHints(element, volVars);
-
- // calculate volume current variables
- volVars.update(problem_(), element, fvGeometry_(), false);
- curVolVarsPtr_ = &volVars;
-
- asImp_().evalStorage_();
- }
-
- /*!
- * \brief Compute the flux term for the current solution.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param curVolVars The volume averaged variables for all
- * sub-contol volumes of the element
- */
- void evalFluxes(const Element &element,
- const ElementVolumeVariables &curVolVars)
- {
- elemPtr_ = &element;
-
- FVElementGeometry fvElemGeom;
- fvElemGeom.update(gridView_(), element);
- fvElemGeomPtr_ = &fvElemGeom;
-
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
-
- residual_.resize(1);
- residual_ = 0;
-
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = 0;
- curVolVarsPtr_ = &curVolVars;
- asImp_().evalFluxes_();
- }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero.
- *
- * \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 prevVolVars The volume averaged variables for all
- * sub-control volumes of the element at the previous
- * time level
- * \param curVolVars The volume averaged variables for all
- * sub-control volumes of the element at the current
- * time level
- * \param bcTypes The types of the boundary conditions for all
- * vertices of the element
- */
- void eval(const Element &element,
- const FVElementGeometry &fvElemGeom,
- const ElementVolumeVariables &prevVolVars,
- const ElementVolumeVariables &curVolVars,
- const ElementBoundaryTypes &bcTypes)
- {
- Valgrind::CheckDefined(prevVolVars);
- Valgrind::CheckDefined(curVolVars);
-
-#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvElemGeom.numNeighbors; i++) {
- prevVolVars[i].checkDefined();
- curVolVars[i].checkDefined();
- }
-#endif // HAVE_VALGRIND
-
- elemPtr_ = &element;
- fvElemGeomPtr_ = &fvElemGeom;
- bcTypesPtr_ = &bcTypes;
- prevVolVarsPtr_ = &prevVolVars;
- curVolVarsPtr_ = &curVolVars;
-
- // resize the vectors for all terms
- residual_.resize(1);
- storageTerm_.resize(1);
-
- residual_ = 0.0;
- storageTerm_ = 0.0;
-
- asImp_().evalFluxes_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- Valgrind::CheckDefined(residual_[0]);
-#endif // HAVE_VALGRIND
-
- asImp_().evalVolumeTerms_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- Valgrind::CheckDefined(residual_[0]);
-#endif // HAVE_VALGRIND
-
- // evaluate the boundary conditions
- asImp_().evalBoundary_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- Valgrind::CheckDefined(residual_[0]);
-#endif // HAVE_VALGRIND
- }
-
- /*!
- * \brief Add Outflow boundary conditions for a single sub-control
- * volume face to the local residual.
- *
- * \note This is so far an empty method doing not more than
- * nothing. A beta version is available for the 2p2c and
- * 2p2cni model. There you can find a sample implementation.
- */
- void evalOutflowSegment(const IntersectionIterator &isIt,
- int scvIdx,
- int boundaryFaceIdx)
- {}
-
- /*!
- * \brief Returns the local residual for all sub-control
- * volumes of the element.
- */
- const ElementSolutionVector &residual() const
- { return residual_; }
-
- /*!
- * \brief Returns the local residual for a given sub-control
- * volume of the element.
- *
- * \param scvIdx The local index of the sub-control volume
- * (i.e. the element's local vertex index)
- */
- const PrimaryVariables &residual(int scvIdx) const
- { return residual_[scvIdx]; }
-
- /*!
- * \brief Returns the storage term for all sub-control volumes of the
- * element.
- */
- const ElementSolutionVector &storageTerm() const
- { return storageTerm_; }
-
- /*!
- * \brief Returns the storage term for a given sub-control volumes
- * of the element.
- */
- const PrimaryVariables &storageTerm(int scvIdx) const
- { return storageTerm_[scvIdx]; }
-
- void evalOutflowSegment(const IntersectionIterator &isIt)
- {
- const BoundaryTypes &bcTypes = this->bcTypes_(0);
-
- // deal with outflow boundaries
- if (bcTypes.hasOutflow())
- {
- PrimaryVariables values(0.0);
- asImp_().computeFlux(values,
- /*boundaryFaceIdx=*/isIt->indexInInside(),
- true);
- Valgrind::CheckDefined(values);
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- if (!bcTypes.isOutflow(eqIdx) )
- continue;
- this->residual_[0][eqIdx] += values[eqIdx];
- }
- }
- }
protected:
- Implementation &asImp_()
- {
- assert(static_cast<Implementation*>(this) != 0);
- return *static_cast<Implementation*>(this);
- }
-
- const Implementation &asImp_() const
- {
- assert(static_cast<const Implementation*>(this) != 0);
- return *static_cast<const Implementation*>(this);
- }
-
- /*!
- * \brief Evaluate the boundary conditions
- * of the current element.
- */
- void evalBoundary_()
- {
- if (bcTypes_().hasNeumann() || bcTypes_().hasOutflow())
- asImp_().evalBoundaryFluxes_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
- Valgrind::CheckDefined(residual_[0]);
-#endif // HAVE_VALGRIND
-
- if (bcTypes_().hasDirichlet())
- asImp_().evalDirichlet_();
- }
/*!
* \brief Set the values of the Dirichlet boundary control volumes
@@ -365,15 +79,15 @@ protected:
*/
void evalDirichlet_()
{
- PrimaryVariables tmp(0);
+ PrimaryVariables dirichletValues(0);
- const BoundaryTypes &bcTypes = bcTypes_(0);
+ const BoundaryTypes &bcTypes = this->bcTypes_(0);
if (! bcTypes.hasDirichlet())
return;
- Valgrind::SetUndefined(tmp);
- // HACK: ask for Dirichlet value at element center
- asImp_().problem_().dirichletAtPos(tmp, element_().geometry().center());
+ Valgrind::SetUndefined(dirichletValues);
+ // HACK: ask for Dirichlet value at element center
+ this->asImp_().problem_().dirichletAtPos(dirichletValues, this->element_().geometry().center());
// set the dirichlet conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
@@ -381,12 +95,12 @@ protected:
continue;
int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
assert(0 <= pvIdx && pvIdx < numEq);
- Valgrind::CheckDefined(tmp[pvIdx]);
+ Valgrind::CheckDefined(dirichletValues[pvIdx]);
- residual_[0][eqIdx] =
- (curPriVar_(0, pvIdx) - tmp[pvIdx]);
+ this->residual_[0][eqIdx] =
+ (this->curPriVar_(0, pvIdx) - dirichletValues[pvIdx]);
- storageTerm_[0][eqIdx] = 0.0;
+ this->storageTerm_[0][eqIdx] = 0.0;
}
}
@@ -396,8 +110,8 @@ protected:
*/
void evalBoundaryFluxes_()
{
- IntersectionIterator isIt = gridView_().ibegin(element_());
- const IntersectionIterator &endIt = gridView_().iend(element_());
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ const IntersectionIterator &endIt = this->gridView_().iend(this->element_());
for (; isIt != endIt; ++isIt)
{
// handle only faces on the boundary
@@ -406,18 +120,15 @@ protected:
// add the residual of all vertices of the boundary
// segment
- asImp_().evalNeumannSegment_(isIt);
-
+ this->asImp_().evalNeumannSegment_(isIt);
+
// 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);
+ this->asImp_().evalOutflowSegment_(isIt);
}
}
/*!
- * \brief Add Neumann boundary conditions for a single sub-control
- * volume face to the local residual.
+ * \brief Add Neumann boundary conditions for a single intersection
*/
void evalNeumannSegment_(const IntersectionIterator &isIt)
{
@@ -425,42 +136,66 @@ protected:
PrimaryVariables values(0.0);
BoundaryTypes bcTypes;
- problem_().boundaryTypes(bcTypes, *isIt);
+ this->problem_().boundaryTypes(bcTypes, *isIt);
// deal with neumann boundaries
if (bcTypes.hasNeumann()) {
Valgrind::SetUndefined(values);
- problem_().boxSDNeumann(values,
- element_(),
- fvGeometry_(),
- *isIt,
- /*scvIdx=*/0,
- /*boundaryFaceIdx=*/isIt->indexInInside(),
- curVolVars_());
+ this->problem_().boxSDNeumann(values,
+ this->element_(),
+ this->fvGeometry_(),
+ *isIt,
+ /*scvIdx=*/0,
+ /*boundaryFaceIdx=*/isIt->indexInInside(),
+ this->curVolVars_());
values *= isIt->geometry().volume()
- * curVolVars_(0).extrusionFactor();
+ * this->curVolVars_(0).extrusionFactor();
Valgrind::CheckDefined(values);
// set the neumann conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
if (!bcTypes.isNeumann(eqIdx))
continue;
- residual_[0][eqIdx] += values[eqIdx];
+ this->residual_[0][eqIdx] += values[eqIdx];
+ }
+ }
+ }
+
+ /*!
+ * \brief Add outflow boundary conditions for a single intersection
+ */
+ void evalOutflowSegment_(const IntersectionIterator &isIt)
+ {
+ const BoundaryTypes &bcTypes = this->bcTypes_(0);
+
+ // deal with outflow boundaries
+ if (bcTypes.hasOutflow())
+ {
+ PrimaryVariables values(0.0);
+ this->asImp_().computeFlux(values,
+ /*boundaryFaceIdx=*/isIt->indexInInside(),
+ true);
+ Valgrind::CheckDefined(values);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (!bcTypes.isOutflow(eqIdx) )
+ continue;
+ this->residual_[0][eqIdx] += values[eqIdx];
}
}
}
/*!
- * \brief Add the flux terms to the local residual of all
- * sub-control volumes of the current element.
+ * \brief Add the flux terms to the local residual of the current element
*/
void evalFluxes_()
{
// calculate the mass flux over the faces and subtract
// it from the local rates
int faceIdx = -1;
- IntersectionIterator isIt = gridView_().ibegin(element_());
- IntersectionIterator isEndIt = gridView_().iend(element_());
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ IntersectionIterator isEndIt = this->gridView_().iend(this->element_());
for (; isIt != isEndIt; ++isIt) {
if (!isIt->neighbor())
continue;
@@ -469,240 +204,14 @@ protected:
PrimaryVariables flux;
Valgrind::SetUndefined(flux);
- asImp_().computeFlux(flux, faceIdx);
+ this->asImp_().computeFlux(flux, faceIdx);
Valgrind::CheckDefined(flux);
- flux *= curVolVars_(0).extrusionFactor();
-
- // The balance equation for a finite volume is:
- //
- // dStorage/dt = Flux + Source
- //
- // where the 'Flux' and the 'Source' terms represent the
- // mass per second which _ENTER_ the finite
- // volume. Re-arranging this, we get
- //
- // dStorage/dt - Source - Flux = 0
- //
- // Since the flux calculated by computeFlux() goes _OUT_
- // of sub-control volume i and _INTO_ sub-control volume
- // j, we need to add the flux to finite volume i and
- // subtract it from finite volume j
- residual_[0] += flux;
- }
- }
-
- /*!
- * \brief Set the local residual to the storage terms of all
- * sub-control volumes of the current element.
- */
- void evalStorage_()
- {
- storageTerm_.resize(1);
- storageTerm_ = 0;
-
- // calculate the amount of conservation each quantity inside
- // all sub control volumes
- Valgrind::SetUndefined(storageTerm_[0]);
- asImp_().computeStorage(storageTerm_[0], 0, /*isOldSol=*/false);
- storageTerm_[0] *= fvGeometry_().subContVol[0].volume
- * curVolVars_(0).extrusionFactor();
- Valgrind::CheckDefined(storageTerm_[0]);
- }
-
- /*!
- * \brief Add the change in the storage terms and the source term
- * to the local residual of all sub-control volumes of the
- * current element.
- */
- void evalVolumeTerms_()
- {
- // evaluate the volume terms (storage + source terms)
- for (int i=0; i < 1; i++)
- {
- Scalar extrusionFactor =
- curVolVars_(i).extrusionFactor();
-
- PrimaryVariables tmp(0.);
-
- // mass balance within the element. this is the
- // \f$\frac{m}{\partial t}\f$ term if using implicit
- // euler as time discretization.
- //
- // 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);
+ flux *= this->curVolVars_(0).extrusionFactor();
- storageTerm_[i] -= tmp;
- storageTerm_[i] *=
- fvGeometry_().subContVol[i].volume
- / problem_().timeManager().timeStepSize()
- * extrusionFactor;
- residual_[i] += storageTerm_[i];
-
- // subtract the source term from the local rate
- Valgrind::SetUndefined(tmp);
- asImp_().computeSource(tmp, i);
- Valgrind::CheckDefined(tmp);
- tmp *= fvGeometry_().subContVol[i].volume * extrusionFactor;
- residual_[i] -= tmp;
-
- // make sure that only defined quantities were used
- // to calculate the residual.
- Valgrind::CheckDefined(residual_[i]);
+ this->residual_[0] += flux;
}
}
-
- /*!
- * \brief Returns a reference to the problem.
- */
- const Problem &problem_() const
- { return *problemPtr_; };
-
- /*!
- * \brief Returns a reference to the model.
- */
- const Model &model_() const
- { return problem_().model(); };
-
- /*!
- * \brief Returns a reference to the vertex mapper.
- */
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); };
-
- /*!
- * \brief Returns a reference to the grid view.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns a reference to the current element.
- */
- const Element &element_() const
- {
- Valgrind::CheckDefined(elemPtr_);
- return *elemPtr_;
- }
-
- /*!
- * \brief Returns a reference to the current element's finite
- * volume geometry.
- */
- const FVElementGeometry &fvGeometry_() const
- {
- Valgrind::CheckDefined(fvElemGeomPtr_);
- return *fvElemGeomPtr_;
- }
-
- /*!
- * \brief Returns a reference to the primary variables of
- * the last time step of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &prevPriVars_(int i) const
- {
- return prevVolVars_(i).priVars();
- }
-
- /*!
- * \brief Returns a reference to the primary variables of the i'th
- * sub-control volume of the current element.
- */
- const PrimaryVariables &curPriVars_(int i) const
- {
- return curVolVars_(i).priVars();
- }
-
- /*!
- * \brief Returns the j'th primary of the i'th sub-control volume
- * of the current element.
- */
- Scalar curPriVar_(int i, int j) const
- {
- return curVolVars_(i).priVar(j);
- }
-
- /*!
- * \brief Returns a reference to the current volume variables of
- * all sub-control volumes of the current element.
- */
- const ElementVolumeVariables &curVolVars_() const
- {
- Valgrind::CheckDefined(curVolVarsPtr_);
- return *curVolVarsPtr_;
- }
-
- /*!
- * \brief Returns a reference to the volume variables of the i-th
- * sub-control volume of the current element.
- */
- const VolumeVariables &curVolVars_(int i) const
- {
- return curVolVars_()[i];
- }
-
- /*!
- * \brief Returns a reference to the previous time step's volume
- * variables of all sub-control volumes of the current
- * element.
- */
- const ElementVolumeVariables &prevVolVars_() const
- {
- Valgrind::CheckDefined(prevVolVarsPtr_);
- return *prevVolVarsPtr_;
- }
-
- /*!
- * \brief Returns a reference to the previous time step's volume
- * variables of the i-th sub-control volume of the current
- * element.
- */
- const VolumeVariables &prevVolVars_(int i) const
- {
- return prevVolVars_()[i];
- }
-
- /*!
- * \brief Returns a reference to the boundary types of all
- * sub-control volumes of the current element.
- */
- const ElementBoundaryTypes &bcTypes_() const
- {
- Valgrind::CheckDefined(bcTypesPtr_);
- return *bcTypesPtr_;
- }
-
- /*!
- * \brief Returns a reference to the boundary types of the i-th
- * sub-control volume of the current element.
- */
- const BoundaryTypes &bcTypes_(int i) const
- {
- return bcTypes_()[i];
- }
-
-protected:
- ElementSolutionVector storageTerm_;
- ElementSolutionVector residual_;
-
- // The problem we would like to solve
- Problem *problemPtr_;
-
- const Element *elemPtr_;
- const FVElementGeometry *fvElemGeomPtr_;
-
- // current and previous secondary variables for the element
- const ElementVolumeVariables *prevVolVarsPtr_;
- const ElementVolumeVariables *curVolVarsPtr_;
-
- const ElementBoundaryTypes *bcTypesPtr_;
};
}
diff --git a/dumux/implicit/common/implicitlocalresidual.hh b/dumux/implicit/common/implicitlocalresidual.hh
index 53eaa0bf554990078effe06e67ccab58387bf4e9..f020e8bbd1a32fd1a387bc4cfda69a1a7a56a0fd 100644
--- a/dumux/implicit/common/implicitlocalresidual.hh
+++ b/dumux/implicit/common/implicitlocalresidual.hh
@@ -20,8 +20,8 @@
* \file
* \brief Calculates the residual of models based on the box scheme element-wise.
*/
-#ifndef DUMUX_BOX_LOCAL_RESIDUAL_HH
-#define DUMUX_BOX_LOCAL_RESIDUAL_HH
+#ifndef DUMUX_IMPLICIT_LOCAL_RESIDUAL_HH
+#define DUMUX_IMPLICIT_LOCAL_RESIDUAL_HH
#include <dune/istl/matrix.hh>
#include <dune/grid/common/geometry.hh>
@@ -33,15 +33,15 @@
namespace Dumux
{
/*!
- * \ingroup BoxModel
- * \ingroup BoxLocalResidual
+ * \ingroup ImplicitModel
+ * \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the residual matrix for models
- * based on the box scheme.
+ * using a fully implicit discretization.
*
* \todo Please doc me more!
*/
template<class TypeTag>
-class BoxLocalResidual
+class ImplicitLocalResidual
{
private:
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
@@ -49,22 +49,10 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, NumEq),
- dim = GridView::dimension
- };
-
typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
- typedef typename GridView::Grid::ctype CoordScalar;
- typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
- typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
-
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
@@ -73,13 +61,13 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
// copying the local residual class is not a good idea
- BoxLocalResidual(const BoxLocalResidual &);
+ ImplicitLocalResidual(const ImplicitLocalResidual &);
public:
- BoxLocalResidual()
+ ImplicitLocalResidual()
{ }
- ~BoxLocalResidual()
+ ~ImplicitLocalResidual()
{ }
/*!
@@ -190,7 +178,7 @@ public:
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeometry_());
- residual_.resize(fvGeometry_().numVertices);
+ residual_.resize(fvGeometry_().numSCV);
residual_ = 0;
bcTypesPtr_ = &bcTypes;
@@ -225,7 +213,7 @@ public:
Valgrind::CheckDefined(curVolVars);
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry.numVertices; i++) {
+ for (int i = 0; i < prevVolVars.size(); i++) {
prevVolVars[i].checkDefined();
curVolVars[i].checkDefined();
}
@@ -238,9 +226,9 @@ public:
curVolVarsPtr_ = &curVolVars;
// resize the vectors for all terms
- int numVerts = fvGeometry_().numVertices;
- residual_.resize(numVerts);
- storageTerm_.resize(numVerts);
+ int numSCV = fvGeometry_().numSCV;
+ residual_.resize(numSCV);
+ storageTerm_.resize(numSCV);
residual_ = 0.0;
storageTerm_ = 0.0;
@@ -248,14 +236,14 @@ public:
asImp_().evalFluxes_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
+ for (int i=0; i < fvGeometry_().numSCV; i++)
Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
asImp_().evalVolumeTerms_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++) {
+ for (int i=0; i < fvGeometry_().numSCV; i++) {
Valgrind::CheckDefined(residual_[i]);
}
#endif // HAVE_VALGRIND
@@ -264,7 +252,7 @@ public:
asImp_().evalBoundary_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
+ for (int i=0; i < fvGeometry_().numSCV; i++)
Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
}
@@ -322,7 +310,7 @@ protected:
if (bcTypes_().hasNeumann() || bcTypes_().hasOutflow())
asImp_().evalBoundaryFluxes_();
#if !defined NDEBUG && HAVE_VALGRIND
- for (int i=0; i < fvGeometry_().numVertices; i++)
+ for (int i=0; i < fvGeometry_().numSCV; i++)
Valgrind::CheckDefined(residual_[i]);
#endif // HAVE_VALGRIND
@@ -330,209 +318,18 @@ protected:
asImp_().evalDirichlet_();
}
- /*!
- * \brief Set the values of the Dirichlet boundary control volumes
- * of the current element.
- */
- void evalDirichlet_()
- {
- PrimaryVariables dirichletValues(0);
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; ++scvIdx) {
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
-
- if (bcTypes.hasDirichlet()) {
- // ask the problem for the dirichlet values
- 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) {
- if (bcTypes.isDirichlet(eqIdx)) {
- int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
- assert(0 <= pvIdx && pvIdx < numEq);
- Valgrind::CheckDefined(dirichletValues[pvIdx]);
-
- residual_[scvIdx][eqIdx] =
- curPriVar_(scvIdx, pvIdx) - dirichletValues[pvIdx];
-
- storageTerm_[scvIdx][eqIdx] = 0.0;
- }
- }
- }
- }
- }
-
- /*!
- * \brief Add all Neumann and outflow boundary conditions to the local
- * residual.
- */
- void evalBoundaryFluxes_()
- {
- Dune::GeometryType geoType = element_().geometry().type();
- const ReferenceElement &refElement = ReferenceElements::general(geoType);
-
- IntersectionIterator isIt = gridView_().ibegin(element_());
- const IntersectionIterator &endIt = gridView_().iend(element_());
- for (; isIt != endIt; ++isIt)
- {
- // handle only faces on the boundary
- if (isIt->boundary()) {
- // Assemble the boundary for all vertices of the current
- // face
- int faceIdx = isIt->indexInInside();
- int numFaceVerts = refElement.size(faceIdx, 1, dim);
- for (int faceVertIdx = 0;
- faceVertIdx < numFaceVerts;
- ++faceVertIdx)
- {
- int scvIdx = refElement.subEntity(faceIdx,
- 1,
- faceVertIdx,
- dim);
-
- int boundaryFaceIdx =
- fvGeometry_().boundaryFaceIndex(faceIdx, faceVertIdx);
-
- // add the residual of all vertices of the boundary
- // segment
- asImp_().evalNeumannSegment_(isIt,
- 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,
- scvIdx,
- boundaryFaceIdx);
- }
- }
- }
- }
-
- /*!
- * \brief Add Neumann boundary conditions for a single sub-control
- * volume face to the local residual.
- */
- void evalNeumannSegment_(const IntersectionIterator &isIt,
- const int scvIdx,
- const int boundaryFaceIdx)
- {
- // temporary vector to store the neumann boundary fluxes
- PrimaryVariables neumannFlux(0.0);
- const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
-
- // deal with neumann boundaries
- if (bcTypes.hasNeumann()) {
- Valgrind::SetUndefined(neumannFlux);
- problem_().boxSDNeumann(neumannFlux,
- element_(),
- fvGeometry_(),
- *isIt,
- scvIdx,
- boundaryFaceIdx,
- curVolVars_());
- neumannFlux *=
- fvGeometry_().boundaryFace[boundaryFaceIdx].area
- * curVolVars_(scvIdx).extrusionFactor();
- Valgrind::CheckDefined(neumannFlux);
-
- // set the neumann conditions
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- if (!bcTypes.isNeumann(eqIdx))
- continue;
- residual_[scvIdx][eqIdx] += neumannFlux[eqIdx];
- }
- }
- }
-
- /*!
- * \brief Add outflow boundary conditions for a single sub-control
- * volume face to the local residual.
- *
- * \param isIt The intersection iterator of current element
- * \param scvIdx The index of the considered face of the sub-control volume
- * \param boundaryFaceIdx The index of the considered boundary face of the sub control volume
- */
- void evalOutflowSegment_(const IntersectionIterator &isIt,
- const int scvIdx,
- const int boundaryFaceIdx)
- {
- const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
- // deal with outflow boundaries
- if (bcTypes.hasOutflow())
- {
- //calculate outflow fluxes
- PrimaryVariables values(0.0);
- asImp_().computeFlux(values, boundaryFaceIdx, true);
- Valgrind::CheckDefined(values);
-
- for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
- {
- if (!bcTypes.isOutflow(equationIdx) )
- continue;
- // deduce outflow
- this->residual_[scvIdx][equationIdx] += values[equationIdx];
- }
- }
- }
-
- /*!
- * \brief Add the flux terms to the local residual of all
- * sub-control volumes of the current element.
- */
- void evalFluxes_()
- {
- // calculate the mass flux over the faces and subtract
- // it from the local rates
- for (int scvfIdx = 0; scvfIdx < fvGeometry_().numEdges; scvfIdx++)
- {
- int i = fvGeometry_().subContVolFace[scvfIdx].i;
- int j = fvGeometry_().subContVolFace[scvfIdx].j;
-
- PrimaryVariables flux;
-
- Valgrind::SetUndefined(flux);
- asImp_().computeFlux(flux, scvfIdx);
- Valgrind::CheckDefined(flux);
-
- Scalar extrusionFactor =
- (curVolVars_(i).extrusionFactor()
- + curVolVars_(j).extrusionFactor())
- / 2;
- flux *= extrusionFactor;
-
- // The balance equation for a finite volume is:
- //
- // dStorage/dt = Flux + Source
- //
- // where the 'Flux' and the 'Source' terms represent the
- // mass per second which _ENTER_ the finite
- // volume. Re-arranging this, we get
- //
- // dStorage/dt - Source - Flux = 0
- //
- // Since the flux calculated by computeFlux() goes _OUT_
- // of sub-control volume i and _INTO_ sub-control volume
- // j, we need to add the flux to finite volume i and
- // subtract it from finite volume j
- residual_[i] += flux;
- residual_[j] -= flux;
- }
- }
-
/*!
* \brief Set the local residual to the storage terms of all
* sub-control volumes of the current element.
*/
void evalStorage_()
{
- storageTerm_.resize(fvGeometry_().numVertices);
+ storageTerm_.resize(fvGeometry_().numSCV);
storageTerm_ = 0;
// calculate the amount of conservation each quantity inside
// all sub control volumes
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++) {
+ for (int scvIdx = 0; scvIdx < fvGeometry_().numSCV; scvIdx++) {
Valgrind::SetUndefined(storageTerm_[scvIdx]);
asImp_().computeStorage(storageTerm_[scvIdx], scvIdx, /*isOldSol=*/false);
storageTerm_[scvIdx] *=
@@ -550,7 +347,7 @@ protected:
void evalVolumeTerms_()
{
// evaluate the volume terms (storage + source terms)
- for (int scvIdx = 0; scvIdx < fvGeometry_().numVertices; scvIdx++)
+ for (int scvIdx = 0; scvIdx < fvGeometry_().numSCV; scvIdx++)
{
Scalar extrusionFactor =
curVolVars_(scvIdx).extrusionFactor();
@@ -602,12 +399,6 @@ protected:
const Model &model_() const
{ return problem_().model(); };
- /*!
- * \brief Returns a reference to the vertex mapper.
- */
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); };
-
/*!
* \brief Returns a reference to the grid view.
*/
diff --git a/dumux/implicit/common/implicitporousmediaproblem.hh b/dumux/implicit/common/implicitporousmediaproblem.hh
index 0fbde7b6edbbef7299712e892a35a8ce6c6128b3..fab231a7b6486b762d62ef2c0ca6cfd9410c0872 100644
--- a/dumux/implicit/common/implicitporousmediaproblem.hh
+++ b/dumux/implicit/common/implicitporousmediaproblem.hh
@@ -1,7 +1,6 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
- * 2012 by Bernd Flemisch *
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *