From bf8f55f58fcc8cb07caf287f7ac5230fde1f4b83 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Tue, 10 Aug 2010 08:40:53 +0000
Subject: [PATCH] box: implement partial reassembly of the jacobian matrix
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4041 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/2p2c/2p2cfluidstate.hh | 6 +
.../common/boxelementboundarytypes.hh | 6 +
dumux/boxmodels/common/boxlocaljacobian.hh | 63 ++++-
dumux/boxmodels/common/boxlocalresidual.hh | 3 +
dumux/boxmodels/common/boxmodel.hh | 10 +
dumux/boxmodels/common/boxproblem.hh | 3 +
dumux/boxmodels/common/boxproperties.hh | 8 +
dumux/boxmodels/common/pdelabboxassembler.hh | 234 ++++++++++++++----
.../common/pdelabboxlocaloperator.hh | 28 ++-
dumux/nonlinear/newtoncontroller.hh | 19 +-
dumux/nonlinear/newtonmethod.hh | 1 +
test/boxmodels/2p/test_2p.cc | 4 +-
test/boxmodels/2p2c/injectionproblem.hh | 3 +
13 files changed, 314 insertions(+), 74 deletions(-)
diff --git a/dumux/boxmodels/2p2c/2p2cfluidstate.hh b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
index a965f104b7..07f80c7001 100644
--- a/dumux/boxmodels/2p2c/2p2cfluidstate.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluidstate.hh
@@ -296,6 +296,12 @@ public:
Scalar phasePressure(int phaseIdx) const
{ return phasePressure_[phaseIdx]; }
+ /*!
+ * \brief Return the fugacity of a component [Pa].
+ */
+ Scalar fugacity(int phaseIdx, int compIdx) const
+ { return moleFrac(phaseIdx, compIdx)*phasePressure(compIdx); };
+
/*!
* \brief Returns the capillary pressure [Pa]
*/
diff --git a/dumux/boxmodels/common/boxelementboundarytypes.hh b/dumux/boxmodels/common/boxelementboundarytypes.hh
index 99050933d3..08815ae866 100644
--- a/dumux/boxmodels/common/boxelementboundarytypes.hh
+++ b/dumux/boxmodels/common/boxelementboundarytypes.hh
@@ -51,6 +51,12 @@ class BoxElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTa
enum { dim = GridView::dimension };
public:
+ // copying a the boundary types of an element should be explicitly
+ // requested
+ explicit BoxElementBoundaryTypes(const BoxElementBoundaryTypes &v)
+ : ParentType(v)
+ {}
+
/*!
* \brief The constructor.
*/
diff --git a/dumux/boxmodels/common/boxlocaljacobian.hh b/dumux/boxmodels/common/boxlocaljacobian.hh
index 905b53755f..a0e9775d48 100644
--- a/dumux/boxmodels/common/boxlocaljacobian.hh
+++ b/dumux/boxmodels/common/boxlocaljacobian.hh
@@ -64,12 +64,17 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
enum {
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
+ dimWorld = GridView::dimensionworld,
+
+ Red = JacobianAssembler::Red,
+ Yellow = JacobianAssembler::Yellow,
+ Green = JacobianAssembler::Green
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
@@ -102,10 +107,14 @@ private:
typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
+ // copying a local jacobian is not a good idea
+ BoxLocalJacobian(const BoxLocalJacobian &);
+
public:
BoxLocalJacobian()
{
Valgrind::SetUndefined(problemPtr_);
+ notReassembled = 0;
}
@@ -117,6 +126,7 @@ public:
A_.setSize(2<<dim, 2<<dim);
}
+ int notReassembled;
/*!
* \brief Assemble the linear system of equations for the
* verts of a element, given a local solution 'localU'.
@@ -126,6 +136,15 @@ public:
// set the current grid element and update the element's
// finite volume geometry
elemPtr_ = &element;
+ reset_();
+
+ int elemColor = jacAsm_().elementColor(element);
+ if (elemColor == Green) {
+ ++notReassembled;
+ // Green elements don't need to be reassembled
+ return;
+ }
+
fvElemGeom_.update(gridView_(), elem_());
bcTypes_.update(problem_(), elem_(), fvElemGeom_);
@@ -136,7 +155,6 @@ public:
// not thread save.) The real solution are context objects!
problem_().updateCouplingParams(elem_());
-
int numVertices = fvElemGeom_.numVertices;
// update the secondary variables for the element at the last
@@ -222,12 +240,31 @@ protected:
const Model &model_() const
{ return problem_().model(); };
+ /*!
+ * \brief Returns a reference to the jacobian assembler.
+ */
+ const JacobianAssembler &jacAsm_() const
+ { return model_().jacobianAssembler(); }
+
/*!
* \brief Returns a reference to the vertex mapper.
*/
const VertexMapper &vertexMapper_() const
{ return problem_().vertexMapper(); };
+ /*!
+ * \brief Reset the local jacobian matrix to 0
+ */
+ void reset_()
+ {
+ int n = elem_().template count<dim>();
+ for (int i = 0; i < n; ++ i) {
+ for (int j = 0; j < n; ++ j) {
+ A_[i][j] = 0.0;
+ }
+ }
+ }
+
/*!
* \brief Compute the partial derivatives to a primary variable at
* an degree of freedom.
@@ -318,15 +355,21 @@ protected:
*/
void updateLocalJacobian_(int scvIdx,
int pvIdx,
- const ElementSolutionVector &stiffness)
+ const ElementSolutionVector &deriv)
{
- for (int i = 0; i < fvElemGeom_.numVertices; i++) {
+ for (int i = 0; i < fvElemGeom_.numVertices; i++)
+ {
+ if (jacAsm_().vertexColor(elem_(), i) == Green) {
+ // Green vertices are not to be changed!
+ continue;
+ }
+
for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- // A[i][scvIdx][eqIdx][pvIdx] is the approximate rate
- // of change of the residual of equation 'eqIdx' at
- // vertex 'i' depending on the primary variable
- // 'pvIdx' at vertex 'scvIdx'.
- this->A_[i][scvIdx][eqIdx][pvIdx] = stiffness[i][eqIdx];
+ // A[i][scvIdx][eqIdx][pvIdx] is the rate of change of
+ // the residual of equation 'eqIdx' at vertex 'i'
+ // depending on the primary variable 'pvIdx' at vertex
+ // 'scvIdx'.
+ this->A_[i][scvIdx][eqIdx][pvIdx] = deriv[i][eqIdx];
Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
}
}
@@ -347,7 +390,7 @@ protected:
LocalResidual localResidual_;
LocalBlockMatrix A_;
-};
+ };
}
#endif
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
index 49cd5f9846..e2d26a4021 100644
--- a/dumux/boxmodels/common/boxlocalresidual.hh
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -95,6 +95,9 @@ private:
typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
+ // copying the local residual class is not a good idea
+ BoxLocalResidual(const BoxLocalResidual &);
+
public:
BoxLocalResidual()
{ }
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
index b5a3b3ce6e..4044bcbf9a 100644
--- a/dumux/boxmodels/common/boxmodel.hh
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -96,6 +96,9 @@ class BoxModel
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ // copying a model is not a good idea
+ BoxModel(const BoxModel &);
+
public:
/*!
* \brief The constructor.
@@ -228,6 +231,12 @@ public:
JacobianAssembler &jacobianAssembler()
{ return *jacAsm_; }
+ /*!
+ * \copydoc jacobianAssembler()
+ */
+ const JacobianAssembler &jacobianAssembler() const
+ { return *jacAsm_; }
+
/*!
* \brief Returns the local jacobian which calculates the local
* stiffness matrix for an arbitrary element.
@@ -321,6 +330,7 @@ public:
// previous time step so that we can start the next
// update at a physically meaningful solution.
uCur_ = uPrev_;
+ jacAsm_->reassembleAll();
};
/*!
diff --git a/dumux/boxmodels/common/boxproblem.hh b/dumux/boxmodels/common/boxproblem.hh
index 2abf4e9b17..a6eb89a676 100644
--- a/dumux/boxmodels/common/boxproblem.hh
+++ b/dumux/boxmodels/common/boxproblem.hh
@@ -65,6 +65,9 @@ private:
typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ // copying a problem is not a good idea
+ BoxProblem(const BoxProblem &);
+
public:
BoxProblem(TimeManager &timeManager, const GridView &gridView)
: gridView_(gridView),
diff --git a/dumux/boxmodels/common/boxproperties.hh b/dumux/boxmodels/common/boxproperties.hh
index 4b0d4574aa..f8a6fca422 100644
--- a/dumux/boxmodels/common/boxproperties.hh
+++ b/dumux/boxmodels/common/boxproperties.hh
@@ -92,6 +92,11 @@ NEW_PROP_TAG(GridOperatorSpace); //!< The used grid operator space
//! of the next time step.
NEW_PROP_TAG(EnableJacobianRecycling);
+//! Specify whether the jacobian matrix should be only reassembled for
+//! elements where at least one vertex is above the specified
+//! tolerance
+NEW_PROP_TAG(EnablePartialReassemble);
+
// mappers from local to global indices
NEW_PROP_TAG(VertexMapper);
NEW_PROP_TAG(ElementMapper);
@@ -365,6 +370,9 @@ SET_PROP(BoxModel, LocalOperator)
// enable jacobian matrix recycling by default
SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, true);
+// enable partial reassembling by default
+SET_BOOL_PROP(BoxModel, EnablePartialReassemble, true);
+
// \}
}
diff --git a/dumux/boxmodels/common/pdelabboxassembler.hh b/dumux/boxmodels/common/pdelabboxassembler.hh
index 4f36805ea0..3bec66b0e4 100644
--- a/dumux/boxmodels/common/pdelabboxassembler.hh
+++ b/dumux/boxmodels/common/pdelabboxassembler.hh
@@ -64,7 +64,21 @@ class BoxAssembler
typedef JacobianMatrix Matrix;
typedef Matrix RepresentationType;
+ enum {
+ enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)),
+ enableJacobianRecycling = GET_PROP_VALUE(TypeTag, PTAG(EnableJacobianRecycling))
+ };
+
+ // copying the jacobian assembler is not a good idea
+ BoxAssembler(const BoxAssembler &);
+
public:
+ enum EntityColor {
+ Red, //!< entity needs to be reassembled because it's above the tolerance
+ Yellow, //!< entity needs to be reassembled because a neighboring element is red
+ Green //!< entity does not need to be reassembled
+ };
+
BoxAssembler()
{
problemPtr_ = 0;
@@ -127,12 +141,30 @@ public:
ghostIndices_.push_back(vIdx);
}
};
+
+ int numVerts = gridView_().size(dim);
+ int numElems = gridView_().size(0);
+ residual_.resize(numVerts);
+
+ // initialize data needed for partial reassembly
+ if (enablePartialReassemble) {
+ vertexColor_.resize(numVerts);
+ elementColor_.resize(numElems);
+ }
+ std::fill(vertexColor_.begin(),
+ vertexColor_.end(),
+ Red);
+ std::fill(elementColor_.begin(),
+ elementColor_.end(),
+ Red);
}
void assemble(SolutionVector &u)
{
+ // assemble the global jacobian matrix
if (!reuseMatrix_) {
- *matrix_ = 0;
+ // we actually need to reassemle!
+ resetMatrix_();
gridOperatorSpace_->jacobian(u, *matrix_);
}
reuseMatrix_ = false;
@@ -148,67 +180,141 @@ public:
(*matrix_)[globI][globI] = Id;
}
- residual_.resize(u.size());
+ // calculate the global residual
residual_ = 0;
gridOperatorSpace_->residual(u, residual_);
+ }
-#if 0
- // rescale jacobian and right hand side to the largest
- // entry on the main diagonal block matrix
- typedef typename Matrix::RowIterator RowIterator;
- typedef typename Matrix::ColIterator ColIterator;
- typedef typename Matrix::block_type BlockType;
- const typename Matrix::block_type::size_type rowsInBlock = Matrix::block_type::rows;
- const typename Matrix::block_type::size_type colsInBlock = Matrix::block_type::cols;
- Scalar diagonalEntry[rowsInBlock];
- Vector diagonalEntries(*f);
- RowIterator endIBlock = matrix_->end();
- for (RowIterator iBlock = matrix_->begin(); iBlock != endIBlock; ++iBlock) {
- BlockType &diagBlock = (*iBlock)[iBlock.index()];
-
- for (int i = 0; i < rowsInBlock; ++i) {
- diagonalEntry[i] = 0;
- for (int j = 0; j < colsInBlock; ++j) {
- diagonalEntry[i] = std::max(diagonalEntry[i],
- std::abs(diagBlock[i][j]));
- }
-
- if (diagonalEntry[i] < 1e-14)
- diagonalEntry[i] = 1.0;
+ void setMatrixReuseable(bool yesno = true)
+ {
+ if (enableJacobianRecycling)
+ reuseMatrix_ = yesno;
+ }
- diagonalEntries[iBlock.index()][i] = diagonalEntry[i];
- }
+ void reassembleAll()
+ {
+ std::fill(vertexColor_.begin(),
+ vertexColor_.end(),
+ Red);
+ std::fill(elementColor_.begin(),
+ elementColor_.end(),
+ Red);
+ }
+
+ void markVertexRed(int globalVertIdx, bool yesno)
+ {
+ if (enablePartialReassemble)
+ vertexColor_[globalVertIdx] = yesno?Red:Green;
+ }
+
+ void computeColors()
+ {
+
+ if (!enablePartialReassemble)
+ return;
+
+ ElementIterator elemIt = gridView_().template begin<0>();
+ ElementIterator elemEndIt = gridView_().template end<0>();
+
+ // Mark all red elements
+ for (; elemIt != elemEndIt; ++elemIt) {
+ bool needReassemble = false;
+ int numVerts = elemIt->template count<dim>();
+ for (int i=0; i < numVerts; ++i) {
+ int globalI = vertexMapper_().map(*elemIt, i, dim);
+ if (vertexColor_[globalI] == Red) {
+ needReassemble = true;
+ break;
+ }
+ };
+
+ int globalElemIdx = elementMapper_().map(*elemIt);
+ elementColor_[globalElemIdx] = needReassemble?Red:Green;
}
-
- Dune::PDELab::AddDataHandle<GridFunctionSpace,Vector> adddh(*gridFunctionSpace_, diagonalEntries);
- if (gridFunctionSpace_->gridview().comm().size()>1)
- gridFunctionSpace_->gridview().communicate(adddh,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
- for (RowIterator iBlock = matrix_->begin(); iBlock != endIBlock; ++iBlock) {
- // divide right-hand side
- for (int i = 0; i < rowsInBlock; i++) {
- (*f)[iBlock.index()][i] /= diagonalEntries[iBlock.index()][i];
+ // Mark all yellow vertices
+ elemIt = gridView_().template begin<0>();
+ for (; elemIt != elemEndIt; ++elemIt) {
+ int elemIdx = this->elementMapper_().map(*elemIt);
+ if (elementColor_[elemIdx] == Green)
+ continue; // green elements do not tint vertices
+ // 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] = Yellow;
+ };
+ }
+
+ // Mark all yellow elements
+ int numReassemble = 0;
+ int numGreen = 0;
+ elemIt = gridView_().template begin<0>();
+ for (; elemIt != elemEndIt; ++elemIt) {
+ int elemIdx = this->elementMapper_().map(*elemIt);
+ if (elementColor_[elemIdx] == Red) {
+ ++ numReassemble;
+ continue; // element is red already!
}
- // divide row of the jacobian
- ColIterator endJBlock = iBlock->end();
- for (ColIterator jBlock = iBlock->begin(); jBlock != endJBlock; ++jBlock) {
- for (int i = 0; i < rowsInBlock; i++) {
- for (int j = 0; j < colsInBlock; j++) {
- (*jBlock)[i][j] /= diagonalEntries[iBlock.index()][i];
- }
+ 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] == Yellow) {
+ ++ numReassemble;
+ isYellow = true;
+ break;
}
+ };
+
+ if (isYellow) {
+ ++ isYellow;
+ elementColor_[elemIdx] = Yellow;
}
+
+ if (elementColor_[elemIdx] == Green)
+ ++ numGreen;
}
-#endif
+
}
-
- void setMatrixReuseable(bool yesno = true)
- { reuseMatrix_ = yesno; }
+ EntityColor vertexColor(const Element &element, int vertIdx) const
+ {
+ if (!enablePartialReassemble)
+ return Red; // reassemble unconditionally!
+
+ int globalIdx = vertexMapper_().map(element, vertIdx, dim);
+ return vertexColor_[globalIdx];
+ }
+
+ EntityColor vertexColor(int globalVertIdx) const
+ {
+ if (!enablePartialReassemble)
+ return Red; // reassemble unconditionally!
+ return vertexColor_[globalVertIdx];
+ }
+
+ EntityColor elementColor(const Element &element) const
+ {
+ if (!enablePartialReassemble)
+ return Red; // reassemble unconditionally!
+
+ int globalIdx = elementMapper_().map(element);
+ return elementColor_[globalIdx];
+ }
+ EntityColor elementColor(int globalElementIdx) const
+ {
+ if (!enablePartialReassemble)
+ return Red; // reassemble unconditionally!
+ return elementColor_[globalElementIdx];
+ }
+
const GridFunctionSpace& gridFunctionSpace() const
{
return *gridFunctionSpace_;
@@ -228,8 +334,38 @@ public:
{ return residual_; }
private:
+ void resetMatrix_()
+ {
+ if (!enablePartialReassemble) {
+ (*matrix_) = 0;
+ return;
+ }
+
+ // reset all entries corrosponding to a red vertex
+ for (int rowIdx = 0; rowIdx < matrix_->N(); ++rowIdx) {
+ if (vertexColor_[rowIdx] == Green)
+ continue; // the equations for this control volume are
+ // already below the treshold
+ // reset 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;
+ }
+ };
+
+ //printSparseMatrix(std::cout, *matrix_, "J", "row");
+ }
+
+ 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
@@ -248,6 +384,8 @@ private:
Matrix *matrix_;
bool reuseMatrix_;
+ std::vector<EntityColor> vertexColor_;
+ std::vector<EntityColor> elementColor_;
std::vector<int> ghostIndices_;
SolutionVector residual_;
diff --git a/dumux/boxmodels/common/pdelabboxlocaloperator.hh b/dumux/boxmodels/common/pdelabboxlocaloperator.hh
index 410d606c4d..e60f24c99f 100644
--- a/dumux/boxmodels/common/pdelabboxlocaloperator.hh
+++ b/dumux/boxmodels/common/pdelabboxlocaloperator.hh
@@ -37,6 +37,15 @@ class BoxLocalOperator
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags
{
+ // copying the local operator for PDELab is not a good idea
+ BoxLocalOperator(const BoxLocalOperator &);
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
+ enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
+
public:
// pattern assembly flags
enum { doPatternVolume = true };
@@ -44,10 +53,6 @@ public:
// residual assembly flags
enum { doAlphaVolume = true };
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
BoxLocalOperator(Model &model)
: model_(model)
@@ -60,6 +65,11 @@ public:
{
typedef typename LFSU::Traits::SizeType size_type;
+ //enum { Green = JacobianAssembler::Green };
+ //if (model_.jacobianAssembler().elementColor(eg.entity()) == Green)
+ // Green elements don't need to be reassembled
+ // return;
+
model_.localResidual().eval(eg.entity());
int numVertices = x.size()/numEq;
@@ -76,13 +86,15 @@ public:
Dune::PDELab::LocalMatrix<R>& mat) const
{
typedef typename LFSU::Traits::SizeType size_type;
-
+
model_.localJacobian().assemble(eg.entity());
int numVertices = x.size()/numEq;
- for (size_type j=0; j<lfsu.size(); j++)
- for (size_type i=0; i<lfsu.size(); i++)
- mat(i,j) = (model_.localJacobian().mat(i%numVertices,j%numVertices))[i/numVertices][j/numVertices];
+ for (size_type j=0; j<lfsu.size(); j++) {
+ for (size_type i=0; i<lfsu.size(); i++) {
+ mat(i,j) = (model_.localJacobian().mat(i%numVertices,j%numVertices))[i/numVertices][j/numVertices];
+ }
+ }
}
private:
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index 7fb4af4559..ff5ab54139 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -278,7 +278,6 @@ public:
int newtonNumSteps()
{ return numSteps_; }
-
/*!
* \brief Update the error of the solution compared to the
* previous iteration.
@@ -293,6 +292,7 @@ public:
int idxI = -1;
int idxJ = -1;
+ int numReassemble = 0;
for (int i = 0; i < int(uOld.size()); ++i) {
bool needReassemble = false;
for (int j = 0; j < FV::size; ++j) {
@@ -300,8 +300,9 @@ public:
=
std::abs(deltaU[i][j])
/ std::max(std::abs(uOld[i][j]), Scalar(1e-4));
- if (tmp > tolerance_/2)
+ if (tmp > tolerance_ / 10) {
needReassemble = true;
+ }
if (tmp > error_)
{
idxI = i;
@@ -310,11 +311,17 @@ public:
}
}
- //model_().jacobianAssembler().setReassembleVertex(i,
- // needReassemble);
+ if (needReassemble) {
+ ++ numReassemble;
+ }
+
+ model_().jacobianAssembler().markVertexRed(i,
+ needReassemble);
}
- //model_().jacobianAssembler().calcElementsToReassemble();
+ model_().jacobianAssembler().computeColors();
+
+ endIterMsg() << ", red vertices: " << numReassemble << "/" << uOld.size();
error_ = gridView_().comm().max(error_);
}
@@ -384,7 +391,7 @@ public:
writeConvergence_(uOld, deltaU);
newtonUpdateRelError(uOld, deltaU);
-
+
deltaU *= -1;
deltaU += uOld;
}
diff --git a/dumux/nonlinear/newtonmethod.hh b/dumux/nonlinear/newtonmethod.hh
index d37acd263c..30e733d2f3 100644
--- a/dumux/nonlinear/newtonmethod.hh
+++ b/dumux/nonlinear/newtonmethod.hh
@@ -28,6 +28,7 @@
#include <dumux/io/vtkmultiwriter.hh>
+#include <dune/istl/io.hh>
namespace Dumux
{
diff --git a/test/boxmodels/2p/test_2p.cc b/test/boxmodels/2p/test_2p.cc
index 50511015b9..39406638ca 100644
--- a/test/boxmodels/2p/test_2p.cc
+++ b/test/boxmodels/2p/test_2p.cc
@@ -186,8 +186,8 @@ int main(int argc, char** argv)
Dune::FieldVector<int,dim> res; // cell resolution
upperRight[0] = 6.0;
upperRight[1] = 4.0;
- res[0] = 48;
- res[1] = 32;
+ res[0] = 96;
+ res[1] = 64;
std::auto_ptr<Grid> grid(CreateGrid<Grid, Scalar>::create(upperRight, res));
diff --git a/test/boxmodels/2p2c/injectionproblem.hh b/test/boxmodels/2p2c/injectionproblem.hh
index 43aa9b6e36..ade58ca9d2 100644
--- a/test/boxmodels/2p2c/injectionproblem.hh
+++ b/test/boxmodels/2p2c/injectionproblem.hh
@@ -89,6 +89,9 @@ SET_BOOL_PROP(InjectionProblem, EnableGravity, true);
// Enable gravity
SET_INT_PROP(InjectionProblem, NewtonLinearSolverVerbosity, 0);
+
+SET_BOOL_PROP(InjectionProblem, EnableJacobianRecycling, true);
+SET_BOOL_PROP(InjectionProblem, EnablePartialReassemble, true);
}
--
GitLab