diff --git a/dune/faultnetworks/dgmgtransfer.hh b/dune/faultnetworks/dgmgtransfer.hh
index 641d3f8cdba07d63535d73a50687811b09603737..7201df5e094cd43045cebb355064cc7b23af1fc7 100644
--- a/dune/faultnetworks/dgmgtransfer.hh
+++ b/dune/faultnetworks/dgmgtransfer.hh
@@ -218,6 +218,24 @@ public:
{
GenericMultigridTransfer::restrict(prolongationMatrix_, f, t, -1);
}
+
+ template< class OperatorType>
+ void galerkinRestrictSetOccupation(const OperatorType& fineMat, OperatorType& coarseMat)
+ {
+ GenericMultigridTransfer::galerkinRestrictSetOccupation(prolongationMatrix_, fineMat, coarseMat);
+ }
+
+ template<class FineMatrixType, class CoarseMatrixType>
+ void galerkinRestrict(const FineMatrixType& fineMat, CoarseMatrixType& coarseMat)
+ {
+ GenericMultigridTransfer::galerkinRestrict(prolongationMatrix_, fineMat, coarseMat);
+ }
+
+ template <class BitVectorType>
+ void restrictToFathers(const BitVectorType& f, BitVectorType& t) const
+ {
+ GenericMultigridTransfer::restrictBitFieldToFathers(prolongationMatrix_, f, t, -1);
+ }
};
#endif
diff --git a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
index c4c2b04fa72c2ee04a0c0e48140516f913575d7c..a14f005352084a4c1dd86aa6a47c6fa60e08b2c0 100644
--- a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
@@ -21,14 +21,14 @@ public:
const Mode mode = Mode::additive) :
Base(levelInterfaceNetwork, patchLevelBasis, localAssembler, localInterfaceAssemblers, mode) {}
- void build() {
+ void preprocess() {
const int patchLevel = this->patchLevelBasis_.faultNetwork().level();
CellPatchFactory<BasisType> patchFactory(this->patchLevelBasis_, this->basis_, patchLevel, this->level_);
const auto& localToGlobal = patchFactory.getLocalToGlobal();
const auto& boundaryDofs = patchFactory.getBoundaryDofs();
VectorType rhs;
- rhs.resize(this->matrix_.N());
+ rhs.resize(this->mat_->N());
rhs = 0;
std::cout << "CellPatchPreconditioner::build() level: " << this->level_ << std::endl;
@@ -38,7 +38,7 @@ public:
for (size_t i=0; i<cellCount; i++) {
using LocalProblemType = typename std::remove_reference<decltype(*(this->localProblems_[i]))>::type;
- this->localProblems_[i] = std::make_shared<LocalProblemType>(this->matrix_, rhs, localToGlobal[i], boundaryDofs[i]);
+ this->localProblems_[i] = std::make_shared<LocalProblemType>(*this->mat_, rhs, localToGlobal[i], boundaryDofs[i]);
}
}
};
diff --git a/dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh b/dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh
index 2b7b4e2bb80f6ad5ec302097654dbb5bc320deb4..b831764cdbadc6bd7d6c55ea7a8c40e949ca8e03 100644
--- a/dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh
@@ -57,7 +57,7 @@ public:
assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
}
- void build() {
+ void preprocess() {
//printBasisDofLocation(basis_);
GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index c9fb33fdb8fce847ece8f595310aa741c60733e6..3d9a9526ad7f09db9d610e7a35cf8ad96625679a 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -44,7 +44,6 @@ protected:
size_t patchDepth_;
BoundaryMode boundaryMode_;
- MatrixType matrix_;
std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
public:
@@ -66,11 +65,10 @@ public:
setPatchDepth();
setBoundaryMode();
setDirection();
- setup();
}
// has to setup localProblems_
- virtual void build() = 0;
+ virtual void preprocess() = 0;
void setPatchDepth(const size_t patchDepth = 0) {
patchDepth_ = patchDepth;
@@ -98,6 +96,8 @@ public:
}
virtual void iterate() {
+ *(this->x_) = 0;
+
if (mode_ == ADDITIVE)
iterateAdd();
else
@@ -121,47 +121,12 @@ public:
}
private:
- void setup() {
- // assemble stiffness matrix for entire level including all faults
- GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
- globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
-
- // set boundary conditions
- Dune::BitSetVector<1> globalBoundaryDofs;
- BoundaryPatch<GridView> boundaryPatch(levelInterfaceNetwork_.levelGridView(), true);
- constructBoundaryDofs(boundaryPatch, basis_, globalBoundaryDofs);
-
- typedef typename MatrixType::row_type RowType;
- typedef typename RowType::ConstIterator ColumnIterator;
-
- for(size_t i=0; i<globalBoundaryDofs.size(); i++) {
- if(!globalBoundaryDofs[i][0])
- continue;
-
- RowType& row = matrix_[i];
-
- ColumnIterator cIt = row.begin();
- ColumnIterator cEndIt = row.end();
-
- for(; cIt!=cEndIt; ++cIt) {
- row[cIt.index()] = 0;
- }
- row[i] = 1;
- }
- }
-
void iterateAdd() {
- //*(this->x_) = 0;
-
VectorType it, x;
for (size_t i=0; i<localProblems_.size(); i++) {
localProblems_[i]->solve(it);
localProblems_[i]->prolong(it, x);
- /*if (i==5) {
- writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
- }*/
-
*(this->x_) += x;
}
}
@@ -178,48 +143,22 @@ private:
}
void iterateStep(size_t i) {
- //*(this->x_) = 0;
-
auto& localProblem = *localProblems_[i];
- VectorType v, v1, v2, v3;
- localProblem.restrict(*(this->x_), v);
-
- v1.resize(localProblem.getLocalMat().N());
- localProblem.getLocalMat().mv(v, v1);
- v3.resize(matrix_.N());
- matrix_.mv(*(this->x_), v3);
- localProblem.restrict(v3, v2);
+ // compute residual
+ VectorType residual = *this->rhs_;
+ this->mat_->mmv(*this->x_, residual);
- //print(v1, "v1");
- //print(v2, "v2");
+ // restrict residual
+ VectorType localResidual;
+ localProblem.restrict(residual, localResidual);
- VectorType r1;
- localProblem.restrict(*(this->rhs_), r1);
- Dune::MatrixVector::subtractProduct(r1, localProblem.getLocalMat(), v);
+ // solve local problem
+ localProblem.updateLocalRhs(localResidual);
- /*localProblem.updateLocalRhs(r);
-
- VectorType x;
+ // prolong
+ VectorType x, v;
localProblem.solve(x);
- v += x;
-
- localProblem.updateVector(v, *(this->x_));*/
-
- VectorType localR;
- VectorType r = *(this->rhs_);
- Dune::MatrixVector::subtractProduct(r, matrix_, *(this->x_));
- localProblem.restrict(r, localR);
-
- // print(r1, "r1");
- // print(localR, "localR");
-
- localProblem.updateLocalRhs(localR);
-
- VectorType x;
- localProblem.solve(x);
-
- //VectorType v;
localProblem.prolong(x, v);
*(this->x_) += v;
}
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner_working.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner_working.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c9fb33fdb8fce847ece8f595310aa741c60733e6
--- /dev/null
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner_working.hh
@@ -0,0 +1,229 @@
+#ifndef LEVEL_PATCH_PRECONDITIONER_HH
+#define LEVEL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+
+#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
+#include <dune/faultnetworks/localproblem.hh>
+#include <dune/faultnetworks/levelinterfacenetwork.hh>
+#include <dune/faultnetworks/utils/debugutils.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+
+template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+
+public:
+ enum Mode {ADDITIVE, MULTIPLICATIVE};
+ enum BoundaryMode {homogeneous, fromIterate};
+ enum Direction {FORWARD, BACKWARD, SYMMETRIC};
+
+protected:
+ typedef typename BasisType::GridView GridView;
+ typedef typename GridView::Grid GridType;
+
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
+ const BasisType& patchLevelBasis_;
+ const LocalAssembler& localAssembler_;
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
+
+ const Mode mode_;
+ Direction multDirection_;
+
+ const GridType& grid_;
+ const int level_;
+ const BasisType basis_;
+
+ size_t patchDepth_;
+ BoundaryMode boundaryMode_;
+
+ MatrixType matrix_;
+ std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
+
+public:
+ LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+ const BasisType& patchLevelBasis,
+ const LocalAssembler& localAssembler,
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+ const Mode mode = LevelPatchPreconditioner::Mode::ADDITIVE) :
+ levelInterfaceNetwork_(levelInterfaceNetwork),
+ patchLevelBasis_(patchLevelBasis),
+ localAssembler_(localAssembler),
+ localInterfaceAssemblers_(localInterfaceAssemblers),
+ mode_(mode),
+ grid_(levelInterfaceNetwork_.grid()),
+ level_(levelInterfaceNetwork_.level()),
+ basis_(levelInterfaceNetwork_)
+ {
+ assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
+ setPatchDepth();
+ setBoundaryMode();
+ setDirection();
+ setup();
+ }
+
+ // has to setup localProblems_
+ virtual void build() = 0;
+
+ void setPatchDepth(const size_t patchDepth = 0) {
+ patchDepth_ = patchDepth;
+ }
+
+ void setBoundaryMode(const BoundaryMode boundaryMode = LevelPatchPreconditioner::BoundaryMode::homogeneous) {
+ boundaryMode_ = boundaryMode;
+ }
+
+ virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
+ this->x_ = &x;
+ this->rhs_ = &rhs;
+ this->mat_ = Dune::stackobject_to_shared_ptr(mat);
+
+ for (size_t i=0; i<localProblems_.size(); i++) {
+ if (boundaryMode_ == BoundaryMode::homogeneous)
+ localProblems_[i]->updateRhs(rhs);
+ else
+ localProblems_[i]->updateRhsAndBoundary(rhs, x);
+ }
+ }
+
+ void setDirection(Direction dir = SYMMETRIC) {
+ multDirection_ = dir;
+ }
+
+ virtual void iterate() {
+ if (mode_ == ADDITIVE)
+ iterateAdd();
+ else
+ iterateMult();
+ }
+
+ const BasisType& basis() const {
+ return basis_;
+ }
+
+ const GridView& gridView() const {
+ return basis_.getGridView();
+ }
+
+ const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork() const {
+ return levelInterfaceNetwork_;
+ }
+
+ size_t size() const {
+ return localProblems_.size();
+ }
+
+private:
+ void setup() {
+ // assemble stiffness matrix for entire level including all faults
+ GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
+ globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
+
+ // set boundary conditions
+ Dune::BitSetVector<1> globalBoundaryDofs;
+ BoundaryPatch<GridView> boundaryPatch(levelInterfaceNetwork_.levelGridView(), true);
+ constructBoundaryDofs(boundaryPatch, basis_, globalBoundaryDofs);
+
+ typedef typename MatrixType::row_type RowType;
+ typedef typename RowType::ConstIterator ColumnIterator;
+
+ for(size_t i=0; i<globalBoundaryDofs.size(); i++) {
+ if(!globalBoundaryDofs[i][0])
+ continue;
+
+ RowType& row = matrix_[i];
+
+ ColumnIterator cIt = row.begin();
+ ColumnIterator cEndIt = row.end();
+
+ for(; cIt!=cEndIt; ++cIt) {
+ row[cIt.index()] = 0;
+ }
+ row[i] = 1;
+ }
+ }
+
+ void iterateAdd() {
+ //*(this->x_) = 0;
+
+ VectorType it, x;
+ for (size_t i=0; i<localProblems_.size(); i++) {
+ localProblems_[i]->solve(it);
+ localProblems_[i]->prolong(it, x);
+
+ /*if (i==5) {
+ writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
+ }*/
+
+ *(this->x_) += x;
+ }
+ }
+
+ void iterateMult() {
+ if (multDirection_ != BACKWARD) {
+ for (size_t i=0; i<localProblems_.size(); i++)
+ iterateStep(i);
+ }
+
+ if (multDirection_ != Direction::FORWARD)
+ for (size_t i=localProblems_.size()-1; i>=0 && i<localProblems_.size(); i--)
+ iterateStep(i);
+ }
+
+ void iterateStep(size_t i) {
+ //*(this->x_) = 0;
+
+ auto& localProblem = *localProblems_[i];
+
+ VectorType v, v1, v2, v3;
+ localProblem.restrict(*(this->x_), v);
+
+ v1.resize(localProblem.getLocalMat().N());
+ localProblem.getLocalMat().mv(v, v1);
+ v3.resize(matrix_.N());
+ matrix_.mv(*(this->x_), v3);
+ localProblem.restrict(v3, v2);
+
+ //print(v1, "v1");
+ //print(v2, "v2");
+
+ VectorType r1;
+ localProblem.restrict(*(this->rhs_), r1);
+ Dune::MatrixVector::subtractProduct(r1, localProblem.getLocalMat(), v);
+
+ /*localProblem.updateLocalRhs(r);
+
+ VectorType x;
+ localProblem.solve(x);
+ v += x;
+
+ localProblem.updateVector(v, *(this->x_));*/
+
+ VectorType localR;
+ VectorType r = *(this->rhs_);
+ Dune::MatrixVector::subtractProduct(r, matrix_, *(this->x_));
+ localProblem.restrict(r, localR);
+
+ // print(r1, "r1");
+ // print(localR, "localR");
+
+ localProblem.updateLocalRhs(localR);
+
+ VectorType x;
+ localProblem.solve(x);
+
+ //VectorType v;
+ localProblem.prolong(x, v);
+ *(this->x_) += v;
+ }
+};
+
+#endif
+
diff --git a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
index e4bf038e631ba11dc8f103e985b782ead0f786b6..2f12d043cca61bef96032619c7d7be320341b462 100644
--- a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
@@ -10,6 +10,9 @@
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+#include <dune/matrix-vector/genericvectortools.hh>
+#include <dune/matrix-vector/resize.hh>
+
#include <dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh>
#include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
#include <dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh>
@@ -45,19 +48,24 @@ private:
using MultDirection = typename LevelPatchPreconditionerType::Direction;
MultDirection multDirection_;
-
int itCounter_;
std::shared_ptr<LevelGlobalPreconditionerType> coarseGlobalPreconditioner_;
std::vector<std::shared_ptr<LevelPatchPreconditionerType>> levelPatchPreconditioners_;
- std::vector<VectorType> levelX_;
- std::vector<VectorType> levelRhs_;
std::shared_ptr<LevelInterfaceNetwork<GridView>> allFaultLevelInterfaceNetwork_;
std::shared_ptr<BasisType> allFaultBasis_;
+ // data hierarchy
+ std::vector<std::shared_ptr<const MatrixType> > matrixHierarchy_;
+ std::vector<BitVector*> ignoreNodesHierarchy_;
+ std::vector<std::shared_ptr<VectorType> > xHierarchy_;
+ std::vector<VectorType> rhsHierarchy_;
+
std::vector<std::shared_ptr<DGMGTransfer<BasisType>>> mgTransfer_;
+ bool preprocessCalled_;
+
public:
MultilevelPatchPreconditioner(const InterfaceNetwork<GridType>& interfaceNetwork,
const BitVector& activeLevels,
@@ -68,7 +76,8 @@ public:
activeLevels_(activeLevels),
localAssemblers_(localAssemblers),
localInterfaceAssemblers_(localInterfaceAssemblers),
- parset_(parset)
+ parset_(parset),
+ preprocessCalled_(false)
{
parseSettings();
@@ -77,14 +86,13 @@ public:
assert(activeLevels.size() == localAssemblers_.size() && activeLevels.size() == localInterfaceAssemblers_.size());
- // init level fault preconditioners and multigrid transfer
+ // init level patch preconditioners
levelPatchPreconditioners_.resize(0);
- mgTransfer_.resize(0);
for (size_t i=0; i<activeLevels_.size(); i++) {
if (activeLevels_[i][0]) {
// init global problem on coarsest level
- const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork.levelInterfaceNetwork(i);
+ const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork_.levelInterfaceNetwork(i);
coarseGlobalPreconditioner_ = std::make_shared<LevelGlobalPreconditionerType>(levelNetwork, *localAssemblers_[i], localInterfaceAssemblers_[i]);
break;
@@ -92,41 +100,93 @@ public:
}
auto maxLevel = setupLevelPreconditioners();
-
- levelX_.resize(levelPatchPreconditioners_.size()+1);
- levelRhs_.resize(levelPatchPreconditioners_.size()+1);
-
allFaultBasis_ = std::make_shared<BasisType>(interfaceNetwork_.levelInterfaceNetwork(maxLevel));
- // init multigrid transfer
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- mgTransfer_.push_back(std::make_shared<DGMGTransfer<BasisType>>(levelPatchPreconditioners_[i]->basis(), *allFaultBasis_));
- }
- mgTransfer_.push_back(std::make_shared<DGMGTransfer<BasisType>>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
-
itCounter_ = 0;
}
+ virtual ~MultilevelPatchPreconditioner() {
+ for (int i=0; i<int(ignoreNodesHierarchy_.size()-1); i++) {
+ if (ignoreNodesHierarchy_[i])
+ delete(ignoreNodesHierarchy_[i]);
+ }
+ }
+
+ virtual void setProblem(const MatrixType& mat,VectorType& x, const VectorType& rhs) {
+ LinearIterationStep<MatrixType, VectorType>::setProblem(mat, x, rhs);
+
+ preprocessCalled_ = false;
+ }
- virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
- this->x_ = &x;
+ void setRhs(const VectorType& rhs) {
this->rhs_ = &rhs;
- this->mat_ = Dune::stackobject_to_shared_ptr(mat);
+ rhsHierarchy_.back() = rhs;
+ }
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- mgTransfer_[i]->restrict(x, levelX_[i]);
- mgTransfer_[i]->restrict(rhs, levelRhs_[i]);
+ virtual void preprocess() {
+ if (preprocessCalled_)
+ std::cout << "Warning: Preprocess has already been called without calling iterate afterwards!" << std::endl;
+
+ preprocessCalled_ = true;
+
+ const int numLevels = levelPatchPreconditioners_.size()+1;
- levelPatchPreconditioners_[i]->setProblem(mat, levelX_[i], levelRhs_[i]);
+ // init transfer operators
+ mgTransfer_.resize(numLevels-1);
+ mgTransfer_[0] = std::make_shared<DGMGTransfer<BasisType>>(coarseGlobalPreconditioner_->basis(), levelPatchPreconditioners_[0]->basis());
+ for (size_t i=1; i<mgTransfer_.size(); i++) {
+ mgTransfer_[i] = std::make_shared<DGMGTransfer<BasisType>>(levelPatchPreconditioners_[i-1]->basis(), levelPatchPreconditioners_[i]->basis());
}
- size_t j = levelPatchPreconditioners_.size();
- mgTransfer_[j]->restrict(x, levelX_[j]);
- mgTransfer_[j]->restrict(rhs, levelRhs_[j]);
+ // resize hierarchy containers to number of levels
+ xHierarchy_.resize(numLevels);
+ rhsHierarchy_.resize(numLevels);
+ matrixHierarchy_.resize(numLevels);
+ ignoreNodesHierarchy_.resize(numLevels);
+
+ for (int i=0; i<numLevels-1; i++) {
+ matrixHierarchy_[i].reset();
+ xHierarchy_[i].reset();
+ ignoreNodesHierarchy_[i] = NULL;
+ }
+
+ matrixHierarchy_.back() = this->mat_;
+ xHierarchy_.back() = Dune::stackobject_to_shared_ptr(*(this->x_));
+ rhsHierarchy_.back() = *(this->rhs_);
+
+ // assemble hierarchy of matrices
+ for (int i=numLevels-2; i>=0; i--) {
+ matrixHierarchy_[i] = std::shared_ptr<MatrixType>(new MatrixType);
+ xHierarchy_[i] = std::shared_ptr<VectorType>(new VectorType);
+
+ mgTransfer_[i]->galerkinRestrictSetOccupation(*matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(matrixHierarchy_[i]));
+ mgTransfer_[i]->galerkinRestrict(*matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(matrixHierarchy_[i]));
+
+ // Set solution vector sizes for the lower levels
+ Dune::MatrixVector::resize(*xHierarchy_[i], *matrixHierarchy_[i]);
+ }
+
+ // set dirichlet bitfields
+ if (this->ignoreNodes_!=0) {
+ ignoreNodesHierarchy_[numLevels-1] = const_cast<BitVector*>(this->ignoreNodes_);
- coarseGlobalPreconditioner_->setProblem(mat, levelX_[j], levelRhs_[j]);
+ for (int i=numLevels-2; i>=0; --i) {
+ ignoreNodesHierarchy_[i] = new BitVector();
+ mgTransfer_[i]->restrictToFathers(*ignoreNodesHierarchy_[i+1], *ignoreNodesHierarchy_[i]);
+ }
+ } else
+ DUNE_THROW(SolverError, "We need a set of nodes to ignore");
+
+ coarseGlobalPreconditioner_->setProblem(*matrixHierarchy_[0], *xHierarchy_[0], rhsHierarchy_[0]);
+ coarseGlobalPreconditioner_->preprocess();
+
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ levelPatchPreconditioners_[i]->setProblem(*matrixHierarchy_[i+1], *xHierarchy_[i+1], rhsHierarchy_[i+1]);
+ levelPatchPreconditioners_[i]->preprocess();
+ }
}
+
virtual void iterate() {
if (mode_ == Mode::ADDITIVE)
iterateAdd();
@@ -134,14 +194,6 @@ public:
iterateMult();
}
- void build() {
- coarseGlobalPreconditioner_->build();
-
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- levelPatchPreconditioners_[i]->build();
- }
- }
-
std::shared_ptr<LevelPatchPreconditionerType> levelPatchPreconditioner(const int level) {
return levelPatchPreconditioners_[level];
}
@@ -178,32 +230,31 @@ private:
}
void iterateAdd() {
- VectorType x;
+ *(this->x_) = 0;
+
+ for (int i=mgTransfer_.size()-1; i>=0; i--) {
+ *xHierarchy_[i] = 0;
+ mgTransfer_[i]->restrict(rhsHierarchy_[i+1], rhsHierarchy_[i]);
+ }
+
+ // coarseGlobalPreconditioner_->setProblem(*matrixHierarchy_[0], *xHierarchy_[0], rhsHierarchy_[0]);
+ coarseGlobalPreconditioner_->iterate();
for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ // levelPatchPreconditioners_[i]->setProblem(*matrixHierarchy_[i+1], *xHierarchy_[i+1], rhsHierarchy_[i+1]);
levelPatchPreconditioners_[i]->iterate();
- const VectorType& it = levelPatchPreconditioners_[i]->getSol();
-
- mgTransfer_[i]->prolong(it, x);
- //writeToVTK(*allFaultBasis_, x, "/home/mi/podlesny/data/faultnetworks/sol/", "preconditioner_step_"+std::to_string(i));
+ VectorType x(xHierarchy_[i+1]->size());
+ mgTransfer_[i]->prolong(*xHierarchy_[i], x);
- *(this->x_) += x;
+ *xHierarchy_[i+1] += x;
}
-
- coarseGlobalPreconditioner_->iterate();
- const VectorType& it = coarseGlobalPreconditioner_->getSol();
-
- mgTransfer_[levelPatchPreconditioners_.size()]->prolong(it, x);
- *(this->x_) += x;
}
void iterateMult() {
- if (multDirection_ != MultDirection::FORWARD)
- for (size_t i=levelPatchPreconditioners_.size()-1; i>=0 && i<levelPatchPreconditioners_.size(); i--)
- iterateStep(i, MultDirection::BACKWARD);
+/*
size_t j = levelPatchPreconditioners_.size();
mgTransfer_[j]->restrict(*(this->x_), levelX_[j]);
@@ -232,31 +283,82 @@ private:
}
itCounter_++;
+
+ preprocessCalled_ = false;
+
+ if (multDirection_ != MultDirection::FORWARD)
+ for (size_t i=levelPatchPreconditioners_.size()-1; i>=0 && i<levelPatchPreconditioners_.size(); i--)
+ iterateStep(i, MultDirection::BACKWARD);
+
+ //presmoother_[level]->setProblem(*(this->matrixHierarchy_[level]), *x[level], rhs[level]);
+ //presmoother_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level];
+
+ // restriction
+
+ // compute residual
+ /* // fineResidual = rhs[level] - mat[level] * x[level];
+ VectorType fineResidual = rhs[level];
+ mat[level]->mmv(*x[level], fineResidual);
+
+ // restrict residual
+ this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]);
+
+
+ // Set Dirichlet values.
+ MatrixVector::Generic::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
+
+ // Choose all zeros as the initial correction
+ *x[level-1] = 0;
+
+ // ///////////////////////////////////////
+ // Recursively solve the coarser system
+ level--;
+ for (int i=0; i<mu_; i++)
+ iterate();
+ level++;
+
+ // ////////////////////////////////////////
+ // Prolong
+
+ // add correction to the presmoothed solution
+ VectorType tmp;
+ this->mgTransfer_[level-1]->prolong(*x[level-1], tmp);
+ *x[level] += tmp;
+ }
+
+ // Postsmoothing
+ postsmoother_[level]->setProblem(*(mat[level]), *x[level], rhs[level]);
+ postsmoother_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level]; */
}
- void iterateStep(const size_t i, const MultDirection dir) {
- //mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
+ void iterateStep(const size_t level, const MultDirection dir) {
+ /* levelPatchPreconditioners_[level]->setProblem(*matrixHierarchy_[level], *xHierarchy_[level], rhsHierarchy_[level]);
+ levelPatchPreconditioners_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level];
+ levelPatchPreconditioners_[level]->setDirection(dir);
- VectorType residual = *(this->rhs_);
- Dune::MatrixVector::subtractProduct(residual, *(this->mat_), *(this->x_));
+ // compute residual
+ // fineResidual = rhs[level] - mat[level] * x[level];
+ VectorType fineResidual = rhsHierarchy_[level];
+ matrixHierarchy_[level]->mmv(*xHierarchy_[level], fineResidual);
- VectorType localResidual;
- mgTransfer_[i]->restrict(residual, localResidual);
- //print(levelRhs_[i], "levelLocalCoarseRhs: ");
- //writeToVTK(levelPatchPreconditioners_[i]->basis(), levelRhs_[i], "/storage/mi/podlesny/data/faultnetworks/rhs/fine", "exactvertexdata_step_"+std::to_string(itCounter_));
+ // restrict residual
+ mgTransfer_[level-1]->restrict(fineResidual, rhsHierarchy_[level-1]);
+
+ // set dirichlet values
+ Dune::MatrixVector::Generic::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
+ // initial correction
+ *xHierarchy_[level-1] = 0;
- levelPatchPreconditioners_[i]->setProblem(*(this->mat_), levelX_[i], localResidual);
- levelPatchPreconditioners_[i]->setDirection(dir);
- levelPatchPreconditioners_[i]->iterate();
- const VectorType& it = levelPatchPreconditioners_[i]->getSol();
+ levelPatchPreconditioners_[level]->iterate();
+ const VectorType& it = levelPatchPreconditioners_[level]->getSol();
VectorType x;
mgTransfer_[i]->prolong(it, x);
*(this->x_) += x;
- //writeToVTK(*allFaultBasis_, x, "/storage/mi/podlesny/data/faultnetworks/fineIterates/multilevel", "exactvertexdata_step_"+std::to_string(itCounter_));
+ //writeToVTK(*allFaultBasis_, x, "/storage/mi/podlesny/data/faultnetworks/fineIterates/multilevel", "exactvertexdata_step_"+std::to_string(itCounter_)); */
}
auto setupLevelPreconditioners() {
diff --git a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner_working.hh b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner_working.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e4bf038e631ba11dc8f103e985b782ead0f786b6
--- /dev/null
+++ b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner_working.hh
@@ -0,0 +1,337 @@
+#ifndef MULTILEVEL_PATCH_PRECONDITIONER_HH
+#define MULTILEVEL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+
+#include <dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh>
+#include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
+#include <dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh>
+#include <dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh>
+#include <dune/faultnetworks/levelinterfacenetwork.hh>
+#include <dune/faultnetworks/interfacenetwork.hh>
+#include <dune/faultnetworks/dgmgtransfer.hh>
+
+#include <dune/faultnetworks/utils/debugutils.hh>
+
+template <class BasisType, class LocalOperatorAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+
+public:
+ using LevelPatchPreconditionerType = LevelPatchPreconditioner<BasisType, LocalOperatorAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+ using LevelGlobalPreconditionerType = LevelGlobalPreconditioner<BasisType, LocalOperatorAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+
+private:
+ using GridView = typename BasisType::GridView;
+ using GridType = typename GridView::Grid ;
+ using BitVector = typename Dune::BitSetVector<1>;
+
+ const InterfaceNetwork<GridType>& interfaceNetwork_;
+ const BitVector& activeLevels_;
+ const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers_;
+ const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers_;
+
+ // preconditioner parset settings
+ const Dune::ParameterTree& parset_;
+ using Mode = typename LevelPatchPreconditionerType::Mode;
+ Mode mode_;
+
+ using MultDirection = typename LevelPatchPreconditionerType::Direction;
+ MultDirection multDirection_;
+
+
+ int itCounter_;
+
+ std::shared_ptr<LevelGlobalPreconditionerType> coarseGlobalPreconditioner_;
+ std::vector<std::shared_ptr<LevelPatchPreconditionerType>> levelPatchPreconditioners_;
+ std::vector<VectorType> levelX_;
+ std::vector<VectorType> levelRhs_;
+
+ std::shared_ptr<LevelInterfaceNetwork<GridView>> allFaultLevelInterfaceNetwork_;
+ std::shared_ptr<BasisType> allFaultBasis_;
+
+ std::vector<std::shared_ptr<DGMGTransfer<BasisType>>> mgTransfer_;
+
+public:
+ MultilevelPatchPreconditioner(const InterfaceNetwork<GridType>& interfaceNetwork,
+ const BitVector& activeLevels,
+ const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers,
+ const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers,
+ const Dune::ParameterTree& parset) :
+ interfaceNetwork_(interfaceNetwork),
+ activeLevels_(activeLevels),
+ localAssemblers_(localAssemblers),
+ localInterfaceAssemblers_(localInterfaceAssemblers),
+ parset_(parset)
+ {
+ parseSettings();
+
+ if (activeLevels_.size() > (size_t) interfaceNetwork.maxLevel() +1)
+ DUNE_THROW(Dune::Exception, "MultilevelFaultPreconditioner: too many active levels; preconditioner supports at most (grid.maxLevel + 1) levels!");
+
+ assert(activeLevels.size() == localAssemblers_.size() && activeLevels.size() == localInterfaceAssemblers_.size());
+
+ // init level fault preconditioners and multigrid transfer
+ levelPatchPreconditioners_.resize(0);
+ mgTransfer_.resize(0);
+
+ for (size_t i=0; i<activeLevels_.size(); i++) {
+ if (activeLevels_[i][0]) {
+ // init global problem on coarsest level
+ const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork.levelInterfaceNetwork(i);
+ coarseGlobalPreconditioner_ = std::make_shared<LevelGlobalPreconditionerType>(levelNetwork, *localAssemblers_[i], localInterfaceAssemblers_[i]);
+
+ break;
+ }
+ }
+
+ auto maxLevel = setupLevelPreconditioners();
+
+ levelX_.resize(levelPatchPreconditioners_.size()+1);
+ levelRhs_.resize(levelPatchPreconditioners_.size()+1);
+
+ allFaultBasis_ = std::make_shared<BasisType>(interfaceNetwork_.levelInterfaceNetwork(maxLevel));
+
+ // init multigrid transfer
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ mgTransfer_.push_back(std::make_shared<DGMGTransfer<BasisType>>(levelPatchPreconditioners_[i]->basis(), *allFaultBasis_));
+ }
+ mgTransfer_.push_back(std::make_shared<DGMGTransfer<BasisType>>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
+
+ itCounter_ = 0;
+ }
+
+
+ virtual void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) {
+ this->x_ = &x;
+ this->rhs_ = &rhs;
+ this->mat_ = Dune::stackobject_to_shared_ptr(mat);
+
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ mgTransfer_[i]->restrict(x, levelX_[i]);
+ mgTransfer_[i]->restrict(rhs, levelRhs_[i]);
+
+ levelPatchPreconditioners_[i]->setProblem(mat, levelX_[i], levelRhs_[i]);
+ }
+
+ size_t j = levelPatchPreconditioners_.size();
+ mgTransfer_[j]->restrict(x, levelX_[j]);
+ mgTransfer_[j]->restrict(rhs, levelRhs_[j]);
+
+ coarseGlobalPreconditioner_->setProblem(mat, levelX_[j], levelRhs_[j]);
+ }
+
+ virtual void iterate() {
+ if (mode_ == Mode::ADDITIVE)
+ iterateAdd();
+ else
+ iterateMult();
+ }
+
+ void build() {
+ coarseGlobalPreconditioner_->build();
+
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ levelPatchPreconditioners_[i]->build();
+ }
+ }
+
+ std::shared_ptr<LevelPatchPreconditionerType> levelPatchPreconditioner(const int level) {
+ return levelPatchPreconditioners_[level];
+ }
+
+ const BasisType& basis() const {
+ return *allFaultBasis_;
+ }
+
+ size_t size() const {
+ return levelPatchPreconditioners_.size()+1;
+ }
+
+private:
+ void parseSettings() {
+ const auto mode = parset_.get<std::string>("mode");
+ if (mode == "ADDITIVE") {
+ mode_ = Mode::ADDITIVE;
+ } else if (mode == "MULTIPLICATIVE") {
+ mode_ = Mode::MULTIPLICATIVE;
+ } else {
+ DUNE_THROW(Dune::Exception, "MultilevelFaultPreconditioner::parseSettings unknown mode! Possible options: ADDITIVE , MULTIPLICATIVE");
+ }
+
+ auto multDirection = parset_.get<std::string>("multDirection");
+ if (multDirection == "FORWARD") {
+ multDirection_ = MultDirection::FORWARD;
+ } else if (multDirection == "BACKWARD") {
+ multDirection_ = MultDirection::BACKWARD;
+ } else if (multDirection == "SYMMETRIC") {
+ multDirection_ = MultDirection::SYMMETRIC;
+ } else {
+ DUNE_THROW(Dune::Exception, "MultilevelFaultPreconditioner::parseSettings unknown multDirection! Possible options: FORWARD , BACKWARD, SYMMETRIC");
+ }
+ }
+
+ void iterateAdd() {
+ VectorType x;
+
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ levelPatchPreconditioners_[i]->iterate();
+ const VectorType& it = levelPatchPreconditioners_[i]->getSol();
+
+ mgTransfer_[i]->prolong(it, x);
+
+ //writeToVTK(*allFaultBasis_, x, "/home/mi/podlesny/data/faultnetworks/sol/", "preconditioner_step_"+std::to_string(i));
+
+ *(this->x_) += x;
+ }
+
+ coarseGlobalPreconditioner_->iterate();
+ const VectorType& it = coarseGlobalPreconditioner_->getSol();
+
+ mgTransfer_[levelPatchPreconditioners_.size()]->prolong(it, x);
+ *(this->x_) += x;
+ }
+
+
+ void iterateMult() {
+
+ if (multDirection_ != MultDirection::FORWARD)
+ for (size_t i=levelPatchPreconditioners_.size()-1; i>=0 && i<levelPatchPreconditioners_.size(); i--)
+ iterateStep(i, MultDirection::BACKWARD);
+
+ size_t j = levelPatchPreconditioners_.size();
+ mgTransfer_[j]->restrict(*(this->x_), levelX_[j]);
+
+ //print(levelRhs_[j], "localCoarseRhs: ");
+ //writeToVTK(coarseGlobalPreconditioner_->basis(), levelRhs_[j], "/storage/mi/podlesny/data/faultnetworks/rhs/coarse", "exactvertexdata_step_"+std::to_string(itCounter_));
+
+ VectorType residual = *(this->rhs_);
+ Dune::MatrixVector::subtractProduct(residual, *(this->mat_), *(this->x_));
+
+ VectorType localResidual;
+ mgTransfer_[j]->restrict(residual, localResidual);
+
+ coarseGlobalPreconditioner_->setProblem(*(this->mat_), levelX_[j], localResidual);
+ coarseGlobalPreconditioner_->iterate();
+ const VectorType& it = coarseGlobalPreconditioner_->getSol();
+
+ VectorType x;
+ mgTransfer_[j]->prolong(it, x);
+ *(this->x_) += x;
+ //writeToVTK(*allFaultBasis_, x, "/storage/mi/podlesny/data/faultnetworks/coarseIterates/multilevel", "exactvertexdata_step_"+std::to_string(itCounter_));
+
+ if (multDirection_ != MultDirection::BACKWARD) {
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++)
+ iterateStep(i, MultDirection::FORWARD);
+ }
+
+ itCounter_++;
+ }
+
+ void iterateStep(const size_t i, const MultDirection dir) {
+ //mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
+
+ VectorType residual = *(this->rhs_);
+ Dune::MatrixVector::subtractProduct(residual, *(this->mat_), *(this->x_));
+
+ VectorType localResidual;
+ mgTransfer_[i]->restrict(residual, localResidual);
+ //print(levelRhs_[i], "levelLocalCoarseRhs: ");
+ //writeToVTK(levelPatchPreconditioners_[i]->basis(), levelRhs_[i], "/storage/mi/podlesny/data/faultnetworks/rhs/fine", "exactvertexdata_step_"+std::to_string(itCounter_));
+
+
+ levelPatchPreconditioners_[i]->setProblem(*(this->mat_), levelX_[i], localResidual);
+
+ levelPatchPreconditioners_[i]->setDirection(dir);
+ levelPatchPreconditioners_[i]->iterate();
+ const VectorType& it = levelPatchPreconditioners_[i]->getSol();
+
+ VectorType x;
+ mgTransfer_[i]->prolong(it, x);
+ *(this->x_) += x;
+
+ //writeToVTK(*allFaultBasis_, x, "/storage/mi/podlesny/data/faultnetworks/fineIterates/multilevel", "exactvertexdata_step_"+std::to_string(itCounter_));
+ }
+
+ auto setupLevelPreconditioners() {
+ const auto preconditionerType = parset_.get<std::string>("patch");
+ if (preconditionerType == "CELL") {
+ return setupCellPreconditioner();
+ } else if (preconditionerType == "SUPPORT") {
+ return setupSupportPreconditioner();
+ } else {
+ DUNE_THROW(Dune::Exception, "MultilevelFaultPreconditioner: unknown levelpatchpreconditioner type! Possible options: CELL , SUPPORT");
+ }
+ }
+
+ auto setupCellPreconditioner() {
+ using CellPreconditioner = CellPatchPreconditioner<BasisType, LocalOperatorAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+
+ int maxLevel = 0;
+
+ bool skip = true;
+ for (size_t i=0; i<activeLevels_.size(); i++) {
+ if (activeLevels_[i][0] && skip) {
+ skip = false;
+ continue;
+ }
+
+ if (activeLevels_[i][0]) {
+ const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork_.levelInterfaceNetwork(i);
+
+ if (levelPatchPreconditioners_.size() == 0) {
+ levelPatchPreconditioners_.push_back(std::make_shared<CellPreconditioner>(levelNetwork, coarseGlobalPreconditioner_->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], mode_));
+ } else {
+ levelPatchPreconditioners_.push_back(std::make_shared<CellPreconditioner>(levelNetwork, levelPatchPreconditioners_.back()->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], mode_));
+ levelPatchPreconditioners_.back()->setBoundaryMode(LevelPatchPreconditionerType::BoundaryMode::homogeneous);
+ }
+
+ maxLevel = i;
+ }
+ }
+
+ return maxLevel;
+ }
+
+ auto setupSupportPreconditioner() {
+ using SupportPreconditioner = SupportPatchPreconditioner<BasisType, LocalOperatorAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+
+ const auto patchDepth = parset_.get<size_t>("patchDepth");
+
+ int maxLevel = 0;
+
+ bool skip = true;
+ for (size_t i=0; i<activeLevels_.size(); i++) {
+ if (skip && activeLevels_[i][0]) {
+ skip = false;
+ continue;
+ }
+
+ if (activeLevels_[i][0]) {
+ // init local patch preconditioners on each level
+ const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork_.levelInterfaceNetwork(i);
+
+ if (levelPatchPreconditioners_.size() == 0) {
+ levelPatchPreconditioners_.push_back(std::make_shared<SupportPreconditioner>(levelNetwork, coarseGlobalPreconditioner_->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], mode_));
+ } else {
+ levelPatchPreconditioners_.push_back(std::make_shared<SupportPreconditioner>(levelNetwork, levelPatchPreconditioners_.back()->basis(), *localAssemblers_[i], localInterfaceAssemblers_[i], mode_));
+ levelPatchPreconditioners_.back()->setPatchDepth(patchDepth);
+ levelPatchPreconditioners_.back()->setBoundaryMode(LevelPatchPreconditionerType::BoundaryMode::homogeneous);
+ }
+
+ maxLevel = i;
+ }
+ }
+
+ return maxLevel;
+ }
+};
+
+#endif
+
diff --git a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
index 68f544cabfcd176529da1164c9b1749f4a74389e..dc4064e81981890575975847ab157993829b3c0f 100644
--- a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
@@ -26,7 +26,7 @@ public:
Base(levelInterfaceNetwork, patchLevelBasis, localAssembler, localInterfaceAssemblers, mode)
{}
- void build() {
+ void preprocess() {
// init vertexInElements
const int dim = GridType::dimension;
using Element = typename GridType::template Codim<0>::Entity;
@@ -60,11 +60,11 @@ public:
const auto& boundaryDofs = patchFactory.getBoundaryDofs();
VectorType rhs;
- rhs.resize(this->matrix_.N());
+ rhs.resize(this->mat_->N());
rhs = 0;
using LocalProblemType = typename std::remove_reference<decltype(*(this->localProblems_[i]))>::type;
- this->localProblems_[i] = std::make_shared<LocalProblemType>(this->matrix_, rhs, localToGlobal, boundaryDofs);
+ this->localProblems_[i] = std::make_shared<LocalProblemType>(*this->mat_, rhs, localToGlobal, boundaryDofs);
}
}
};
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.cc b/src/cantorfaultnetworks/cantorfaultnetwork.cc
index e2ed377d029c7229d6f8ff5125a8c3b088b26ecf..90a9b2b9a97e37b4e5e625b4c1311e15ab5fc85d 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.cc
@@ -317,7 +317,6 @@ int main(int argc, char** argv) { try
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
}*/
- preconditioner.build();
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.parset b/src/cantorfaultnetworks/cantorfaultnetwork.parset
index fe8182cf83ad0868dc2354c936de20b271a04bac..0af54c6260ea293ff2e6e68bb3eb56d8c2c31744 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.parset
@@ -3,19 +3,19 @@ resultPath = ../cantorfaultnetworks/results/
[preconditioner]
patch = SUPPORT # CELL , SUPPORT
-mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
+mode = ADDITIVE # ADDITIVE , MULTIPLICATIVE
multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1
###########################################
[problem0]
-oscDataFile = oscDataLaplace32.mat
+oscDataFile = oscDataLaplace16.mat
# level resolution in 2^(-...)
-coarseResolution = 0
-fineResolution = 5
-exactResolution = 8
+coarseResolution = 3
+fineResolution = 4
+exactResolution = 5
minCantorResolution = 0
penaltyFactor = 1
diff --git a/src/cantorfaultnetworks/results/sparse/compare.txt b/src/cantorfaultnetworks/results/sparse/compare.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7202b5a59ffa6ce289ec4968d8f27b4d49415495
--- /dev/null
+++ b/src/cantorfaultnetworks/results/sparse/compare.txt
@@ -0,0 +1,114 @@
+[problem0]
+
+local Cantor-type network
+
+coarseCantorLevel = 1
+fineCantorLevel = 2
+maxCantorLevel = 3
+
+grid resolutions: 4^{-k}, k = 1, 2, 3
+
+K = 2 (fineCantorLevel)
+iterate in step \nu: u_K^{(\nu)}
+exact solutions: u_{K+1}, u_{K}
+
+in iteration step \nu:
+----------------------
+discError: || u_{K+1} - u_{K}^{(\nu)} || / || u_{K+1} - u_{K}^{(*)} ||, where u_{K}^{(*)} is initial iterate given by u_{K-1}
+rate: \rho^{(\nu)} = || u_K - u_K^{(\nu)} || / || u_K - u_K^{(\nu-1)} ||
+criterion: || u_{K} - u_{K}^{(\nu)} || \leq || u_{K+1} - u_{K} ||
+
+totalConvRate in step \nu: \delta^{(\nu)} = (\Pi_{i=0}^{\nu} \rho^{(i)})^{1/ \nu}
+
+maxTotalConvRate: max{\delta^{(\nu)} : 0 \leq \nu \leq 7}
+totalConvRate: \delta^{(7)}
+
+
+2-level methods:
+
+
+--------------------------------------------------------------------------------------
+[preconditioner]
+patch = SUPPORT
+mode = ADDITIVE
+patchDepth = 1
+
+--- CGSolver ---
+ iter discError rate criterion
+----------------------------------------
+ 0 4.4748610e-01 0.38363 0
+ 1 2.6385601e-01 0.23138 1
+ 2 2.5082680e-01 0.30402 1
+ 3 2.4957306e-01 0.28518 1
+ 4 2.4947130e-01 0.29271 1
+ 5 2.4946257e-01 0.29074 1
+ 6 2.4946183e-01 0.29931 1
+ 7 2.4946177e-01 0.25816 1
+maxTotalConvRate: 0.3836291, totalConvRate: 0.2904196, 8 iterations performed
+--------------------------------------------------------------------------------------
+
+
+
+--------------------------------------------------------------------------------------
+[preconditioner]
+patch = SUPPORT
+mode = MULTIPLICATIVE
+multDirection = SYMMETRIC
+
+--- CGSolver ---
+ iter discError rate criterion
+----------------------------------------
+ 0 2.6563220e-01 0.09424 1
+ 1 2.4973859e-01 0.12881 1
+ 2 2.4947325e-01 0.20368 1
+ 3 2.4946253e-01 0.25855 1
+ 4 2.4946184e-01 0.32400 1
+ 5 2.4946177e-01 0.37774 1
+ 6 2.4946176e-01 0.42970 1
+ 7 2.4946176e-01 0.47934 1
+maxTotalConvRate: 0.2517121, totalConvRate: 0.2517121, 8 iterations performed
+--------------------------------------------------------------------------------------
+
+
+
+--------------------------------------------------------------------------------------
+[preconditioner]
+patch = CELL
+mode = ADDITIVE
+patchDepth = 1
+
+--- CGSolver ---
+ iter discError rate criterion
+----------------------------------------
+ 0 3.6789829e-01 0.27923 0
+ 1 2.9049278e-01 0.55046 1
+ 2 2.6631671e-01 0.62641 1
+ 3 2.5254058e-01 0.42165 1
+ 4 2.5012727e-01 0.46381 1
+ 5 2.4953864e-01 0.33968 1
+ 6 2.4948824e-01 0.58679 1
+ 7 2.4946635e-01 0.41623 1
+maxTotalConvRate: 0.4583344, totalConvRate: 0.4458779, 8 iterations performed
+--------------------------------------------------------------------------------------
+
+
+
+--------------------------------------------------------------------------------------
+[preconditioner]
+patch = CELL
+mode = MULTIPLICATIVE
+multDirection = SYMMETRIC
+
+--- CGSolver ---
+ iter discError rate criterion
+----------------------------------------
+ 0 2.8217734e-01 0.13619 1
+ 1 2.5034252e-01 0.15909 1
+ 2 2.4948254e-01 0.15347 1
+ 3 2.4946213e-01 0.13378 1
+ 4 2.4946178e-01 0.19682 1
+ 5 2.4946176e-01 0.10302 1
+ 6 2.4946176e-01 0.13350 1
+ 7 2.4946176e-01 0.15741 1
+maxTotalConvRate: 0.1543307, totalConvRate: 0.1444492, 8 iterations performed
+--------------------------------------------------------------------------------------
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
index 8c1da5ba1b467ba14f6d0460745484a40cfd8bf9..8e7ece9cb6a1510149ae0cbab36b8f078e6b83d9 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
@@ -314,7 +314,6 @@ int main(int argc, char** argv) { try
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
}*/
- preconditioner.build();
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
index 456c5bfc2935601ddc4e07e80fd4070fd2f9aff7..95857a7b236331928f4af765f5326e537a5e5d01 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
@@ -3,7 +3,7 @@ resultPath = ../cantorfaultnetworks/results/sparse/
[preconditioner]
patch = SUPPORT # CELL , SUPPORT
-mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
+mode = ADDITIVE # ADDITIVE , MULTIPLICATIVE
multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1
diff --git a/src/geofaultnetworks/geofaultnetwork.cc b/src/geofaultnetworks/geofaultnetwork.cc
index ddad719dd1e090bb824d67e40fcd2e11c6d71755..03f1b84e7debcfbbaf6007cb68450f3e9878d159 100644
--- a/src/geofaultnetworks/geofaultnetwork.cc
+++ b/src/geofaultnetworks/geofaultnetwork.cc
@@ -350,7 +350,6 @@ int main(int argc, char** argv) { try
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
}*/
- preconditioner.build();
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
diff --git a/src/geofaultnetworks/rockfaultnetwork.cc b/src/geofaultnetworks/rockfaultnetwork.cc
index cf0c2c3f91d020daeff3d5a0f66e42caf5d847f5..777b2e11cbfbb32089f955dfd0dee024c8fbd301 100644
--- a/src/geofaultnetworks/rockfaultnetwork.cc
+++ b/src/geofaultnetworks/rockfaultnetwork.cc
@@ -324,7 +324,6 @@ int main(int argc, char** argv) { try
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
}*/
- preconditioner.build();
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
std::cout << std::endl << std::endl;
diff --git a/src/geofaultnetworks/rockfaultnetwork.parset b/src/geofaultnetworks/rockfaultnetwork.parset
index a977724af639b1d266215f4e4b1775f592515c43..7c93947185ff266288b8af4fd7ef2ce274843c89 100644
--- a/src/geofaultnetworks/rockfaultnetwork.parset
+++ b/src/geofaultnetworks/rockfaultnetwork.parset
@@ -17,9 +17,9 @@ patchDepth = 1
oscDataFile = oscDataLaplace4.mat
# level resolution in 2^(-...)
-coarseResolution = 5
-fineResolution = 8
-exactResolution = 9
+coarseResolution = 3
+fineResolution = 4
+exactResolution = 5
penaltyFactor = 1