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

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
parent f74a044d
......@@ -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]
*/
......
......@@ -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.
*/
......
......@@ -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
......@@ -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()
{ }
......
......@@ -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();
};
/*!
......
......@@ -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),
......
......@@ -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);
// \}
}
......
......@@ -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_;
......
......@@ -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:
......
......@@ -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;