diff --git a/dune/faultnetworks/dgmgtransfer.hh b/dune/faultnetworks/dgmgtransfer.hh
index 7201df5e094cd43045cebb355064cc7b23af1fc7..641d3f8cdba07d63535d73a50687811b09603737 100644
--- a/dune/faultnetworks/dgmgtransfer.hh
+++ b/dune/faultnetworks/dgmgtransfer.hh
@@ -218,24 +218,6 @@ 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 a14f005352084a4c1dd86aa6a47c6fa60e08b2c0..c4c2b04fa72c2ee04a0c0e48140516f913575d7c 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 preprocess() {
+ void build() {
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->mat_->N());
+ rhs.resize(this->matrix_.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->mat_, rhs, localToGlobal[i], boundaryDofs[i]);
+ this->localProblems_[i] = std::make_shared<LocalProblemType>(this->matrix_, rhs, localToGlobal[i], boundaryDofs[i]);
}
}
};
diff --git a/dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh b/dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh
index b831764cdbadc6bd7d6c55ea7a8c40e949ca8e03..2b7b4e2bb80f6ad5ec302097654dbb5bc320deb4 100644
--- a/dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelglobalpreconditioner.hh
@@ -57,7 +57,7 @@ public:
assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
}
- void preprocess() {
+ void build() {
//printBasisDofLocation(basis_);
GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index 3d9a9526ad7f09db9d610e7a35cf8ad96625679a..c9fb33fdb8fce847ece8f595310aa741c60733e6 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -44,6 +44,7 @@ protected:
size_t patchDepth_;
BoundaryMode boundaryMode_;
+ MatrixType matrix_;
std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
public:
@@ -65,10 +66,11 @@ public:
setPatchDepth();
setBoundaryMode();
setDirection();
+ setup();
}
// has to setup localProblems_
- virtual void preprocess() = 0;
+ virtual void build() = 0;
void setPatchDepth(const size_t patchDepth = 0) {
patchDepth_ = patchDepth;
@@ -96,8 +98,6 @@ public:
}
virtual void iterate() {
- *(this->x_) = 0;
-
if (mode_ == ADDITIVE)
iterateAdd();
else
@@ -121,12 +121,47 @@ 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;
}
}
@@ -143,22 +178,48 @@ private:
}
void iterateStep(size_t i) {
+ //*(this->x_) = 0;
+
auto& localProblem = *localProblems_[i];
- // compute residual
- VectorType residual = *this->rhs_;
- this->mat_->mmv(*this->x_, residual);
+ 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);
- // restrict residual
- VectorType localResidual;
- localProblem.restrict(residual, localResidual);
+ //print(v1, "v1");
+ //print(v2, "v2");
- // solve local problem
- localProblem.updateLocalRhs(localResidual);
+ VectorType r1;
+ localProblem.restrict(*(this->rhs_), r1);
+ Dune::MatrixVector::subtractProduct(r1, localProblem.getLocalMat(), v);
- // prolong
- VectorType x, 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;
}
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner_working.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner_working.hh
deleted file mode 100644
index c9fb33fdb8fce847ece8f595310aa741c60733e6..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner_working.hh
+++ /dev/null
@@ -1,229 +0,0 @@
-#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 2f12d043cca61bef96032619c7d7be320341b462..e4bf038e631ba11dc8f103e985b782ead0f786b6 100644
--- a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
@@ -10,9 +10,6 @@
#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>
@@ -48,24 +45,19 @@ 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,
@@ -76,8 +68,7 @@ public:
activeLevels_(activeLevels),
localAssemblers_(localAssemblers),
localInterfaceAssemblers_(localInterfaceAssemblers),
- parset_(parset),
- preprocessCalled_(false)
+ parset_(parset)
{
parseSettings();
@@ -86,13 +77,14 @@ public:
assert(activeLevels.size() == localAssemblers_.size() && activeLevels.size() == localInterfaceAssemblers_.size());
- // init level patch preconditioners
+ // 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);
+ const LevelInterfaceNetwork<GridView>& levelNetwork = interfaceNetwork.levelInterfaceNetwork(i);
coarseGlobalPreconditioner_ = std::make_shared<LevelGlobalPreconditionerType>(levelNetwork, *localAssemblers_[i], localInterfaceAssemblers_[i]);
break;
@@ -100,93 +92,41 @@ public:
}
auto maxLevel = setupLevelPreconditioners();
- allFaultBasis_ = std::make_shared<BasisType>(interfaceNetwork_.levelInterfaceNetwork(maxLevel));
-
- 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;
- }
-
- void setRhs(const VectorType& rhs) {
- this->rhs_ = &rhs;
- rhsHierarchy_.back() = rhs;
- }
-
- virtual void preprocess() {
- if (preprocessCalled_)
- std::cout << "Warning: Preprocess has already been called without calling iterate afterwards!" << std::endl;
- preprocessCalled_ = true;
+ levelX_.resize(levelPatchPreconditioners_.size()+1);
+ levelRhs_.resize(levelPatchPreconditioners_.size()+1);
- const int numLevels = levelPatchPreconditioners_.size()+1;
+ allFaultBasis_ = std::make_shared<BasisType>(interfaceNetwork_.levelInterfaceNetwork(maxLevel));
- // 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());
+ // 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_));
- // 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;
- }
+ itCounter_ = 0;
+ }
- 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);
+ 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);
- 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]));
+ for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ mgTransfer_[i]->restrict(x, levelX_[i]);
+ mgTransfer_[i]->restrict(rhs, levelRhs_[i]);
- // Set solution vector sizes for the lower levels
- Dune::MatrixVector::resize(*xHierarchy_[i], *matrixHierarchy_[i]);
+ levelPatchPreconditioners_[i]->setProblem(mat, levelX_[i], levelRhs_[i]);
}
- // set dirichlet bitfields
- if (this->ignoreNodes_!=0) {
- ignoreNodesHierarchy_[numLevels-1] = const_cast<BitVector*>(this->ignoreNodes_);
-
- 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();
+ size_t j = levelPatchPreconditioners_.size();
+ mgTransfer_[j]->restrict(x, levelX_[j]);
+ mgTransfer_[j]->restrict(rhs, levelRhs_[j]);
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- levelPatchPreconditioners_[i]->setProblem(*matrixHierarchy_[i+1], *xHierarchy_[i+1], rhsHierarchy_[i+1]);
- levelPatchPreconditioners_[i]->preprocess();
- }
+ coarseGlobalPreconditioner_->setProblem(mat, levelX_[j], levelRhs_[j]);
}
-
virtual void iterate() {
if (mode_ == Mode::ADDITIVE)
iterateAdd();
@@ -194,6 +134,14 @@ 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];
}
@@ -230,31 +178,32 @@ private:
}
void iterateAdd() {
- *(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();
+ VectorType x;
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);
- VectorType x(xHierarchy_[i+1]->size());
- mgTransfer_[i]->prolong(*xHierarchy_[i], x);
+ //writeToVTK(*allFaultBasis_, x, "/home/mi/podlesny/data/faultnetworks/sol/", "preconditioner_step_"+std::to_string(i));
- *xHierarchy_[i+1] += x;
+ *(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]);
@@ -283,82 +232,31 @@ 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 level, const MultDirection dir) {
- /* levelPatchPreconditioners_[level]->setProblem(*matrixHierarchy_[level], *xHierarchy_[level], rhsHierarchy_[level]);
- levelPatchPreconditioners_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level];
- levelPatchPreconditioners_[level]->setDirection(dir);
+ void iterateStep(const size_t i, const MultDirection dir) {
+ //mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
- // compute residual
- // fineResidual = rhs[level] - mat[level] * x[level];
- VectorType fineResidual = rhsHierarchy_[level];
- matrixHierarchy_[level]->mmv(*xHierarchy_[level], fineResidual);
-
- // restrict residual
- mgTransfer_[level-1]->restrict(fineResidual, rhsHierarchy_[level-1]);
+ VectorType residual = *(this->rhs_);
+ Dune::MatrixVector::subtractProduct(residual, *(this->mat_), *(this->x_));
- // set dirichlet values
- Dune::MatrixVector::Generic::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
+ 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_));
- // initial correction
- *xHierarchy_[level-1] = 0;
+ levelPatchPreconditioners_[i]->setProblem(*(this->mat_), levelX_[i], localResidual);
- levelPatchPreconditioners_[level]->iterate();
- const VectorType& it = levelPatchPreconditioners_[level]->getSol();
+ 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_)); */
+ //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
deleted file mode 100644
index e4bf038e631ba11dc8f103e985b782ead0f786b6..0000000000000000000000000000000000000000
--- a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner_working.hh
+++ /dev/null
@@ -1,337 +0,0 @@
-#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 dc4064e81981890575975847ab157993829b3c0f..68f544cabfcd176529da1164c9b1749f4a74389e 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 preprocess() {
+ void build() {
// 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->mat_->N());
+ rhs.resize(this->matrix_.N());
rhs = 0;
using LocalProblemType = typename std::remove_reference<decltype(*(this->localProblems_[i]))>::type;
- this->localProblems_[i] = std::make_shared<LocalProblemType>(*this->mat_, rhs, localToGlobal, boundaryDofs);
+ this->localProblems_[i] = std::make_shared<LocalProblemType>(this->matrix_, rhs, localToGlobal, boundaryDofs);
}
}
};
diff --git a/src/cantorfaultnetworks/cantorfaultnetwork.cc b/src/cantorfaultnetworks/cantorfaultnetwork.cc
index 90a9b2b9a97e37b4e5e625b4c1311e15ab5fc85d..e2ed377d029c7229d6f8ff5125a8c3b088b26ecf 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.cc
@@ -317,6 +317,7 @@ 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 0af54c6260ea293ff2e6e68bb3eb56d8c2c31744..fe8182cf83ad0868dc2354c936de20b271a04bac 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.parset
@@ -3,19 +3,19 @@ resultPath = ../cantorfaultnetworks/results/
[preconditioner]
patch = SUPPORT # CELL , SUPPORT
-mode = ADDITIVE # ADDITIVE , MULTIPLICATIVE
+mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1
###########################################
[problem0]
-oscDataFile = oscDataLaplace16.mat
+oscDataFile = oscDataLaplace32.mat
# level resolution in 2^(-...)
-coarseResolution = 3
-fineResolution = 4
-exactResolution = 5
+coarseResolution = 0
+fineResolution = 5
+exactResolution = 8
minCantorResolution = 0
penaltyFactor = 1
diff --git a/src/cantorfaultnetworks/results/sparse/compare.txt b/src/cantorfaultnetworks/results/sparse/compare.txt
deleted file mode 100644
index 7202b5a59ffa6ce289ec4968d8f27b4d49415495..0000000000000000000000000000000000000000
--- a/src/cantorfaultnetworks/results/sparse/compare.txt
+++ /dev/null
@@ -1,114 +0,0 @@
-[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 8e7ece9cb6a1510149ae0cbab36b8f078e6b83d9..8c1da5ba1b467ba14f6d0460745484a40cfd8bf9 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
@@ -314,6 +314,7 @@ 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 95857a7b236331928f4af765f5326e537a5e5d01..456c5bfc2935601ddc4e07e80fd4070fd2f9aff7 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 = ADDITIVE # ADDITIVE , MULTIPLICATIVE
+mode = MULTIPLICATIVE # ADDITIVE , MULTIPLICATIVE
multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1
diff --git a/src/geofaultnetworks/geofaultnetwork.cc b/src/geofaultnetworks/geofaultnetwork.cc
index 03f1b84e7debcfbbaf6007cb68450f3e9878d159..ddad719dd1e090bb824d67e40fcd2e11c6d71755 100644
--- a/src/geofaultnetworks/geofaultnetwork.cc
+++ b/src/geofaultnetworks/geofaultnetwork.cc
@@ -350,6 +350,7 @@ 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 777b2e11cbfbb32089f955dfd0dee024c8fbd301..cf0c2c3f91d020daeff3d5a0f66e42caf5d847f5 100644
--- a/src/geofaultnetworks/rockfaultnetwork.cc
+++ b/src/geofaultnetworks/rockfaultnetwork.cc
@@ -324,6 +324,7 @@ 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 7c93947185ff266288b8af4fd7ef2ce274843c89..a977724af639b1d266215f4e4b1775f592515c43 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 = 3
-fineResolution = 4
-exactResolution = 5
+coarseResolution = 5
+fineResolution = 8
+exactResolution = 9
penaltyFactor = 1