From 50340c30dbc217205f93392d64b1bcc6c0368e68 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Sun, 31 May 2020 17:13:53 +0200
Subject: [PATCH] .
---
.../patchfactories/cellpatchfactory.hh | 6 +-
.../patchfactories/supportpatchfactory.hh | 24 +--
.../preconditioners/CMakeLists.txt | 4 +
.../cellpatchpreconditioner.hh | 43 ++++
.../levelpatchpreconditioner.hh | 190 ++++++------------
.../supportpatchpreconditioner.hh | 70 +++++++
.../sparsecantorfaultnetwork.cc | 2 +-
7 files changed, 186 insertions(+), 153 deletions(-)
create mode 100644 dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
create mode 100644 dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
diff --git a/dune/faultnetworks/patchfactories/cellpatchfactory.hh b/dune/faultnetworks/patchfactories/cellpatchfactory.hh
index 530c462..7b2379e 100644
--- a/dune/faultnetworks/patchfactories/cellpatchfactory.hh
+++ b/dune/faultnetworks/patchfactories/cellpatchfactory.hh
@@ -23,9 +23,7 @@ class CellPatchFactory
typedef typename GridType::ctype ctype;
static const int dim = GridType::dimension;
typedef typename GridType::template Codim<0>::Entity Element;
- typedef typename GridType::template Codim<dim>::Entity Vertex;
- typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType, Dune::MCMGElementLayout > ElementMapper;
private:
@@ -42,13 +40,13 @@ class CellPatchFactory
ElementMapper mapper;
private:
- void print(const Dune::BitSetVector<1>& x, std::string message){
+ /*void print(const Dune::BitSetVector<1>& x, std::string message){
std::cout << message << std::endl;
for (size_t i=0; i<x.size(); i++) {
std::cout << x[i][0] << " ";
}
std::cout << std::endl << std::endl;
- }
+ }*/
bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) {
return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
diff --git a/dune/faultnetworks/patchfactories/supportpatchfactory.hh b/dune/faultnetworks/patchfactories/supportpatchfactory.hh
index 20754c0..07a297d 100644
--- a/dune/faultnetworks/patchfactories/supportpatchfactory.hh
+++ b/dune/faultnetworks/patchfactories/supportpatchfactory.hh
@@ -1,7 +1,8 @@
#ifndef SUPPORT_PATCH_FACTORY_HH
#define SUPPORT_PATCH_FACTORY_HH
-#include<queue>
+#include <queue>
+#include <set>
#include <dune/common/bitsetvector.hh>
#include <dune/common/fvector.hh>
@@ -227,27 +228,6 @@ private:
}
}
}
-
- /* add interior coarse level vertices to boundary dofs
- for(int i=0; i<coarseGeometry.corners(); ++i) {
- const auto& local = fineGeometry.local(coarseGeometry.corner(i));
-
- if (!ref.checkInside(local))
- continue;
-
- const int localBasisSize = fineLFE.localBasis().size();
- std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
- fineLFE.localBasis().evaluateFunction(local, localBasisRep);
-
- for(int k=0; k<localBasisSize; ++k) {
- if (localBasisRep[k]==1) {
- int dofIndex = fineBasis_.index(fineElem, k);
- localBoundaryDofs.insert(dofIndex);
- break;
- }
- }
- } */
-
}
};
#endif
diff --git a/dune/faultnetworks/preconditioners/CMakeLists.txt b/dune/faultnetworks/preconditioners/CMakeLists.txt
index f427498..4063392 100644
--- a/dune/faultnetworks/preconditioners/CMakeLists.txt
+++ b/dune/faultnetworks/preconditioners/CMakeLists.txt
@@ -2,12 +2,16 @@
add_custom_target(faultnetworks_preconditioners_src SOURCES
levelglobalpreconditioner.hh
+ cellpatchpreconditioner.hh
+ supportpatchpreconditioner.hh
levelpatchpreconditioner.hh
multilevelpatchpreconditioner.hh
)
install(FILES
levelglobalpreconditioner.hh
+ cellpatchpreconditioner.hh
+ supportpatchpreconditioner.hh
levelpatchpreconditioner.hh
multilevelpatchpreconditioner.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/preconditioners)
diff --git a/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
new file mode 100644
index 0000000..e938979
--- /dev/null
+++ b/dune/faultnetworks/preconditioners/cellpatchpreconditioner.hh
@@ -0,0 +1,43 @@
+#ifndef CELL_PATCH_PRECONDITIONER_HH
+#define CELL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
+#include <dune/faultnetworks/patchfactories/cellpatchfactory.hh>
+
+template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class CellPatchPreconditioner : public LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType> {
+private:
+ using Base = LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+ using GridView = typename Base::GridView;
+ using Mode = typename Base::Mode;
+
+public:
+ CellPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+ const LocalAssembler& localAssembler,
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+ const Mode mode = Mode::additive) :
+ Base(levelInterfaceNetwork, localAssembler, localInterfaceAssemblers, mode) {}
+
+ void build() {
+ CellPatchFactory<BasisType> patchFactory(this->basis_, this->level_);
+ const auto& localToGlobal = patchFactory.getLocalToGlobal();
+ const auto& boundaryDofs = patchFactory.getBoundaryDofs();
+
+ VectorType rhs;
+ rhs.resize(this->matrix_.N());
+ rhs = 0;
+
+ std::cout << "CellPatchPreconditioner::build() level: " << this->level_ << std::endl;
+
+ 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]);
+ }
+ }
+};
+
+#endif
+
diff --git a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
index bbe464f..60a06a3 100644
--- a/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
+++ b/dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh
@@ -10,7 +10,6 @@
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/faultnetworks/assemblers/globalfaultassembler.hh>
-#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
#include <dune/faultnetworks/localproblem.hh>
#include <dune/faultnetworks/levelinterfacenetwork.hh>
#include <dune/faultnetworks/utils/debugutils.hh>
@@ -18,20 +17,19 @@
#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 Mode {ADDITIVE, MULTIPLICATIVE};
enum BoundaryMode {homogeneous, fromIterate};
+ enum Direction {FORWARD, BACKWARD, SYMMETRIC};
private:
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_;
@@ -46,30 +44,8 @@ private:
MatrixType matrix_;
std::vector<std::shared_ptr<OscLocalProblem<MatrixType, VectorType>> > localProblems_;
-public:
-
- // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
- 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_)
- {
- setPatchDepth();
- setBoundaryMode();
-
- assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
- }
-
- void build() {
+private:
+ void setup() {
// assemble stiffness matrix for entire level including all faults
GlobalFaultAssembler<BasisType, BasisType> globalFaultAssembler(basis_, basis_, levelInterfaceNetwork_);
globalFaultAssembler.assembleOperator(localAssembler_, localInterfaceAssemblers_, matrix_);
@@ -96,70 +72,31 @@ public:
}
row[i] = 1;
}
+ }
- // init vertexInElements
- const int dim = GridType::dimension;
- typedef typename GridType::template Codim<0>::Entity EntityType;
- std::vector<std::vector<EntityType>> vertexInElements;
-
- const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
- vertexInElements.resize(patchLevelGridView.size(dim));
- for (size_t i=0; i<vertexInElements.size(); i++) {
- vertexInElements[i].resize(0);
- }
-
- typedef typename BasisType::LocalFiniteElement FE;
- typedef typename GridView::template Codim <0>::Iterator ElementLevelIterator;
- ElementLevelIterator endElemIt = patchLevelGridView.template end <0>();
- for (ElementLevelIterator it = patchLevelGridView. template begin <0>(); it!=endElemIt; ++it) {
-
- // compute coarseGrid vertexInElements
- const FE& coarseFE = patchLevelBasis_.getLocalFiniteElement(*it);
-
- for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
- int dofIndex = patchLevelBasis_.indexInGridView(*it, i);
- vertexInElements[dofIndex].push_back(*it);
- }
- }
-
- localProblems_.resize(vertexInElements.size());
-
- std::cout << std::endl;
- std::cout << "---------------------------------------------" << std::endl;
- std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
-
- // init local fine level corrections
- Dune::Timer timer;
- timer.reset();
- timer.start();
-
- const int patchLevel = patchLevelBasis_.faultNetwork().level();
- for (size_t i=0; i<vertexInElements.size(); i++) {
- //std::cout << i << std::endl;
- //std::cout << "---------------" << std::endl;
-
- SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, basis_, patchLevel, level_, vertexInElements, i, patchDepth_);
- std::vector<int>& localToGlobal = patchFactory.getLocalToGlobal();
- Dune::BitSetVector<1>& boundaryDofs = patchFactory.getBoundaryDofs();
- VectorType rhs;
- rhs.resize(matrix_.N());
- rhs = 0;
-
- //print(localToGlobal, "localToGlobal: ");
- //print(boundaryDofs, "boundaryDofs: ");
+public:
+ // has to setup localProblems_
+ virtual void build();
- localProblems_[i] = std::make_shared<OscLocalProblem<MatrixType, VectorType>>(matrix_, rhs, localToGlobal, boundaryDofs);
+public:
+ LevelPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+ const LocalAssembler& localAssembler,
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+ const Mode mode = LevelPatchPreconditioner::Mode::additive) :
+ levelInterfaceNetwork_(levelInterfaceNetwork),
+ localAssembler_(localAssembler),
+ localInterfaceAssemblers_(localInterfaceAssemblers),
+ mode_(mode),
+ grid_(levelInterfaceNetwork_.grid()),
+ level_(levelInterfaceNetwork_.level()),
+ basis_(levelInterfaceNetwork_)
+ {
+ setPatchDepth();
+ setBoundaryMode();
- /*
- if ((i+1) % 10 == 0) {
- std::cout << '\r' << std::floor((i+1.0)/patchLevelBasis_.size()*100) << " % done. Elapsed time: " << timer.elapsed() << " seconds. Predicted total setup time: " << timer.elapsed()/(i+1.0)*patchLevelBasis_.size() << " seconds." << std::flush;
- }*/
- }
+ assert(localInterfaceAssemblers_.size() == levelInterfaceNetwork_.size());
- std::cout << std::endl;
- std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
- std::cout << "---------------------------------------------" << std::endl;
- timer.stop();
+ setup();
}
void setPatchDepth(const size_t patchDepth = 0) {
@@ -183,13 +120,30 @@ public:
}
}
- virtual void iterate() {
- if (mode_ == additive)
+ virtual void iterate(Direction dir = SYMMETRIC) {
+ if (mode_ == ADDITIVE)
iterateAdd();
else
- iterateMult();
+ iterateMult(dir);
+ }
+
+ 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 iterateAdd() {
*(this->x_) = 0;
@@ -206,54 +160,38 @@ public:
}
}
- void iterateMult() {
+ void iterateMult(Direction dir) {
+ if (dir != BACKWARD) {
+ for (size_t i=0; i<localProblems_.size(); i++)
+ iterateStep(i);
+ }
+
+ if (dir != 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;
- VectorType it, x;
- for (size_t i=0; i<localProblems_.size(); i++) {
- VectorType updatedRhs(*(this->rhs_));
- matrix_.mmv(*(this->x_), updatedRhs);
+ auto rhs = *(this->rhs_);
- //step(i);
- //print(updatedRhs, "localRhs: ");
- //writeToVTK(basis_, updatedRhs, "/storage/mi/podlesny/data/faultnetworks/rhs/local", "exactvertexdata_step_"+std::to_string(i));
+ VectorType it, x;
if (boundaryMode_ == BoundaryMode::homogeneous)
- localProblems_[i]->updateRhs(updatedRhs);
+ localProblems_[i]->updateRhs(rhs);
else
- localProblems_[i]->updateRhsAndBoundary(updatedRhs, *(this->x_));
+ localProblems_[i]->updateRhsAndBoundary(rhs, *(this->x_));
localProblems_[i]->solve(it);
localProblems_[i]->prolong(it, x);
+ VectorType residual;
+ localProblems_[i]->localResidual(it, residual);
+ localProblems_[i]->updateVector(residual, rhs);
- //print(it, "it: ");
- //print(x, "x: ");
-
- //writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/sol/local", "exactvertexdata_step_"+std::to_string(i));
-
- /*if (i==5) {
- writeToVTK(basis_, x, "/storage/mi/podlesny/data/faultnetworks/iterates/", "exactvertexdata_patchDepth_"+std::to_string(patchDepth_));
- }*/
*(this->x_) += x;
- }
- }
-
- 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();
}
};
diff --git a/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
new file mode 100644
index 0000000..baad90f
--- /dev/null
+++ b/dune/faultnetworks/preconditioners/supportpatchpreconditioner.hh
@@ -0,0 +1,70 @@
+#ifndef SUPPORT_PATCH_PRECONDITIONER_HH
+#define SUPPORT_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/faultnetworks/preconditioners/levelpatchpreconditioner.hh>
+#include <dune/faultnetworks/patchfactories/supportpatchfactory.hh>
+
+template <class BasisType, class LocalAssembler, class LocalInterfaceAssembler, class MatrixType, class VectorType>
+class SupportPatchPreconditioner : public LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType> {
+private:
+ using Base = LevelPatchPreconditioner<BasisType, LocalAssembler, LocalInterfaceAssembler, MatrixType, VectorType>;
+ using GridView = typename Base::GridView;
+ using GridType = typename Base::GridType;
+ using Mode = typename Base::Mode;
+
+ const BasisType& patchLevelBasis_;
+
+public:
+ // for each coarse patch given by patchLevelBasis: set up local problem, compute local correction
+ SupportPatchPreconditioner(const LevelInterfaceNetwork<GridView>& levelInterfaceNetwork,
+ const BasisType& patchLevelBasis,
+ const LocalAssembler& localAssembler,
+ const std::vector<std::shared_ptr<LocalInterfaceAssembler>>& localInterfaceAssemblers,
+ const Mode mode = Mode::additive) :
+ Base(levelInterfaceNetwork, localAssembler, localInterfaceAssemblers, mode),
+ patchLevelBasis_(patchLevelBasis) {}
+
+ void build() {
+ // init vertexInElements
+ const int dim = GridType::dimension;
+ using Element = typename GridType::template Codim<0>::Entity;
+ std::vector<std::vector<Element>> vertexInElements;
+
+ const GridView& patchLevelGridView = patchLevelBasis_.getGridView();
+ vertexInElements.resize(patchLevelGridView.size(dim));
+ for (size_t i=0; i<vertexInElements.size(); i++) {
+ vertexInElements[i].resize(0);
+ }
+
+ for (const auto& elem : elements(patchLevelGridView)) {
+ // compute coarseGrid vertexInElements
+ const auto& coarseFE = patchLevelBasis_.getLocalFiniteElement(elem);
+
+ for (size_t i=0; i<coarseFE.localBasis().size(); ++i) {
+ int dofIndex = patchLevelBasis_.indexInGridView(elem, i);
+ vertexInElements[dofIndex].push_back(elem);
+ }
+ }
+
+ std::cout << "SupportPatchPreconditioner::build() level: " << this->level_ << std::endl;
+
+ this->localProblems_.resize(vertexInElements.size());
+ const int patchLevel = patchLevelBasis_.faultNetwork().level();
+ for (size_t i=0; i<vertexInElements.size(); i++) {
+ SupportPatchFactory<BasisType> patchFactory(patchLevelBasis_, this->basis_, patchLevel, this->level_, vertexInElements, i, this->patchDepth_);
+ const auto& localToGlobal = patchFactory.getLocalToGlobal();
+ const auto& boundaryDofs = patchFactory.getBoundaryDofs();
+
+ VectorType rhs;
+ rhs.resize(this->matrix_.N());
+ rhs = 0;
+
+ this->localProblems_[i] = std::make_shared<decltype(*(this->localProblems_[i]))>(this->matrix_, rhs, localToGlobal, boundaryDofs);
+ }
+ }
+};
+
+#endif
+
diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
index 3441ee8..2dbe0de 100644
--- a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
+++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc
@@ -302,7 +302,7 @@ 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::multiplicative;
+ MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::additive;
MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode);
for (size_t i=0; i<preconditioner.size()-1; i++) {
--
GitLab