From 238545729034208bb372afc3976d57cd7100b98a Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Fri, 30 Nov 2012 12:30:46 +0000
Subject: [PATCH] implicit branch: unify assemblers.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9739 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/box/boxassembler.hh | 668 +++------------------
dumux/implicit/cellcentered/ccassembler.hh | 557 +++--------------
dumux/implicit/common/implicitassembler.hh | 405 ++-----------
3 files changed, 215 insertions(+), 1415 deletions(-)
diff --git a/dumux/implicit/box/boxassembler.hh b/dumux/implicit/box/boxassembler.hh
index d585d360a2..6bbcebb9f7 100644
--- a/dumux/implicit/box/boxassembler.hh
+++ b/dumux/implicit/box/boxassembler.hh
@@ -24,11 +24,7 @@
#ifndef DUMUX_BOX_ASSEMBLER_HH
#define DUMUX_BOX_ASSEMBLER_HH
-#include <dune/grid/common/gridenums.hh>
-
-#include <dumux/implicit/box/boxproperties.hh>
-#include <dumux/linear/vertexborderlistfromgrid.hh>
-#include <dumux/linear/foreignoverlapfrombcrsmatrix.hh>
+#include <dumux/implicit/common/implicitassembler.hh>
#include <dumux/parallel/vertexhandles.hh>
namespace Dumux {
@@ -37,282 +33,26 @@ namespace Dumux {
* \brief An assembler for the global Jacobian matrix for models using the box discretization.
*/
template<class TypeTag>
-class BoxAssembler
+class BoxAssembler : public ImplicitAssembler<TypeTag>
{
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef ImplicitAssembler<TypeTag> ParentType;
+ friend class ImplicitAssembler<TypeTag>;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
-
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-
- enum{ dim = GridView::dimension };
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ enum{ dim = GridView::dimension };
typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
- typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
- typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
-
+public:
+ BoxAssembler(): ParentType() {}
+
+private:
// copying the jacobian assembler is not a good idea
BoxAssembler(const BoxAssembler &);
-public:
- /*!
- * \brief The colors of elements and vertices required for partial
- * Jacobian reassembly.
- */
- enum EntityColor {
- /*!
- * Vertex/element that needs to be reassembled because some
- * relative error is above the tolerance
- */
- Red = 0,
-
- /*!
- * Vertex/element that needs to be reassembled because a
- * neighboring element/vertex is red
- */
- Yellow = 1,
-
- /*!
- * Yellow vertex has only non-green neighbor elements.
- *
- * This means that its relative error is below the tolerance,
- * but its defect can be linearized without any additional
- * cost. This is just an "internal" color which is not used
- * ouside of the jacobian assembler.
- */
- Orange = 2,
-
- /*!
- * Vertex/element that does not need to be reassembled
- */
- Green = 3
- };
-
- BoxAssembler()
- {
- problemPtr_ = 0;
- matrix_ = 0;
-
- // set reassemble accuracy to 0, so that if partial reassembly
- // of the jacobian matrix is disabled, the reassemble accuracy
- // is always smaller than the current relative tolerance
- reassembleAccuracy_ = 0.0;
- }
-
- ~BoxAssembler()
- {
- delete matrix_;
- }
-
- /*!
- * \brief Initialize the jacobian assembler.
- *
- * At this point we can assume that all objects in the problem and
- * the model have been allocated. We can not assume that they are
- * fully initialized, though.
- *
- * \param problem The problem object
- */
- void init(Problem& problem)
- {
- problemPtr_ = &problem;
-
- // initialize the BCRS matrix
- createMatrix_();
-
- // initialize the jacobian matrix and the right hand side
- // vector
- *matrix_ = 0;
- reuseMatrix_ = false;
-
- int numVerts = gridView_().size(dim);
- int numElems = gridView_().size(0);
-
- residual_.resize(numVerts);
-
- // initialize the storage part of the Jacobian matrix. Since
- // we only need this if Jacobian matrix recycling is enabled,
- // we do not waste space if it is disabled
- if (enableJacobianRecycling_()) {
- storageJacobian_.resize(numVerts);
- storageTerm_.resize(numVerts);
- }
-
- if (gridView_().comm().size() > 1)
- totalElems_ = gridView_().comm().sum(numElems);
- else
- totalElems_ = numElems;
-
- // initialize data needed for partial reassembly
- if (enablePartialReassemble_()) {
- vertexColor_.resize(numVerts);
- vertexDelta_.resize(numVerts);
- elementColor_.resize(numElems);
- }
- reassembleAll();
- }
-
- /*!
- * \brief Assemble the global Jacobian of the residual and the residual for the current solution.
- *
- * The current state of affairs (esp. the previous and the current
- * solutions) is represented by the model object.
- */
- void assemble()
- {
- bool printReassembleStatistics = enablePartialReassemble_() && !reuseMatrix_;
- int succeeded;
- try {
- assemble_();
- succeeded = 1;
- if (gridView_().comm().size() > 1)
- succeeded = gridView_().comm().min(succeeded);
- }
- catch (Dumux::NumericalProblem &e)
- {
- std::cout << "rank " << problem_().gridView().comm().rank()
- << " caught an exception while assembling:" << e.what()
- << "\n";
- succeeded = 0;
- if (gridView_().comm().size() > 1)
- succeeded = gridView_().comm().min(succeeded);
- }
-
- if (!succeeded) {
- DUNE_THROW(NumericalProblem,
- "A process did not succeed in linearizing the system");
- }
-
- if (printReassembleStatistics)
- {
- if (gridView_().comm().size() > 1)
- {
- greenElems_ = gridView_().comm().sum(greenElems_);
- reassembleAccuracy_ = gridView_().comm().max(nextReassembleAccuracy_);
- }
- else
- {
- reassembleAccuracy_ = nextReassembleAccuracy_;
- }
-
- problem_().newtonController().endIterMsg()
- << ", reassembled "
- << totalElems_ - greenElems_ << "/" << totalElems_
- << " (" << 100*Scalar(totalElems_ - greenElems_)/totalElems_ << "%) elems @accuracy="
- << reassembleAccuracy_;
- }
-
- // reset all vertex colors to green
- for (unsigned int i = 0; i < vertexColor_.size(); ++i) {
- vertexColor_[i] = Green;
- }
- }
-
- /*!
- * \brief If Jacobian matrix recycling is enabled, this method
- * specifies whether the next call to assemble() just
- * rescales the storage term or does a full reassembly
- *
- * \param yesno If true, only rescale; else do full Jacobian assembly.
- */
- void setMatrixReuseable(const bool yesno = true)
- {
- if (enableJacobianRecycling_())
- reuseMatrix_ = yesno;
- }
-
- /*!
- * \brief If partial Jacobian matrix reassembly is enabled, this
- * method causes all elements to be reassembled in the next
- * assemble() call.
- */
- void reassembleAll()
- {
- // do not reuse the current linearization
- reuseMatrix_ = false;
-
- // do not use partial reassembly for the next iteration
- nextReassembleAccuracy_ = 0.0;
- if (enablePartialReassemble_()) {
- std::fill(vertexColor_.begin(),
- vertexColor_.end(),
- Red);
- std::fill(elementColor_.begin(),
- elementColor_.end(),
- Red);
- std::fill(vertexDelta_.begin(),
- vertexDelta_.end(),
- 0.0);
- }
- }
-
- /*!
- * \brief Returns the largest relative error of a "green" vertex
- * for the most recent call of the assemble() method.
- *
- * This only has an effect if partial Jacobian reassembly is
- * enabled. If it is disabled, then this method always returns 0.
- *
- * This returns the _actual_ relative computed seen by
- * computeColors(), not the tolerance which it was given.
- */
- Scalar reassembleAccuracy() const
- { return reassembleAccuracy_; }
-
- /*!
- * \brief Update the distance where the non-linear system was
- * originally insistently linearized and the point where it
- * will be linerized the next time.
- *
- * This only has an effect if partial reassemble is enabled.
- */
- void updateDiscrepancy(const SolutionVector &u,
- const SolutionVector &uDelta)
- {
- if (!enablePartialReassemble_())
- return;
-
- // update the vector with the distances of the current
- // evaluation point used for linearization from the original
- // evaluation point
- for (unsigned int i = 0; i < vertexDelta_.size(); ++i) {
- PrimaryVariables currentPriVars(u[i]);
- PrimaryVariables nextPriVars(currentPriVars);
- nextPriVars -= uDelta[i];
-
- // we need to add the distance the solution was moved for
- // this vertex
- Scalar dist = model_().relativeErrorVertex(i,
- currentPriVars,
- nextPriVars);
- vertexDelta_[i] += std::abs(dist);
- }
-
- }
-
- /*!
- * \brief Force to reassemble a given vertex next time the
- * assemble() method is called.
- *
- * \param globalVertIdx The global index of the vertex which ought
- * to be red.
- */
- void markVertexRed(const int globalVertIdx)
- {
- if (!enablePartialReassemble_())
- return;
-
- vertexColor_[globalVertIdx] = Red;
- }
-
/*!
* \brief Determine the colors of vertices and elements for partial
* reassembly given a relative tolerance.
@@ -334,25 +74,25 @@ public:
* linearization point and the current solution and
* _not_ the delta vector of the Newton iteration!
*/
- void computeColors(const Scalar relTol)
+ void computeColors_(const Scalar relTol)
{
- if (!enablePartialReassemble_())
+ if (!this->enablePartialReassemble_())
return;
- ElementIterator elemIt = gridView_().template begin<0>();
- ElementIterator elemEndIt = gridView_().template end<0>();
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator elemEndIt = this->gridView_().template end<0>();
// mark the red vertices and update the tolerance of the
// linearization which actually will get achieved
- nextReassembleAccuracy_ = 0;
- for (unsigned int i = 0; i < vertexColor_.size(); ++i) {
- if (vertexDelta_[i] > relTol)
+ this->nextReassembleAccuracy_ = 0;
+ for (unsigned int i = 0; i < this->vertexColor_.size(); ++i) {
+ if (this->delta_[i] > relTol)
// mark vertex as red if discrepancy is larger than
// the relative tolerance
- vertexColor_[i] = Red;
+ this->vertexColor_[i] = ParentType::Red;
else
- nextReassembleAccuracy_ =
- std::max(nextReassembleAccuracy_, vertexDelta_[i]);
+ this->nextReassembleAccuracy_ =
+ std::max(this->nextReassembleAccuracy_, this->delta_[i]);
}
// Mark all red elements
@@ -362,8 +102,8 @@ public:
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) {
+ int globalI = this->vertexMapper_().map(*elemIt, i, dim);
+ if (this->vertexColor_[globalI] == ParentType::Red) {
isRed = true;
break;
}
@@ -371,28 +111,28 @@ public:
// 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);
+ int globalElemIdx = this->elementMapper_().map(*elemIt);
if (isRed)
- elementColor_[globalElemIdx] = Red;
+ this->elementColor_[globalElemIdx] = ParentType::Red;
else
- elementColor_[globalElemIdx] = Green;
+ this->elementColor_[globalElemIdx] = ParentType::Green;
}
// Mark yellow vertices (as orange for the mean time)
- elemIt = gridView_().template begin<0>();
+ elemIt = this->gridView_().template begin<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = this->elementMapper_().map(*elemIt);
- if (elementColor_[elemIdx] != Red)
+ if (this->elementColor_[elemIdx] != ParentType::Red)
continue; // non-red elements do not tint vertices
// yellow!
int numVerts = elemIt->template count<dim>();
- for (int i=0; i < numVerts; ++i) {
- int globalI = vertexMapper_().map(*elemIt, i, dim);
+ for (int i = 0; i < numVerts; ++i) {
+ int globalI = this->vertexMapper_().map(*elemIt, i, dim);
// if a vertex is already red, don't recolor it to
// yellow!
- if (vertexColor_[globalI] != Red) {
- vertexColor_[globalI] = Orange;
+ if (this->vertexColor_[globalI] != ParentType::Red) {
+ this->vertexColor_[globalI] = ParentType::Orange;
}
}
}
@@ -400,17 +140,18 @@ public:
// at this point we communicate the yellow vertices to the
// neighboring processes because a neigbor process may not see
// the red vertex for yellow border vertices
+ typedef typename ParentType::EntityColor EntityColor;
VertexHandleMin<EntityColor, std::vector<EntityColor>, VertexMapper>
- minHandle(vertexColor_, vertexMapper_());
- gridView_().communicate(minHandle,
+ minHandle(this->vertexColor_, this->vertexMapper_());
+ this->gridView_().communicate(minHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
// Mark yellow elements
- elemIt = gridView_().template begin<0>();
+ elemIt = this->gridView_().template begin<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = this->elementMapper_().map(*elemIt);
- if (elementColor_[elemIdx] == Red) {
+ if (this->elementColor_[elemIdx] == ParentType::Red) {
continue; // element is red already!
}
@@ -418,149 +159,74 @@ public:
// (resp. orange at this point) vertex
bool isYellow = false;
int numVerts = elemIt->template count<dim>();
- for (int i=0; i < numVerts; ++i) {
- int globalI = vertexMapper_().map(*elemIt, i, dim);
- if (vertexColor_[globalI] == Orange) {
+ for (int i = 0; i < numVerts; ++i) {
+ int globalI = this->vertexMapper_().map(*elemIt, i, dim);
+ if (this->vertexColor_[globalI] == ParentType::Orange) {
isYellow = true;
break;
}
}
if (isYellow)
- elementColor_[elemIdx] = Yellow;
+ this->elementColor_[elemIdx] = ParentType::Yellow;
}
// Demote orange vertices to yellow ones if it has at least
// one green element as a neighbor.
- elemIt = gridView_().template begin<0>();
+ elemIt = this->gridView_().template begin<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = this->elementMapper_().map(*elemIt);
- if (elementColor_[elemIdx] != Green)
+ if (this->elementColor_[elemIdx] != ParentType::Green)
continue; // yellow and red elements do not make
// orange vertices yellow!
int numVerts = elemIt->template count<dim>();
- for (int i=0; i < numVerts; ++i) {
- int globalI = vertexMapper_().map(*elemIt, i, dim);
+ for (int i = 0; i < numVerts; ++i) {
+ int globalI = this->vertexMapper_().map(*elemIt, i, dim);
// if a vertex is orange, recolor it to yellow!
- if (vertexColor_[globalI] == Orange)
- vertexColor_[globalI] = Yellow;
+ if (this->vertexColor_[globalI] == ParentType::Orange)
+ this->vertexColor_[globalI] = ParentType::Yellow;
}
}
// demote the border orange vertices
VertexHandleMax<EntityColor, std::vector<EntityColor>, VertexMapper>
- maxHandle(vertexColor_,
- vertexMapper_());
- gridView_().communicate(maxHandle,
+ maxHandle(this->vertexColor_,
+ this->vertexMapper_());
+ this->gridView_().communicate(maxHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
// promote the remaining orange vertices to red
- for (unsigned int i=0; i < vertexColor_.size(); ++i) {
+ for (unsigned int i=0; i < this->vertexColor_.size(); ++i) {
// if a vertex is green or yellow don't do anything!
- if (vertexColor_[i] == Green || vertexColor_[i] == Yellow)
+ if (this->vertexColor_[i] == ParentType::Green || this->vertexColor_[i] == ParentType::Yellow)
continue;
// make sure the vertex is red (this is a no-op vertices
// which are already red!)
- vertexColor_[i] = Red;
+ this->vertexColor_[i] = ParentType::Red;
// set the error of this vertex to 0 because the system
// will be consistently linearized at this vertex
- vertexDelta_[i] = 0.0;
+ this->delta_[i] = 0.0;
}
}
- /*!
- * \brief Returns the reassemble color of a vertex
- *
- * \param element An element which contains the vertex
- * \param vertIdx The local index of the vertex in the element.
- */
- int vertexColor(const Element &element, const int vertIdx) const
- {
- if (!enablePartialReassemble_())
- return Red; // reassemble unconditionally!
-
- int globalIdx = vertexMapper_().map(element, vertIdx, dim);
- return vertexColor_[globalIdx];
- }
-
- /*!
- * \brief Returns the reassemble color of a vertex
- *
- * \param globalVertIdx The global index of the vertex.
- */
- int vertexColor(const int globalVertIdx) const
- {
- if (!enablePartialReassemble_())
- return Red; // reassemble unconditionally!
- return vertexColor_[globalVertIdx];
- }
-
- /*!
- * \brief Returns the Jacobian reassemble color of an element
- *
- * \param element The Codim-0 DUNE entity
- */
- int elementColor(const Element &element) const
- {
- if (!enablePartialReassemble_())
- return Red; // reassemble unconditionally!
-
- int globalIdx = elementMapper_().map(element);
- return elementColor_[globalIdx];
- }
-
- /*!
- * \brief Returns the Jacobian reassemble color of an element
- *
- * \param globalElementIdx The global index of the element.
- */
- int elementColor(const int globalElementIdx) const
- {
- if (!enablePartialReassemble_())
- return Red; // reassemble unconditionally!
- return elementColor_[globalElementIdx];
- }
-
- /*!
- * \brief Return constant reference to global Jacobian matrix.
- */
- const JacobianMatrix& matrix() const
- { return *matrix_; }
- JacobianMatrix& matrix()
- { return *matrix_; }
-
- /*!
- * \brief Return constant reference to global residual vector.
- */
- const SolutionVector& residual() const
- { return residual_; }
- SolutionVector& residual()
- { return residual_; }
-
-private:
- static bool enableJacobianRecycling_()
- { return GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableJacobianRecycling); }
- static bool enablePartialReassemble_()
- { return GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnablePartialReassemble); }
-
// Construct the BCRS matrix for the global jacobian
void createMatrix_()
{
- int numVerticesGlobal = gridView_().size(dim);
+ int numVerticesGlobal = this->gridView_().size(dim);
// allocate raw matrix
- matrix_ = new JacobianMatrix(numVerticesGlobal, numVerticesGlobal, JacobianMatrix::random);
+ this->matrix_ = new JacobianMatrix(numVerticesGlobal, numVerticesGlobal, JacobianMatrix::random);
// find out the global indices of the neighboring vertices of
// each vertex
typedef std::set<int> NeighborSet;
std::vector<NeighborSet> neighbors(numVerticesGlobal);
- ElementIterator eIt = gridView_().template begin<0>();
- const ElementIterator eEndIt = gridView_().template end<0>();
+ ElementIterator eIt = this->gridView_().template begin<0>();
+ const ElementIterator eEndIt = this->gridView_().template end<0>();
for (; eIt != eEndIt; ++eIt) {
const Element &elem = *eIt;
int numVerticesLocal = elem.template count<dim>();
@@ -571,7 +237,7 @@ private:
elem.partitionType() != Dune::BorderEntity)
{
for (int i = 0; i < numVerticesLocal; ++i) {
- int globalI = vertexMapper_().map(*eIt, i, dim);
+ int globalI = this->vertexMapper_().map(*eIt, i, dim);
neighbors[globalI].insert(globalI);
}
}
@@ -579,9 +245,9 @@ private:
{
// loop over all element vertices
for (int i = 0; i < numVerticesLocal - 1; ++i) {
- int globalI = vertexMapper_().map(*eIt, i, dim);
+ int globalI = this->vertexMapper_().map(*eIt, i, dim);
for (int j = i + 1; j < numVerticesLocal; ++j) {
- int globalJ = vertexMapper_().map(*eIt, j, dim);
+ int globalJ = this->vertexMapper_().map(*eIt, j, dim);
// make sure that vertex j is in the neighbor set
// of vertex i and vice-versa
neighbors[globalI].insert(globalJ);
@@ -598,9 +264,9 @@ private:
// allocate space for the rows of the matrix
for (int i = 0; i < numVerticesGlobal; ++i) {
- matrix_->setrowsize(i, neighbors[i].size());
+ this->matrix_->setrowsize(i, neighbors[i].size());
}
- matrix_->endrowsizes();
+ this->matrix_->endrowsizes();
// fill the rows with indices. each vertex talks to all of its
// neighbors. (it also talks to itself since vertices are
@@ -609,158 +275,51 @@ private:
typename NeighborSet::iterator nIt = neighbors[i].begin();
typename NeighborSet::iterator nEndIt = neighbors[i].end();
for (; nIt != nEndIt; ++nIt) {
- matrix_->addindex(i, *nIt);
- }
- }
- matrix_->endindices();
- }
-
- // reset the global linear system of equations. if partial
- // reassemble is enabled, this means that the jacobian matrix must
- // only be erased partially!
- void resetSystem_()
- {
- // do not do anything if we can re-use the current linearization
- if (reuseMatrix_)
- return;
-
- // reset the right hand side.
- residual_ = 0.0;
-
- if (!enablePartialReassemble_()) {
- // If partial reassembly of the jacobian is not enabled,
- // we can just reset everything!
- (*matrix_) = 0;
-
- // reset the parts needed for Jacobian recycling
- if (enableJacobianRecycling_()) {
- int numVerticesGlobal = matrix_->N();
- for (int i=0; i < numVerticesGlobal; ++ i) {
- storageJacobian_[i] = 0;
- storageTerm_[i] = 0;
- }
- }
-
- return;
- }
-
- // reset all entries corrosponding to a red or yellow vertex
- for (unsigned int rowIdx = 0; rowIdx < matrix_->N(); ++rowIdx) {
- if (vertexColor_[rowIdx] == Green)
- continue; // the equations for this control volume are
- // already below the treshold
-
- // here we have yellow or red vertices...
-
- // reset the parts needed for Jacobian recycling
- if (enableJacobianRecycling_()) {
- storageJacobian_[rowIdx] = 0;
- storageTerm_[rowIdx] = 0;
- }
-
- // set all matrix entries in the row to 0
- typedef typename JacobianMatrix::ColIterator ColIterator;
- ColIterator colIt = (*matrix_)[rowIdx].begin();
- const ColIterator &colEndIt = (*matrix_)[rowIdx].end();
- for (; colIt != colEndIt; ++colIt) {
- (*colIt) = 0.0;
- }
- }
- }
-
- // linearize the whole system
- void assemble_()
- {
- resetSystem_();
-
- // if we can "recycle" the current linearization, we do it
- // here and be done with it...
- Scalar curDt = problem_().timeManager().timeStepSize();
- if (reuseMatrix_) {
- int numVerticesGlobal = storageJacobian_.size();
- for (int i = 0; i < numVerticesGlobal; ++i) {
- // rescale the mass term of the jacobian matrix
- MatrixBlock &J_i_i = (*matrix_)[i][i];
-
- J_i_i -= storageJacobian_[i];
- storageJacobian_[i] *= oldDt_/curDt;
- J_i_i += storageJacobian_[i];
-
- // use the flux term plus the source term as the new
- // residual (since the delta in the d(storage)/dt is 0
- // for the first iteration and the residual is
- // approximately 0 in the last iteration, the flux
- // term plus the source term must be equal to the
- // negative change of the storage term of the last
- // iteration of the last time step...)
- residual_[i] = storageTerm_[i];
- residual_[i] *= -1;
- }
-
- reuseMatrix_ = false;
- oldDt_ = curDt;
- return;
- }
-
- oldDt_ = curDt;
- greenElems_ = 0;
-
- // reassemble the elements...
- ElementIterator elemIt = gridView_().template begin<0>();
- ElementIterator elemEndIt = gridView_().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- const Element &elem = *elemIt;
- if (elem.partitionType() != Dune::InteriorEntity &&
- elem.partitionType() != Dune::BorderEntity)
- {
- assembleGhostElement_(elem);
- }
- else
- {
- assembleElement_(elem);
+ this->matrix_->addindex(i, *nIt);
}
}
+ this->matrix_->endindices();
}
// assemble a non-ghost element
void assembleElement_(const Element &elem)
{
- if (enablePartialReassemble_()) {
- int globalElemIdx = model_().elementMapper().map(elem);
- if (elementColor_[globalElemIdx] == Green) {
- ++greenElems_;
+ if (this->enablePartialReassemble_()) {
+ int globalElemIdx = this->model_().elementMapper().map(elem);
+ if (this->elementColor_[globalElemIdx] == ParentType::Green) {
+ ++this->greenElems_;
assembleGreenElement_(elem);
return;
}
}
- model_().localJacobian().assemble(elem);
+ this->model_().localJacobian().assemble(elem);
int numVerticesLocal = elem.template count<dim>();
for (int i=0; i < numVerticesLocal; ++ i) {
- int globI = vertexMapper_().map(elem, i, dim);
+ int globI = this->vertexMapper_().map(elem, i, dim);
// update the right hand side
- residual_[globI] += model_().localJacobian().residual(i);
- for (int j = 0; j < residual_[globI].dimension; ++j)
- assert(std::isfinite(residual_[globI][j]));
- if (enableJacobianRecycling_()) {
- storageTerm_[globI] +=
- model_().localJacobian().storageTerm(i);
+ this->residual_[globI] += this->model_().localJacobian().residual(i);
+ for (int j = 0; j < this->residual_[globI].dimension; ++j)
+ assert(std::isfinite(this->residual_[globI][j]));
+ if (this->enableJacobianRecycling_()) {
+ this->storageTerm_[globI] +=
+ this->model_().localJacobian().storageTerm(i);
}
// only update the jacobian matrix for non-green vertices
- if (vertexColor(globI) != Green) {
- if (enableJacobianRecycling_())
- storageJacobian_[globI] +=
- model_().localJacobian().storageJacobian(i);
+ if (this->vertexColor(globI) != ParentType::Green) {
+ if (this->enableJacobianRecycling_())
+ this->storageJacobian_[globI] +=
+ this->model_().localJacobian().storageJacobian(i);
// update the jacobian matrix
- for (int j=0; j < numVerticesLocal; ++ j) {
- int globJ = vertexMapper_().map(elem, j, dim);
- (*matrix_)[globI][globJ] +=
- model_().localJacobian().mat(i,j);
+ for (int j = 0; j < numVerticesLocal; ++ j) {
+ int globJ = this->vertexMapper_().map(elem, j, dim);
+ (*this->matrix_)[globI][globJ] +=
+ this->model_().localJacobian().mat(i,j);
}
}
}
@@ -770,16 +329,16 @@ private:
// residual updated, but the jacobian is left alone...
void assembleGreenElement_(const Element &elem)
{
- model_().localResidual().eval(elem);
+ this->model_().localResidual().eval(elem);
int numVerticesLocal = elem.template count<dim>();
- for (int i=0; i < numVerticesLocal; ++ i) {
- int globI = vertexMapper_().map(elem, i, dim);
+ for (int i = 0; i < numVerticesLocal; ++ i) {
+ int globI = this->vertexMapper_().map(elem, i, dim);
// update the right hand side
- residual_[globI] += model_().localResidual().residual(i);
- if (enableJacobianRecycling_())
- storageTerm_[globI] += model_().localResidual().storageTerm(i);
+ this->residual_[globI] += this->model_().localResidual().residual(i);
+ if (this->enableJacobianRecycling_())
+ this->storageTerm_[globI] += this->model_().localResidual().storageTerm(i);
}
}
@@ -798,59 +357,16 @@ private:
}
// set main diagonal entries for the vertex
- int vIdx = vertexMapper_().map(*vp);
+ int vIdx = this->vertexMapper_().map(*vp);
typedef typename JacobianMatrix::block_type BlockType;
- BlockType &J = (*matrix_)[vIdx][vIdx];
+ BlockType &J = (*this->matrix_)[vIdx][vIdx];
for (int j = 0; j < BlockType::rows; ++j)
J[j][j] = 1.0;
// set residual for the vertex
- residual_[vIdx] = 0;
+ this->residual_[vIdx] = 0;
}
}
-
-
- Problem &problem_()
- { return *problemPtr_; }
- const Problem &problem_() const
- { return *problemPtr_; }
- const Model &model_() const
- { return problem_().model(); }
- Model &model_()
- { return problem_().model(); }
- const GridView &gridView_() const
- { return problem_().gridView(); }
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); }
- const ElementMapper &elementMapper_() const
- { return problem_().elementMapper(); }
-
- Problem *problemPtr_;
-
- // the jacobian matrix
- JacobianMatrix *matrix_;
- // the right-hand side
- SolutionVector residual_;
-
- // attributes required for jacobian matrix recycling
- bool reuseMatrix_;
- // The storage part of the local Jacobian
- std::vector<MatrixBlock> storageJacobian_;
- std::vector<VectorBlock> storageTerm_;
- // time step size of last assembly
- Scalar oldDt_;
-
-
- // attributes required for partial jacobian reassembly
- std::vector<EntityColor> vertexColor_;
- std::vector<EntityColor> elementColor_;
- std::vector<Scalar> vertexDelta_;
-
- int totalElems_;
- int greenElems_;
-
- Scalar nextReassembleAccuracy_;
- Scalar reassembleAccuracy_;
};
} // namespace Dumux
diff --git a/dumux/implicit/cellcentered/ccassembler.hh b/dumux/implicit/cellcentered/ccassembler.hh
index 8969c898fb..f4285df634 100644
--- a/dumux/implicit/cellcentered/ccassembler.hh
+++ b/dumux/implicit/cellcentered/ccassembler.hh
@@ -23,261 +23,39 @@
/*!
* \file
*
- * \brief An assembler for the global Jacobian matrix for models using the box discretization.
+ * \brief An assembler for the global Jacobian matrix for models using the cell centered discretization.
*/
#ifndef DUMUX_CC_ASSEMBLER_HH
#define DUMUX_CC_ASSEMBLER_HH
-#include <dune/grid/common/gridenums.hh>
-
-#include <dumux/implicit/cellcentered/ccproperties.hh>
-#include <dumux/linear/vertexborderlistfromgrid.hh>
-#include <dumux/linear/foreignoverlapfrombcrsmatrix.hh>
-#include <dumux/parallel/vertexhandles.hh>
+#include <dumux/implicit/common/implicitassembler.hh>
namespace Dumux {
/*!
- * \brief An assembler for the global Jacobian matrix for models using the box discretization.
+ * \brief An assembler for the global Jacobian matrix for models using the cell centered discretization.
*/
template<class TypeTag>
-class CCAssembler
+class CCAssembler : public ImplicitAssembler<TypeTag>
{
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef ImplicitAssembler<TypeTag> ParentType;
+ friend class ImplicitAssembler<TypeTag>;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
-
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-
- enum{ dim = GridView::dimension };
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
+ enum{ dim = GridView::dimension };
typedef typename GridView::IntersectionIterator IntersectionIterator;
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
- typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
- typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
+public:
+ CCAssembler(): ParentType() {}
+private:
// copying the jacobian assembler is not a good idea
CCAssembler(const CCAssembler &);
-public:
- /*!
- * \brief The colors of elements required for partial
- * Jacobian reassembly.
- */
- enum EntityColor {
- /*!
- * Element that needs to be reassembled because some
- * relative error is above the tolerance
- */
- Red = 0,
-
- /*!
- * Element that does not need to be reassembled
- */
- Green = 3
- };
-
- CCAssembler()
- {
- problemPtr_ = 0;
- matrix_ = 0;
-
- // set reassemble accuracy to 0, so that if partial reassembly
- // of the jacobian matrix is disabled, the reassemble accuracy
- // is always smaller than the current relative tolerance
- reassembleAccuracy_ = 0.0;
- }
-
- ~CCAssembler()
- {
- delete matrix_;
- }
-
- /*!
- * \brief Initialize the jacobian assembler.
- *
- * At this point we can assume that all objects in the problem and
- * the model have been allocated. We can not assume that they are
- * fully initialized, though.
- *
- * \param problem The problem object
- */
- void init(Problem& problem)
- {
- problemPtr_ = &problem;
-
- // initialize the BCRS matrix
- createMatrix_();
-
- // initialize the jacobian matrix and the right hand side
- // vector
- *matrix_ = 0;
- reuseMatrix_ = false;
-
- int numElems = gridView_().size(0);
-
- residual_.resize(numElems);
-
- // initialize the storage part of the Jacobian matrix. Since
- // we only need this if Jacobian matrix recycling is enabled,
- // we do not waste space if it is disabled
- if (enableJacobianRecycling_()) {
- storageJacobian_.resize(numElems);
- storageTerm_.resize(numElems);
- }
-
- if (gridView_().comm().size() > 1)
- totalElems_ = gridView_().comm().sum(numElems);
- else
- totalElems_ = numElems;
-
- // initialize data needed for partial reassembly
- if (enablePartialReassemble_())
- {
- elementColor_.resize(numElems);
- elementDelta_.resize(numElems);
- }
-
- reassembleAll();
- }
-
- /*!
- * \brief Assemble the global Jacobian of the residual and the residual for the current solution.
- *
- * The current state of affairs (esp. the previous and the current
- * solutions) is represented by the model object.
- */
- void assemble()
- {
- bool printReassembleStatistics = enablePartialReassemble_() && !reuseMatrix_;
- int succeeded;
- try {
- assemble_();
- succeeded = 1;
- if (gridView_().comm().size() > 1)
- succeeded = gridView_().comm().min(succeeded);
- }
- catch (Dumux::NumericalProblem &e)
- {
- std::cout << "rank " << problem_().gridView().comm().rank()
- << " caught an exception while assembling:" << e.what()
- << "\n";
- succeeded = 0;
- if (gridView_().comm().size() > 1)
- succeeded = gridView_().comm().min(succeeded);
- }
-
- if (!succeeded) {
- DUNE_THROW(NumericalProblem,
- "A process did not succeed in linearizing the system");
- }
-
- if (printReassembleStatistics)
- {
- if (gridView_().comm().size() > 1)
- {
- greenElems_ = gridView_().comm().sum(greenElems_);
- reassembleAccuracy_ = gridView_().comm().max(nextReassembleAccuracy_);
- }
- else
- {
- reassembleAccuracy_ = nextReassembleAccuracy_;
- }
-
- problem_().newtonController().endIterMsg()
- << ", reassembled "
- << totalElems_ - greenElems_ << "/" << totalElems_
- << " (" << 100*Scalar(totalElems_ - greenElems_)/totalElems_ << "%) elems @accuracy="
- << reassembleAccuracy_;
- }
- }
-
- /*!
- * \brief If Jacobian matrix recycling is enabled, this method
- * specifies whether the next call to assemble() just
- * rescales the storage term or does a full reassembly
- *
- * \param yesno If true, only rescale; else do full Jacobian assembly.
- */
- void setMatrixReuseable(bool yesno = true)
- {
- if (enableJacobianRecycling_())
- reuseMatrix_ = yesno;
- }
-
- /*!
- * \brief If partial Jacobian matrix reassembly is enabled, this
- * method causes all elements to be reassembled in the next
- * assemble() call.
- */
- void reassembleAll()
- {
- // do not reuse the current linearization
- reuseMatrix_ = false;
-
- // do not use partial reassembly for the next iteration
- nextReassembleAccuracy_ = 0.0;
- if (enablePartialReassemble_()) {
- std::fill(elementColor_.begin(),
- elementColor_.end(),
- Red);
- std::fill(elementDelta_.begin(),
- elementDelta_.end(),
- 0.0);
- }
- }
-
- /*!
- * \brief Returns the largest relative error of a "green" vertex
- * for the most recent call of the assemble() method.
- *
- * This only has an effect if partial Jacobian reassembly is
- * enabled. If it is disabled, then this method always returns 0.
- *
- * This returns the _actual_ relative computed seen by
- * computeColors(), not the tolerance which it was given.
- */
- Scalar reassembleAccuracy() const
- { return reassembleAccuracy_; }
-
- /*!
- * \brief Update the distance where the non-linear system was
- * originally insistently linearized and the point where it
- * will be linerized the next time.
- *
- * This only has an effect if partial reassemble is enabled.
- */
- void updateDiscrepancy(const SolutionVector &u,
- const SolutionVector &uDelta)
- {
- if (!enablePartialReassemble_())
- return;
-
- // update the vector with the distances of the current
- // evaluation point used for linearization from the original
- // evaluation point
- for (unsigned int i = 0; i < elementDelta_.size(); ++i) {
- PrimaryVariables currentPriVars(u[i]);
- PrimaryVariables nextPriVars(currentPriVars);
- nextPriVars -= uDelta[i];
-
- // we need to add the distance the solution was moved for
- // this vertex
- Scalar dist = model_().relativeErrorVertex(i,
- currentPriVars,
- nextPriVars);
- elementDelta_[i] += std::abs(dist);
- }
- }
-
/*!
* \brief Determine the colors of elements for partial
* reassembly given a relative tolerance.
@@ -295,130 +73,79 @@ public:
* linearization point and the current solution and
* _not_ the delta vector of the Newton iteration!
*/
- void computeColors(Scalar relTol)
+ void computeColors_(Scalar relTol)
{
- if (!enablePartialReassemble_())
+ if (!this->enablePartialReassemble_())
return;
- ElementIterator elemIt = gridView_().template begin<0>();
- ElementIterator elemEndIt = gridView_().template end<0>();
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator elemEndIt = this->gridView_().template end<0>();
// mark the red elements and update the tolerance of the
// linearization which actually will get achieved
- nextReassembleAccuracy_ = 0;
+ this->nextReassembleAccuracy_ = 0;
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = this->elementMapper_().map(*elemIt);
- if (elementDelta_[elemIdx] > relTol)
+ if (this->delta_[elemIdx] > relTol)
{
// mark element as red if discrepancy is larger than
// the relative tolerance
- elementColor_[elemIdx] = Red;
+ this->elementColor_[elemIdx] = ParentType::Red;
}
else
{
- elementColor_[elemIdx] = Green;
- nextReassembleAccuracy_ =
- std::max(nextReassembleAccuracy_, elementDelta_[elemIdx]);
+ this->elementColor_[elemIdx] = ParentType::Green;
+ this->nextReassembleAccuracy_ =
+ std::max(this->nextReassembleAccuracy_, this->delta_[elemIdx]);
}
}
// mark the neighbors also red
- elemIt = gridView_().template begin<0>();
+ elemIt = this->gridView_().template begin<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = this->elementMapper_().map(*elemIt);
- if (elementColor_[elemIdx] == Red)
+ if (this->elementColor_[elemIdx] == ParentType::Red)
continue; // element is red already!
- if (elementDelta_[elemIdx] > relTol)
+ if (this->delta_[elemIdx] > relTol)
{
// also mark the neighbors
- IntersectionIterator endIsIt = gridView_().iend(*elemIt);
- for (IntersectionIterator isIt = gridView_().ibegin(*elemIt); isIt != endIsIt; ++isIt)
+ IntersectionIterator endIsIt = this->gridView_().iend(*elemIt);
+ for (IntersectionIterator isIt = this->gridView_().ibegin(*elemIt); isIt != endIsIt; ++isIt)
{
if (isIt->neighbor())
{
int neighborIdx = this->elementMapper_().map(*isIt->outside());
- elementColor_[neighborIdx] = Red;
+ this->elementColor_[neighborIdx] = ParentType::Red;
}
}
}
}
// set the discrepancy of the red elements to zero
- for (int i = 0; i < elementDelta_.size(); i++)
- if (elementColor_[i] == Red)
- elementDelta_[i] = 0;
- }
-
- /*!
- * \brief Returns the Jacobian reassemble color of an element
- *
- * \param element The Codim-0 DUNE entity
- */
- int elementColor(const Element &element) const
- {
- if (!enablePartialReassemble_())
- return Red; // reassemble unconditionally!
-
- int globalIdx = elementMapper_().map(element);
- return elementColor_[globalIdx];
- }
-
- /*!
- * \brief Returns the Jacobian reassemble color of an element
- *
- * \param globalElementIdx The global index of the element.
- */
- int elementColor(const int globalElementIdx) const
- {
- if (!enablePartialReassemble_())
- return Red; // reassemble unconditionally!
- return elementColor_[globalElementIdx];
+ for (int i = 0; i < this->delta_.size(); i++)
+ if (this->elementColor_[i] == ParentType::Red)
+ this->delta_[i] = 0;
}
-
- /*!
- * \brief Return constant reference to global Jacobian matrix.
- */
- const JacobianMatrix& matrix() const
- { return *matrix_; }
- JacobianMatrix& matrix()
- { return *matrix_; }
-
- /*!
- * \brief Return constant reference to global residual vector.
- */
- const SolutionVector& residual() const
- { return residual_; }
- SolutionVector& residual()
- { return residual_; }
- // functions needed for interface consistence
- void markVertexRed(const int globalVertIdx) {}
-
-private:
- static bool enableJacobianRecycling_()
- { return GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableJacobianRecycling); }
- static bool enablePartialReassemble_()
- { return GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnablePartialReassemble); }
-
// Construct the BCRS matrix for the global jacobian
void createMatrix_()
{
- int nElems = gridView_().size(0);
+ int nElems = this->gridView_().size(0);
// allocate raw matrix
- matrix_ = new JacobianMatrix(nElems, nElems, JacobianMatrix::random);
+ this->matrix_ = new JacobianMatrix(nElems, nElems, JacobianMatrix::random);
// find out the global indices of the neighboring elements of
// each element
typedef std::set<int> NeighborSet;
std::vector<NeighborSet> neighbors(nElems);
- ElementIterator eIt = gridView_().template begin<0>();
- const ElementIterator eEndIt = gridView_().template end<0>();
+ ElementIterator eIt = this->gridView_().template begin<0>();
+ const ElementIterator eEndIt = this->gridView_().template end<0>();
for (; eIt != eEndIt; ++eIt) {
const Element &elem = *eIt;
- int globalI = elementMapper_().map(elem);
+ int globalI = this->elementMapper_().map(elem);
neighbors[globalI].insert(globalI);
// if the element is ghost,
@@ -427,13 +154,13 @@ private:
// continue;
// loop over all neighbors
- IntersectionIterator isIt = gridView_().ibegin(elem);
- const IntersectionIterator &endIt = gridView_().iend(elem);
+ IntersectionIterator isIt = this->gridView_().ibegin(elem);
+ const IntersectionIterator &endIt = this->gridView_().iend(elem);
for (; isIt != endIt; ++isIt)
{
if (isIt->neighbor())
{
- int globalJ = elementMapper_().map(*(isIt->outside()));
+ int globalJ = this->elementMapper_().map(*(isIt->outside()));
neighbors[globalI].insert(globalJ);
}
}
@@ -441,9 +168,9 @@ private:
// allocate space for the rows of the matrix
for (int i = 0; i < nElems; ++i) {
- matrix_->setrowsize(i, neighbors[i].size());
+ this->matrix_->setrowsize(i, neighbors[i].size());
}
- matrix_->endrowsizes();
+ this->matrix_->endrowsizes();
// fill the rows with indices. each element talks to all of its
// neighbors and itself.
@@ -451,157 +178,53 @@ private:
typename NeighborSet::iterator nIt = neighbors[i].begin();
typename NeighborSet::iterator nEndIt = neighbors[i].end();
for (; nIt != nEndIt; ++nIt) {
- matrix_->addindex(i, *nIt);
+ this->matrix_->addindex(i, *nIt);
}
}
- matrix_->endindices();
+ this->matrix_->endindices();
}
- // reset the global linear system of equations. if partial
- // reassemble is enabled, this means that the jacobian matrix must
- // only be erased partially!
- void resetSystem_()
- {
- // do not do anything if we can re-use the current linearization
- if (reuseMatrix_)
- return;
-
- // reset the right hand side.
- residual_ = 0.0;
-
- if (!enablePartialReassemble_()) {
- // If partial reassembly of the jacobian is not enabled,
- // we can just reset everything!
- (*matrix_) = 0;
-
- // reset the parts needed for Jacobian recycling
- if (enableJacobianRecycling_()) {
- int numElementsGlobal = matrix_->N();
- for (int i=0; i < numElementsGlobal; ++ i) {
- storageJacobian_[i] = 0;
- storageTerm_[i] = 0;
- }
- }
-
- return;
- }
-
- // reset all entries corrosponding to a red or yellow vertex
- for (unsigned int rowIdx = 0; rowIdx < matrix_->N(); ++rowIdx) {
- if (elementColor_[rowIdx] == Green)
- continue; // the equations for this control volume are
- // already below the treshold
-
- // reset the parts needed for Jacobian recycling
- if (enableJacobianRecycling_()) {
- storageJacobian_[rowIdx] = 0;
- storageTerm_[rowIdx] = 0;
- }
-
- // set all matrix entries in the row to 0
- typedef typename JacobianMatrix::ColIterator ColIterator;
- ColIterator colIt = (*matrix_)[rowIdx].begin();
- const ColIterator &colEndIt = (*matrix_)[rowIdx].end();
- for (; colIt != colEndIt; ++colIt) {
- (*colIt) = 0.0;
- }
- }
- }
-
- // linearize the whole system
- void assemble_()
- {
- resetSystem_();
-
- // if we can "recycle" the current linearization, we do it
- // here and be done with it...
- Scalar curDt = problem_().timeManager().timeStepSize();
- if (reuseMatrix_) {
- int numElementsGlobal = storageJacobian_.size();
- for (int i = 0; i < numElementsGlobal; ++i) {
- // rescale the mass term of the jacobian matrix
- MatrixBlock &J_i_i = (*matrix_)[i][i];
-
- J_i_i -= storageJacobian_[i];
- storageJacobian_[i] *= oldDt_/curDt;
- J_i_i += storageJacobian_[i];
-
- // use the flux term plus the source term as the new
- // residual (since the delta in the d(storage)/dt is 0
- // for the first iteration and the residual is
- // approximately 0 in the last iteration, the flux
- // term plus the source term must be equal to the
- // negative change of the storage term of the last
- // iteration of the last time step...)
- residual_[i] = storageTerm_[i];
- residual_[i] *= -1;
- }
-
- reuseMatrix_ = false;
- oldDt_ = curDt;
- return;
- }
-
- oldDt_ = curDt;
- greenElems_ = 0;
-
- // reassemble the elements...
- ElementIterator elemIt = gridView_().template begin<0>();
- ElementIterator elemEndIt = gridView_().template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- const Element &elem = *elemIt;
- if (elem.partitionType() == Dune::GhostEntity)
- {
- assembleGhostElement_(elem);
- }
- else
- {
- assembleElement_(elem);
- }
- }
- }
-
// assemble a non-ghost element
void assembleElement_(const Element &elem)
{
- if (enablePartialReassemble_()) {
- int globalElemIdx = model_().elementMapper().map(elem);
- if (elementColor_[globalElemIdx] == Green) {
- ++greenElems_;
+ if (this->enablePartialReassemble_()) {
+ int globalElemIdx = this->model_().elementMapper().map(elem);
+ if (this->elementColor_[globalElemIdx] == ParentType::Green) {
+ ++this->greenElems_;
assembleGreenElement_(elem);
return;
}
}
- model_().localJacobian().assemble(elem);
+ this->model_().localJacobian().assemble(elem);
- int globalI = elementMapper_().map(elem);
+ int globalI = this->elementMapper_().map(elem);
// update the right hand side
- residual_[globalI] = model_().localJacobian().residual(0);
- for (int j = 0; j < residual_[globalI].dimension; ++j)
- assert(std::isfinite(residual_[globalI][j]));
- if (enableJacobianRecycling_()) {
- storageTerm_[globalI] +=
- model_().localJacobian().storageTerm(0);
+ this->residual_[globalI] = this->model_().localJacobian().residual(0);
+ for (int j = 0; j < this->residual_[globalI].dimension; ++j)
+ assert(std::isfinite(this->residual_[globalI][j]));
+ if (this->enableJacobianRecycling_()) {
+ this->storageTerm_[globalI] +=
+ this->model_().localJacobian().storageTerm(0);
}
- if (enableJacobianRecycling_())
- storageJacobian_[globalI] +=
- model_().localJacobian().storageJacobian(0);
+ if (this->enableJacobianRecycling_())
+ this->storageJacobian_[globalI] +=
+ this->model_().localJacobian().storageJacobian(0);
// update the diagonal entry
- (*matrix_)[globalI][globalI] = model_().localJacobian().mat(0,0);
+ (*this->matrix_)[globalI][globalI] = this->model_().localJacobian().mat(0,0);
- IntersectionIterator isIt = gridView_().ibegin(elem);
- const IntersectionIterator &endIt = gridView_().iend(elem);
+ IntersectionIterator isIt = this->gridView_().ibegin(elem);
+ const IntersectionIterator &endIt = this->gridView_().iend(elem);
for (int j = 0; isIt != endIt; ++isIt)
{
if (isIt->neighbor())
{
- int globalJ = elementMapper_().map(*(isIt->outside()));
- (*matrix_)[globalI][globalJ] = model_().localJacobian().mat(0,++j);
+ int globalJ = this->elementMapper_().map(*(isIt->outside()));
+ (*this->matrix_)[globalI][globalJ] = this->model_().localJacobian().mat(0,++j);
}
}
}
@@ -610,72 +233,30 @@ private:
// residual updated, but the jacobian is left alone...
void assembleGreenElement_(const Element &elem)
{
- model_().localResidual().eval(elem);
+ this->model_().localResidual().eval(elem);
- int globalI = elementMapper_().map(elem);
+ int globalI = this->elementMapper_().map(elem);
// update the right hand side
- residual_[globalI] += model_().localResidual().residual(0);
- if (enableJacobianRecycling_())
- storageTerm_[globalI] += model_().localResidual().storageTerm(0);
+ this->residual_[globalI] += this->model_().localResidual().residual(0);
+ if (this->enableJacobianRecycling_())
+ this->storageTerm_[globalI] += this->model_().localResidual().storageTerm(0);
}
// "assemble" a ghost element
void assembleGhostElement_(const Element &elem)
{
- int globalI = elementMapper_().map(elem);
+ int globalI = this->elementMapper_().map(elem);
// update the right hand side
- residual_[globalI] = 0.0;
+ this->residual_[globalI] = 0.0;
// update the diagonal entry
typedef typename JacobianMatrix::block_type BlockType;
- BlockType &J = (*matrix_)[globalI][globalI];
+ BlockType &J = (*this->matrix_)[globalI][globalI];
for (int j = 0; j < BlockType::rows; ++j)
J[j][j] = 1.0;
}
-
-
- Problem &problem_()
- { return *problemPtr_; }
- const Problem &problem_() const
- { return *problemPtr_; }
- const Model &model_() const
- { return problem_().model(); }
- Model &model_()
- { return problem_().model(); }
- const GridView &gridView_() const
- { return problem_().gridView(); }
- const VertexMapper &vertexMapper_() const
- { return problem_().vertexMapper(); }
- const ElementMapper &elementMapper_() const
- { return problem_().elementMapper(); }
-
- Problem *problemPtr_;
-
- // the jacobian matrix
- JacobianMatrix *matrix_;
- // the right-hand side
- SolutionVector residual_;
-
- // attributes required for jacobian matrix recycling
- bool reuseMatrix_;
- // The storage part of the local Jacobian
- std::vector<MatrixBlock> storageJacobian_;
- std::vector<VectorBlock> storageTerm_;
- // time step size of last assembly
- Scalar oldDt_;
-
-
- // attributes required for partial jacobian reassembly
- std::vector<EntityColor> elementColor_;
- std::vector<Scalar> elementDelta_;
-
- int totalElems_;
- int greenElems_;
-
- Scalar nextReassembleAccuracy_;
- Scalar reassembleAccuracy_;
};
} // namespace Dumux
diff --git a/dumux/implicit/common/implicitassembler.hh b/dumux/implicit/common/implicitassembler.hh
index d585d360a2..5ca3d76c91 100644
--- a/dumux/implicit/common/implicitassembler.hh
+++ b/dumux/implicit/common/implicitassembler.hh
@@ -19,26 +19,22 @@
/*!
* \file
*
- * \brief An assembler for the global Jacobian matrix for models using the box discretization.
+ * \brief An assembler for the global Jacobian matrix for fully implicit models.
*/
-#ifndef DUMUX_BOX_ASSEMBLER_HH
-#define DUMUX_BOX_ASSEMBLER_HH
+#ifndef DUMUX_IMPLICIT_ASSEMBLER_HH
+#define DUMUX_IMPLICIT_ASSEMBLER_HH
-#include <dune/grid/common/gridenums.hh>
-
-#include <dumux/implicit/box/boxproperties.hh>
-#include <dumux/linear/vertexborderlistfromgrid.hh>
-#include <dumux/linear/foreignoverlapfrombcrsmatrix.hh>
-#include <dumux/parallel/vertexhandles.hh>
+#include "implicitproperties.hh"
namespace Dumux {
/*!
- * \brief An assembler for the global Jacobian matrix for models using the box discretization.
+ * \brief An assembler for the global Jacobian matrix for fully implicit models.
*/
template<class TypeTag>
-class BoxAssembler
+class ImplicitAssembler
{
+ typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -59,8 +55,10 @@ class BoxAssembler
typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
// copying the jacobian assembler is not a good idea
- BoxAssembler(const BoxAssembler &);
+ ImplicitAssembler(const ImplicitAssembler &);
public:
/*!
@@ -96,7 +94,7 @@ public:
Green = 3
};
- BoxAssembler()
+ ImplicitAssembler()
{
problemPtr_ = 0;
matrix_ = 0;
@@ -107,7 +105,7 @@ public:
reassembleAccuracy_ = 0.0;
}
- ~BoxAssembler()
+ ~ImplicitAssembler()
{
delete matrix_;
}
@@ -126,7 +124,7 @@ public:
problemPtr_ = &problem;
// initialize the BCRS matrix
- createMatrix_();
+ asImp_().createMatrix_();
// initialize the jacobian matrix and the right hand side
// vector
@@ -135,15 +133,16 @@ public:
int numVerts = gridView_().size(dim);
int numElems = gridView_().size(0);
+ int numDofs = problem.model().numDofs();
- residual_.resize(numVerts);
+ residual_.resize(numDofs);
// initialize the storage part of the Jacobian matrix. Since
// we only need this if Jacobian matrix recycling is enabled,
// we do not waste space if it is disabled
if (enableJacobianRecycling_()) {
- storageJacobian_.resize(numVerts);
- storageTerm_.resize(numVerts);
+ storageJacobian_.resize(numDofs);
+ storageTerm_.resize(numDofs);
}
if (gridView_().comm().size() > 1)
@@ -153,9 +152,10 @@ public:
// initialize data needed for partial reassembly
if (enablePartialReassemble_()) {
- vertexColor_.resize(numVerts);
- vertexDelta_.resize(numVerts);
+ delta_.resize(numDofs);
elementColor_.resize(numElems);
+ if (isBox)
+ vertexColor_.resize(numVerts);
}
reassembleAll();
}
@@ -248,8 +248,8 @@ public:
std::fill(elementColor_.begin(),
elementColor_.end(),
Red);
- std::fill(vertexDelta_.begin(),
- vertexDelta_.end(),
+ std::fill(delta_.begin(),
+ delta_.end(),
0.0);
}
}
@@ -283,7 +283,7 @@ public:
// update the vector with the distances of the current
// evaluation point used for linearization from the original
// evaluation point
- for (unsigned int i = 0; i < vertexDelta_.size(); ++i) {
+ for (unsigned int i = 0; i < delta_.size(); ++i) {
PrimaryVariables currentPriVars(u[i]);
PrimaryVariables nextPriVars(currentPriVars);
nextPriVars -= uDelta[i];
@@ -293,7 +293,7 @@ public:
Scalar dist = model_().relativeErrorVertex(i,
currentPriVars,
nextPriVars);
- vertexDelta_[i] += std::abs(dist);
+ delta_[i] += std::abs(dist);
}
}
@@ -317,17 +317,6 @@ public:
* \brief Determine the colors of vertices and elements for partial
* reassembly given a relative tolerance.
*
- * The following approach is used:
- *
- * - Set all vertices and elements to 'green'
- * - Mark all vertices as 'red' which exhibit an relative error above
- * the tolerance
- * - Mark all elements which feature a 'red' vetex as 'red'
- * - Mark all vertices which are not 'red' and are part of a
- * 'red' element as 'yellow'
- * - Mark all elements which are not 'red' and contain a
- * 'yellow' vertex as 'yellow'
- *
* \param relTol The relative error below which a vertex won't be
* reassembled. Note that this specifies the
* worst-case relative error between the last
@@ -336,140 +325,7 @@ public:
*/
void computeColors(const Scalar relTol)
{
- if (!enablePartialReassemble_())
- return;
-
- ElementIterator elemIt = gridView_().template begin<0>();
- ElementIterator elemEndIt = gridView_().template end<0>();
-
- // mark the red vertices and update the tolerance of the
- // linearization which actually will get achieved
- nextReassembleAccuracy_ = 0;
- for (unsigned int i = 0; i < vertexColor_.size(); ++i) {
- if (vertexDelta_[i] > relTol)
- // mark vertex as red if discrepancy is larger than
- // the relative tolerance
- vertexColor_[i] = Red;
- else
- nextReassembleAccuracy_ =
- std::max(nextReassembleAccuracy_, vertexDelta_[i]);
- }
-
- // Mark all red elements
- for (; elemIt != elemEndIt; ++elemIt) {
- // 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) {
- 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);
- if (isRed)
- elementColor_[globalElemIdx] = Red;
- else
- elementColor_[globalElemIdx] = Green;
- }
-
- // 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] != Red)
- continue; // non-red elements do not tint vertices
- // yellow!
-
- int numVerts = elemIt->template count<dim>();
- for (int i=0; i < numVerts; ++i) {
- int globalI = vertexMapper_().map(*elemIt, i, dim);
- // if a vertex is already red, don't recolor it to
- // yellow!
- if (vertexColor_[globalI] != Red) {
- vertexColor_[globalI] = Orange;
- }
- }
- }
-
- // at this point we communicate the yellow vertices to the
- // neighboring processes because a neigbor process may not see
- // the red vertex for yellow border vertices
- VertexHandleMin<EntityColor, std::vector<EntityColor>, VertexMapper>
- minHandle(vertexColor_, vertexMapper_());
- gridView_().communicate(minHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
-
- // Mark yellow elements
- elemIt = gridView_().template begin<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- int elemIdx = this->elementMapper_().map(*elemIt);
- if (elementColor_[elemIdx] == Red) {
- 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) {
- int globalI = vertexMapper_().map(*elemIt, i, dim);
- if (vertexColor_[globalI] == Orange) {
- isYellow = true;
- break;
- }
- }
-
- if (isYellow)
- elementColor_[elemIdx] = Yellow;
- }
-
- // Demote orange vertices to yellow ones if it has at least
- // one green element as a neighbor.
- elemIt = gridView_().template begin<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- int elemIdx = this->elementMapper_().map(*elemIt);
- if (elementColor_[elemIdx] != Green)
- continue; // yellow and red elements do not make
- // orange vertices yellow!
-
- int numVerts = elemIt->template count<dim>();
- for (int i=0; i < numVerts; ++i) {
- int globalI = vertexMapper_().map(*elemIt, i, dim);
- // if a vertex is orange, recolor it to yellow!
- if (vertexColor_[globalI] == Orange)
- vertexColor_[globalI] = Yellow;
- }
- }
-
- // demote the border orange vertices
- VertexHandleMax<EntityColor, std::vector<EntityColor>, VertexMapper>
- maxHandle(vertexColor_,
- vertexMapper_());
- gridView_().communicate(maxHandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
-
- // promote the remaining orange vertices to red
- for (unsigned int i=0; i < vertexColor_.size(); ++i) {
- // 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 the error of this vertex to 0 because the system
- // will be consistently linearized at this vertex
- vertexDelta_[i] = 0.0;
- }
+ asImp_().computeColors_(relTol);
}
/*!
@@ -541,79 +397,12 @@ public:
SolutionVector& residual()
{ return residual_; }
-private:
+protected:
static bool enableJacobianRecycling_()
{ return GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnableJacobianRecycling); }
static bool enablePartialReassemble_()
{ return GET_PARAM_FROM_GROUP(TypeTag, bool, Implicit, EnablePartialReassemble); }
- // Construct the BCRS matrix for the global jacobian
- void createMatrix_()
- {
- int numVerticesGlobal = gridView_().size(dim);
-
- // allocate raw matrix
- matrix_ = new JacobianMatrix(numVerticesGlobal, numVerticesGlobal, JacobianMatrix::random);
-
- // find out the global indices of the neighboring vertices of
- // each vertex
- typedef std::set<int> NeighborSet;
- std::vector<NeighborSet> neighbors(numVerticesGlobal);
- ElementIterator eIt = gridView_().template begin<0>();
- const ElementIterator eEndIt = gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt) {
- const Element &elem = *eIt;
- int numVerticesLocal = elem.template count<dim>();
-
- // if the element is not in the interior or the process
- // border, all dofs just contain main-diagonal entries
- if (elem.partitionType() != Dune::InteriorEntity &&
- elem.partitionType() != Dune::BorderEntity)
- {
- for (int i = 0; i < numVerticesLocal; ++i) {
- int globalI = vertexMapper_().map(*eIt, i, dim);
- neighbors[globalI].insert(globalI);
- }
- }
- else
- {
- // loop over all element vertices
- for (int i = 0; i < numVerticesLocal - 1; ++i) {
- int globalI = vertexMapper_().map(*eIt, i, dim);
- for (int j = i + 1; j < numVerticesLocal; ++j) {
- int globalJ = vertexMapper_().map(*eIt, j, dim);
- // make sure that vertex j is in the neighbor set
- // of vertex i and vice-versa
- neighbors[globalI].insert(globalJ);
- neighbors[globalJ].insert(globalI);
- }
- }
- }
- }
-
- // make vertices neighbors to themselfs
- for (int i = 0; i < numVerticesGlobal; ++i) {
- neighbors[i].insert(i);
- }
-
- // allocate space for the rows of the matrix
- for (int i = 0; i < numVerticesGlobal; ++i) {
- matrix_->setrowsize(i, neighbors[i].size());
- }
- matrix_->endrowsizes();
-
- // 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 < numVerticesGlobal; ++i) {
- typename NeighborSet::iterator nIt = neighbors[i].begin();
- typename NeighborSet::iterator nEndIt = neighbors[i].end();
- for (; nIt != nEndIt; ++nIt) {
- matrix_->addindex(i, *nIt);
- }
- }
- matrix_->endindices();
- }
// reset the global linear system of equations. if partial
// reassemble is enabled, this means that the jacobian matrix must
@@ -634,8 +423,7 @@ private:
// reset the parts needed for Jacobian recycling
if (enableJacobianRecycling_()) {
- int numVerticesGlobal = matrix_->N();
- for (int i=0; i < numVerticesGlobal; ++ i) {
+ for (int i = 0; i < matrix_->N(); i++) {
storageJacobian_[i] = 0;
storageTerm_[i] = 0;
}
@@ -646,24 +434,22 @@ private:
// reset all entries corrosponding to a red or yellow vertex
for (unsigned int rowIdx = 0; rowIdx < matrix_->N(); ++rowIdx) {
- if (vertexColor_[rowIdx] == Green)
- continue; // the equations for this control volume are
- // already below the treshold
-
- // here we have yellow or red vertices...
-
- // reset the parts needed for Jacobian recycling
- if (enableJacobianRecycling_()) {
- storageJacobian_[rowIdx] = 0;
- storageTerm_[rowIdx] = 0;
- }
+ if ((isBox && vertexColor_[rowIdx] != Green)
+ || (!isBox && elementColor_[rowIdx] != Green))
+ {
+ // reset the parts needed for Jacobian recycling
+ if (enableJacobianRecycling_()) {
+ storageJacobian_[rowIdx] = 0;
+ storageTerm_[rowIdx] = 0;
+ }
- // set all matrix entries in the row to 0
- typedef typename JacobianMatrix::ColIterator ColIterator;
- ColIterator colIt = (*matrix_)[rowIdx].begin();
- const ColIterator &colEndIt = (*matrix_)[rowIdx].end();
- for (; colIt != colEndIt; ++colIt) {
- (*colIt) = 0.0;
+ // set all matrix entries in the row to 0
+ typedef typename JacobianMatrix::ColIterator ColIterator;
+ ColIterator colIt = (*matrix_)[rowIdx].begin();
+ const ColIterator &colEndIt = (*matrix_)[rowIdx].end();
+ for (; colIt != colEndIt; ++colIt) {
+ (*colIt) = 0.0;
+ }
}
}
}
@@ -677,8 +463,7 @@ private:
// here and be done with it...
Scalar curDt = problem_().timeManager().timeStepSize();
if (reuseMatrix_) {
- int numVerticesGlobal = storageJacobian_.size();
- for (int i = 0; i < numVerticesGlobal; ++i) {
+ for (int i = 0; i < matrix_->N(); i++) {
// rescale the mass term of the jacobian matrix
MatrixBlock &J_i_i = (*matrix_)[i][i];
@@ -710,106 +495,18 @@ private:
ElementIterator elemEndIt = gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element &elem = *elemIt;
- if (elem.partitionType() != Dune::InteriorEntity &&
- elem.partitionType() != Dune::BorderEntity)
+ if (elem.partitionType() == Dune::GhostEntity)
{
- assembleGhostElement_(elem);
+ asImp_().assembleGhostElement_(elem);
}
else
{
- assembleElement_(elem);
- }
- }
- }
-
- // assemble a non-ghost element
- void assembleElement_(const Element &elem)
- {
- if (enablePartialReassemble_()) {
- int globalElemIdx = model_().elementMapper().map(elem);
- if (elementColor_[globalElemIdx] == Green) {
- ++greenElems_;
-
- assembleGreenElement_(elem);
- return;
- }
- }
-
- model_().localJacobian().assemble(elem);
-
- int numVerticesLocal = elem.template count<dim>();
- for (int i=0; i < numVerticesLocal; ++ i) {
- int globI = vertexMapper_().map(elem, i, dim);
-
- // update the right hand side
- residual_[globI] += model_().localJacobian().residual(i);
- for (int j = 0; j < residual_[globI].dimension; ++j)
- assert(std::isfinite(residual_[globI][j]));
- if (enableJacobianRecycling_()) {
- storageTerm_[globI] +=
- model_().localJacobian().storageTerm(i);
- }
-
- // only update the jacobian matrix for non-green vertices
- if (vertexColor(globI) != Green) {
- if (enableJacobianRecycling_())
- storageJacobian_[globI] +=
- model_().localJacobian().storageJacobian(i);
-
- // update the jacobian matrix
- for (int j=0; j < numVerticesLocal; ++ j) {
- int globJ = vertexMapper_().map(elem, j, dim);
- (*matrix_)[globI][globJ] +=
- model_().localJacobian().mat(i,j);
- }
- }
- }
- }
-
- // "assemble" a green element. green elements only get the
- // residual updated, but the jacobian is left alone...
- void assembleGreenElement_(const Element &elem)
- {
- model_().localResidual().eval(elem);
-
- int numVerticesLocal = elem.template count<dim>();
- for (int i=0; i < numVerticesLocal; ++ i) {
- int globI = vertexMapper_().map(elem, i, dim);
-
- // update the right hand side
- residual_[globI] += model_().localResidual().residual(i);
- if (enableJacobianRecycling_())
- storageTerm_[globI] += model_().localResidual().storageTerm(i);
- }
- }
-
- // "assemble" a ghost element
- void assembleGhostElement_(const Element &elem)
- {
- int numVerticesLocal = elem.template count<dim>();
- for (int i=0; i < numVerticesLocal; ++i) {
- const VertexPointer vp = elem.template subEntity<dim>(i);
-
- if (vp->partitionType() == Dune::InteriorEntity ||
- vp->partitionType() == Dune::BorderEntity)
- {
- // do not change the non-ghost vertices
- continue;
+ asImp_().assembleElement_(elem);
}
-
- // set main diagonal entries for the vertex
- int vIdx = vertexMapper_().map(*vp);
- typedef typename JacobianMatrix::block_type BlockType;
- BlockType &J = (*matrix_)[vIdx][vIdx];
- for (int j = 0; j < BlockType::rows; ++j)
- J[j][j] = 1.0;
-
- // set residual for the vertex
- residual_[vIdx] = 0;
}
}
-
+protected:
Problem &problem_()
{ return *problemPtr_; }
const Problem &problem_() const
@@ -844,13 +541,19 @@ private:
// attributes required for partial jacobian reassembly
std::vector<EntityColor> vertexColor_;
std::vector<EntityColor> elementColor_;
- std::vector<Scalar> vertexDelta_;
+ std::vector<Scalar> delta_;
int totalElems_;
int greenElems_;
Scalar nextReassembleAccuracy_;
Scalar reassembleAccuracy_;
+
+private:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
};
} // namespace Dumux
--
GitLab