From 4c67f92379ca1bec8baa43acddc19d53712af843 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Tue, 2 Jun 2020 13:07:45 +0200
Subject: [PATCH] .
---
.../cellpatchpreconditioner.hh | 4 +-
.../levelpatchpreconditioner.hh | 85 +++++++-------
.../multilevelpatchpreconditioner.hh | 111 +++++++++++-------
.../supportpatchpreconditioner.hh | 3 +-
src/cantorfaultnetworks/cantorfaultnetwork.cc | 13 +-
.../sparsecantorfaultnetwork.cc | 3 +-
.../sparsecantorfaultnetwork.parset | 4 +-
src/geofaultnetworks/geofaultnetwork.cc | 11 +-
8 files changed, 130 insertions(+), 104 deletions(-)
diff --git a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
index e938979..8fc86fc 100644
--- a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
@@ -33,8 +33,10 @@ public:
size_t cellCount = localToGlobal.size();
this->localProblems_.resize(cellCount);
+
for (size_t i=0; i<cellCount; i++) {
- this->localProblems_[i] = std::make_shared<decltype(*(this->localProblems_[i]))>(this->matrix_, rhs, localToGlobal[i], boundaryDofs[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]);
}
}
};
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index 331885d..d8c3090 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -25,14 +25,16 @@ public:
enum BoundaryMode {homogeneous, fromIterate};
enum Direction {FORWARD, BACKWARD, SYMMETRIC};
-private:
+protected:
typedef typename BasisType::GridView GridView;
typedef typename GridView::Grid GridType;
const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork_;
const LocalAssembler& localAssembler_;
const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers_;
+
const Mode mode_;
+ Direction multDirection_;
const GridType& grid_;
const int level_;
@@ -44,40 +46,6 @@ private:
MatrixType matrix_;
std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
-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;
- }
- }
-
-public:
- // has to setup localProblems_
- virtual void build();
-
public:
LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
const LocalAssembler& localAssembler,
@@ -94,9 +62,13 @@ public:
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;
}
@@ -118,11 +90,15 @@ public:
}
}
- virtual void iterate(Direction dir = SYMMETRIC) {
+ void setDirection(Direction dir = SYMMETRIC) {
+ multDirection_ = dir;
+ }
+
+ virtual void iterate() {
if (mode_ == ADDITIVE)
iterateAdd();
else
- iterateMult(dir);
+ iterateMult();
}
const BasisType& basis() const {
@@ -142,6 +118,35 @@ 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;
@@ -158,13 +163,13 @@ private:
}
}
- void iterateMult(Direction dir) {
- if (dir != BACKWARD) {
+ void iterateMult() {
+ if (multDirection_ != BACKWARD) {
for (size_t i=0; i<localProblems_.size(); i++)
iterateStep(i);
}
- if (dir != Direction::FORWARD)
+ if (multDirection_ != Direction::FORWARD)
for (size_t i=localProblems_.size()-1; i>=0 && i<localProblems_.size(); i--)
iterateStep(i);
}
diff --git a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
index 87204de..106edb1 100644
--- a/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh
@@ -32,15 +32,19 @@ private:
using GridType = typename GridView::Grid ;
using BitVector = typename Dune::BitSetVector<1>;
- const Dune::ParameterTree& parset_;
-
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;
- const Mode mode_;
+ Mode mode_;
+
+ using MultDirection = typename LevelPatchPreconditionerType::Direction;
+ MultDirection multDirection_;
+
int itCounter_;
@@ -60,19 +64,13 @@ public:
const std::vector<std::shared_ptr<LocalOperatorAssembler>>& localAssemblers,
const std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>>& localInterfaceAssemblers,
const Dune::ParameterTree& parset) :
- parset_(parset),
interfaceNetwork_(interfaceNetwork),
activeLevels_(activeLevels),
localAssemblers_(localAssemblers),
- localInterfaceAssemblers_(localInterfaceAssemblers)
+ localInterfaceAssemblers_(localInterfaceAssemblers),
+ parset_(parset)
{
- const auto mode = parset_.get<std::string>("mode");
- if (mode == "ADDITIVE") {
- mode_ = Mode::ADDITIVE;
- } else {
- mode_ = Mode::MULTIPLICATIVE;
- }
-
+ 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!");
@@ -157,6 +155,28 @@ public:
}
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() {
//*(this->x_) = 0;
@@ -207,41 +227,13 @@ private:
void iterateMult() {
- *(this->x_) = 0;
-
- VectorType x;
-
- for (int i=(levelPatchPreconditioners_.size()-1); i>=0; i--) {
- VectorType updatedRhs(*(this->rhs_));
- this->mat_->mmv(*(this->x_), updatedRhs);
-
- //print(updatedRhs, "globalRhs: ");
-
- mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
- mgTransfer_[i]->restrict(updatedRhs, levelRhs_[i]);
- //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], levelRhs_[i]);
-
- levelPatchPreconditioners_[i]->iterate();
- const VectorType& it = levelPatchPreconditioners_[i]->getSol();
-
- mgTransfer_[i]->prolong(it, x);
-
- //writeToVTK(*allFaultBasis_, x, "/storage/mi/podlesny/data/faultnetworks/fineIterates/multilevel", "exactvertexdata_step_"+std::to_string(itCounter_));
-
-
- *(this->x_) += x;
- }
-
- VectorType updatedCoarseRhs(*(this->rhs_));
- this->mat_->mmv(*(this->x_), updatedCoarseRhs);
+ 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]);
- mgTransfer_[j]->restrict(updatedCoarseRhs, levelRhs_[j]);
//print(levelRhs_[j], "localCoarseRhs: ");
//writeToVTK(coarseGlobalPreconditioner_->basis(), levelRhs_[j], "/storage/mi/podlesny/data/faultnetworks/rhs/coarse", "exactvertexdata_step_"+std::to_string(itCounter_));
@@ -250,15 +242,38 @@ private:
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_));
- *(this->x_) += x;
+ 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]);
+
+ //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], levelRhs_[i]);
+
+ 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") {
@@ -284,6 +299,8 @@ private:
if (activeLevels_[i][0]) {
levelPatchPreconditioners_.push_back(std::make_shared<CellPreconditioner>(interfaceNetwork_.levelInterfaceNetwork(i), *localAssemblers_[i], localInterfaceAssemblers_[i], mode_));
+ levelPatchPreconditioners_.back()->setBoundaryMode(LevelPatchPreconditionerType::BoundaryMode::homogeneous);
+
maxLevel = i;
}
}
@@ -294,6 +311,8 @@ private:
auto setupSupportPreconditioner() {
using SupportPreconditioner = SupportPatchPreconditioner<BasisType, LocalOperatorAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+ const auto patchDepth = parset_.get<size_t>("patchDepth");
+
int maxLevel = 0;
bool skip = true;
@@ -311,6 +330,8 @@ private:
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;
diff --git a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
index ea23ccd..7acb2c1 100644
--- a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
@@ -61,7 +61,8 @@ public:
rhs.resize(this->matrix_.N());
rhs = 0;
- this->localProblems_[i] = std::make_shared<decltype(*(this->localProblems_[i]))>(this->matrix_, rhs, localToGlobal, boundaryDofs);
+ using LocalProblemType = typename std::remove_reference<decltype(*(this->localProblems_[i]))>::type;
+ 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 cf3fd4e..1dd5701 100644
--- a/src/cantorfaultnetworks/cantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/cantorfaultnetwork.cc
@@ -306,17 +306,16 @@ int main(int argc, char** argv) { try
}
}
- typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner;
- typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner;
- MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;
+ const auto& preconditionerParset = parameterSet.sub("preconditioner");
+ using MultilevelPatchPreconditioner = MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
- MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
- for (size_t i=0; i<preconditioner.size()-1; i++) {
- LevelPatchPreconditioner& levelFaultPreconditioner = *preconditioner.levelPatchPreconditioner(i);
+ MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerParset);
+ /*for (size_t i=0; i<preconditioner.size()-1; i++) {
+ auto& levelFaultPreconditioner = *preconditioner.levelPatchPreconditioner(i);
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
- }
+ }*/
preconditioner.build();
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
index cc13b28..4ecac4b 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
@@ -109,7 +109,6 @@ int main(int argc, char** argv) { try
const size_t maxCantorLevel = problemParameters.get<size_t>("maxCantorLevel");
const double c = problemParameters.get<double>("penaltyFactor");
- const int patchDepth = problemParameters.get<int>("patchDepth");
const int maxIterations = problemParameters.get<int>("maxIterations");
const double solverTolerance = problemParameters.get<double>("tol");
@@ -117,7 +116,7 @@ int main(int argc, char** argv) { try
const int fineGridN = std::pow(4, fineCantorLevel);
- std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minCantorLevel) + "_" + std::to_string(fineCantorLevel) + "_" + std::to_string(maxCantorLevel) + "_mult.log");
+ std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minCantorLevel) + "_" + std::to_string(fineCantorLevel) + "_" + std::to_string(maxCantorLevel) + ".log");
std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
index 354f92a..1f8e671 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset
@@ -2,9 +2,9 @@ path = ../data/
resultPath = ../cantorfaultnetworks/results/sparse/
[preconditioner]
-patch = CELL # CELL , SUPPORT
+patch = SUPPORT # CELL , SUPPORT
mode = ADDITIVE # ADDITIVE , MULTIPLICATIVE
-multMode = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
+multDirection = SYMMETRIC # SYMMETRIC , FORWARD , BACKWARD
patchDepth = 1
###########################################
diff --git a/src/geofaultnetworks/geofaultnetwork.cc b/src/geofaultnetworks/geofaultnetwork.cc
index 1cd45f9..ddad719 100644
--- a/src/geofaultnetworks/geofaultnetwork.cc
+++ b/src/geofaultnetworks/geofaultnetwork.cc
@@ -340,17 +340,16 @@ int main(int argc, char** argv) { try
}
}
- typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner;
- typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner;
- MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;
+ const auto& preconditionerParset = parameterSet.sub("preconditioner");
+ using MultilevelPatchPreconditioner = MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
- MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
- for (size_t i=0; i<preconditioner.size()-1; i++) {
+ MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerParset);
+ /*for (size_t i=0; i<preconditioner.size()-1; i++) {
LevelPatchPreconditioner& levelFaultPreconditioner = *preconditioner.levelPatchPreconditioner(i);
levelFaultPreconditioner.setPatchDepth(patchDepth);
levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous);
- }
+ }*/
preconditioner.build();
std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl;
--
GitLab