From fbd31be585f14394d18d8e5ce6a1ceeb81745a3e Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Tue, 9 Jul 2019 17:34:01 +0200
Subject: [PATCH] .
---
src/CMakeLists.txt | 2 +-
src/data-structures/program_state.hh | 20 +---
src/multi-body-problem.cc | 2 +-
.../levelpatchpreconditioner.hh | 29 ++++--
.../multilevelpatchpreconditioner.hh | 94 +++++++++++--------
.../preconditioners/nbodycontacttransfer.cc | 12 +--
src/spatial-solving/solverfactory.cc | 14 ++-
src/spatial-solving/solverfactory.hh | 3 +
src/spatial-solving/solverfactory_tmpl.cc | 12 ++-
9 files changed, 104 insertions(+), 84 deletions(-)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 26276d41..e7cf62d6 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -47,7 +47,7 @@ set(MSW_SOURCE_FILES
multi-body-problem-data/grid/mygrids.cc
multi-body-problem-data/grid/simplexmanager.cc
multi-body-problem.cc
- spatial-solving/solverfactory.cc
+ #spatial-solving/solverfactory.cc
spatial-solving/fixedpointiterator.cc
#spatial-solving/solverfactory_old.cc
time-stepping/adaptivetimestepper.cc
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index 22a2ddfe..e9db4feb 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -13,7 +13,6 @@
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/cgstep.hh>
-#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include <dune/tectonic/bodydata.hh>
@@ -97,24 +96,11 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
using Matrix = typename ContactNetwork::Matrix;
const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
- // set patch preconditioner for linear correction in TNNMG method
- using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
- using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;
-
const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
- auto gsStep = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
- PatchSolver patchSolver(gsStep,
- preconditionerParset.get<size_t>("maximumIterations"),
- preconditionerParset.get<double>("tolerance"),
- nullptr,
- preconditionerParset.get<Solver::VerbosityMode>("verbosity"),
- false); // absolute error
-
+ using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
- Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<MPPMode>("mode"));
- preconditioner.setPatchSolver(std::forward<PatchSolver>(patchSolver));
- preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+ Preconditioner preconditioner(preconditionerParset, contactNetwork, activeLevels);
preconditioner.build();
using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
@@ -124,7 +110,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
linearSolverStep.setPreconditioner(preconditioner);
EnergyNorm<Matrix, Vector> energyNorm(linearSolverStep);
- LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), energyNorm, Solver::QUIET);
+ LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), energyNorm, Solver::QUIET);
// Solving a linear problem with a multigrid solver
auto const solveLinearProblem = [&](
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 6682374d..4215f278 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -268,7 +268,7 @@ int main(int argc, char *argv[]) {
} */
std::cout << "Done! :)" << std::endl;
- return 1;
+ //return 1;
//printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
diff --git a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
index 0b0f6b59..bd5ff1a2 100644
--- a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
@@ -19,12 +19,27 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
+enum MPPMode {additive, multiplicative};
+inline
+std::istream& operator>>(std::istream& lhs, MPPMode& e)
+{
+ std::string s;
+ lhs >> s;
+
+ if (s == "additive" || s == "1")
+ e = MPPMode::additive;
+ else if (s == "multiplicative" || s == "0")
+ e = MPPMode::multiplicative;
+ else
+ lhs.setstate(std::ios_base::failbit);
+
+ return lhs;
+}
template <class LevelContactNetwork, class PatchSolver, class MatrixType, class VectorType>
class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
public:
- enum Mode {additive, multiplicative};
static const int dim = LevelContactNetwork::dim;
private:
@@ -35,7 +50,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
using PatchFactory = SupportPatchFactory<LevelContactNetwork>;
using Patch = typename PatchFactory::Patch;
- const Mode mode_;
+ const MPPMode mode_;
const LevelContactNetwork& levelContactNetwork_;
const LevelContactNetwork& fineContactNetwork_;
@@ -53,7 +68,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
// for each coarse patch given by levelContactNetwork: set up local problem, compute local correction
LevelPatchPreconditioner(const LevelContactNetwork& levelContactNetwork,
const LevelContactNetwork& fineContactNetwork,
- const Mode mode = LevelPatchPreconditioner::Mode::additive) :
+ const MPPMode mode = additive) :
mode_(mode),
levelContactNetwork_(levelContactNetwork),
fineContactNetwork_(fineContactNetwork),
@@ -119,8 +134,8 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
}
}
- void setPatchSolver(PatchSolver&& patchSolver) {
- patchSolver_ = Dune::Solvers::wrap_own_share<PatchSolver>(std::forward<PatchSolver>(patchSolver));
+ void setPatchSolver(std::shared_ptr<PatchSolver> patchSolver) {
+ patchSolver_ = patchSolver;
}
void setPatchDepth(const size_t patchDepth = 0) {
@@ -146,14 +161,14 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
x = 0;
auto& step = patchSolver_->getIterationStep();
- step.setProblem(this->mat_, x, this->rhs_);
+ dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(step).setProblem(*this->mat_, x, *this->rhs_);
step.setIgnore(p);
patchSolver_->check();
patchSolver_->preprocess();
patchSolver_->solve();
- *(this->x_) += patchSolver_->getSol();
+ *(this->x_) += step.getSol();
}
}
diff --git a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index 6e4ef3bf..d16c2a65 100644
--- a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -6,20 +6,25 @@
#include <dune/common/timer.hh>
#include <dune/common/fvector.hh>
#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include "nbodycontacttransfer.hh"
#include "levelpatchpreconditioner.hh"
-enum MPPMode {additive, multiplicative};
-
-template <class ContactNetwork, class PatchSolver, class MatrixType, class VectorType>
+template <class ContactNetwork, class MatrixType, class VectorType>
class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
public:
//struct Mode { enum ModeType {additive, multiplicative}; };
private:
+ using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+
using LevelContactNetwork = typename ContactNetwork::LevelContactNetwork;
using LevelPatchPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
using BitVector = typename LinearIterationStep<MatrixType, VectorType>::BitVector;
@@ -38,21 +43,20 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
std::vector<MatrixType> levelMat_;
std::vector<VectorType> levelX_;
std::vector<VectorType> levelRhs_;
- std::vector<BitVector> levelIgnore_;
- std::unique_ptr<PatchSolver> coarseSolver_;
+ std::vector<std::shared_ptr<PatchSolver>> levelSolvers_;
std::vector<BitVector> recompute_;
std::vector<BitVector> truncation_;
std::vector<std::shared_ptr<MGTransfer>> mgTransfer_; //TODO
public:
- MultilevelPatchPreconditioner(const ContactNetwork& contactNetwork,
- const Dune::BitSetVector<1>& activeLevels,
- const MPPMode mode = additive) :
+ MultilevelPatchPreconditioner(const Dune::ParameterTree& parset,
+ const ContactNetwork& contactNetwork,
+ const Dune::BitSetVector<1>& activeLevels) :
contactNetwork_(contactNetwork),
activeLevels_(activeLevels),
- mode_(mode) {
+ mode_(parset.get<MPPMode>("mode")){
assert(activeLevels.size() == contactNetwork_.nLevels());
@@ -70,17 +74,12 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
const int coarsestLevel = maxLevel;
- MPPMode levelMode = additive;
- if (mode_ != additive){
- levelMode = multiplicative;
- }
-
for (size_t i=maxLevel+1; i<activeLevels_.size(); i++) {
if (activeLevels_[i][0]) {
// init local patch preconditioners on each level
- const auto& fineNetwork = contactNetwork_.level(i);
- levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditioner>(*contactNetwork_.level(maxLevel), fineNetwork, levelMode));
-
+ const auto& fineNetwork = *contactNetwork_.level(i);
+ levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditioner>(*contactNetwork_.level(maxLevel), fineNetwork, mode_));
+ levelPatchPreconditioners_.back()->setPatchDepth(parset.get<size_t>("patchDepth"));
maxLevel = i;
}
}
@@ -88,6 +87,30 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
levelX_.resize(size());
levelRhs_.resize(size());
+ // init patch solvers
+ levelSolvers_.resize(size());
+ auto coarseStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+ EnergyNorm<MatrixType, VectorType> coarseEnergyNorm(coarseStep);
+ levelSolvers_[0] = std::make_shared<PatchSolver>(coarseStep,
+ parset.get<size_t>("maximumIterations"),
+ parset.get<double>("tolerance"),
+ coarseEnergyNorm,
+ parset.get<Solver::VerbosityMode>("verbosity"),
+ false); // absolute error
+
+ for (size_t i=1; i<levelSolvers_.size(); i++) {
+ auto gsStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+ EnergyNorm<MatrixType, VectorType> energyNorm(gsStep);
+ levelSolvers_[i] = std::make_shared<PatchSolver>(gsStep,
+ parset.get<size_t>("maximumIterations"),
+ parset.get<double>("tolerance"),
+ energyNorm,
+ parset.get<Solver::VerbosityMode>("verbosity"),
+ false); // absolute error
+ levelPatchPreconditioners_[i]->setPatchSolver(levelSolvers_[i]);
+ }
+
+
// init multigrid transfer
mgTransfer_.resize(levelPatchPreconditioners_.size()-1);
for (size_t i=0; i<mgTransfer_.size(); i++) {
@@ -149,7 +172,10 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
coarseTransfer.restrict(levelRhs_[1], levelRhs_[0]);
coarseTransfer.galerkinRestrict(levelMat_[1], levelMat_[0]);
- coarseSolver_.getIterationStep().setProblem(levelMat_[0], levelX_[0], levelRhs_[0]);
+ dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(levelSolvers_[0]->getIterationStep()).setProblem(levelMat_[0], levelX_[0], levelRhs_[0]);
+
+ // EnergyNorm<MatrixType, VectorType> coarseErrorNorm(levelSolvers_[0]->getIterationStep());
+ // levelSolvers_[0]->setErrorNorm(coarseErrorNorm);
}
@@ -166,24 +192,20 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
VectorType x;
// solve coarse global problem
- auto& coarseIgnore = levelIgnore_.back();
- mgTransfer_.back()->restrict(this->ignore(), coarseIgnore);
- coarseSolver_.getIterationStep().setIgnore(coarseIgnore);
+ auto& coarseSolver = levelSolvers_[0];
+ coarseSolver->getIterationStep().setIgnore(truncation_[0]);
- coarseSolver_.check();
- coarseSolver_.preprocess();
- coarseSolver_.solve();
+ coarseSolver->check();
+ coarseSolver->preprocess();
+ coarseSolver->solve();
- mgTransfer_.back()->prolong(coarseSolver_.getSol(), x);
+ mgTransfer_[0]->prolong(coarseSolver->getIterationStep().getSol(), x);
- for (size_t i=levelPatchPreconditioners_.size()-2; i>=0; i--) {
+ for (size_t i=1; i<levelPatchPreconditioners_.size(); i++) {
const auto& transfer = *mgTransfer_[i];
auto& preconditioner = *levelPatchPreconditioners_[i];
- auto& _ignore = levelIgnore_[i];
- transfer.restrict(this->ignore(), _ignore);
- preconditioner.setIgnore(_ignore);
-
+ preconditioner.setIgnore(truncation_[i]);
preconditioner.iterate();
x += preconditioner.getSol();
@@ -238,21 +260,13 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
}
void build() {
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ for (size_t i=1; i<levelPatchPreconditioners_.size(); i++) {
levelPatchPreconditioners_[i]->build();
}
}
- void setPatchSolver(PatchSolver&& patchSolver) { //TODO create copies of patch solver for each level
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- levelPatchPreconditioners_[i]->setPatchSolver(std::forward<PatchSolver>(patchSolver));
- }
-
- coarseSolver_ = Dune::Solvers::wrap_own_share<PatchSolver>(std::forward<PatchSolver>(patchSolver));
- }
-
void setPatchDepth(const size_t patchDepth = 0) {
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+ for (size_t i=1; i<levelPatchPreconditioners_.size(); i++) {
levelPatchPreconditioners_[i]->setPatchDepth(patchDepth);
}
}
diff --git a/src/spatial-solving/preconditioners/nbodycontacttransfer.cc b/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
index 8fba1420..1bb9ce37 100644
--- a/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
+++ b/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
@@ -16,8 +16,8 @@ void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwor
const size_t nBodies = contactNetwork.nBodies();
const size_t nCouplings = contactNetwork.nCouplings();
- const auto& coarseContactNetwork = contactNetwork.level(coarseLevel);
- const auto& fineContactNetwork = contactNetwork.level(fineLevel);
+ const auto& coarseContactNetwork = *contactNetwork.level(coarseLevel);
+ const auto& fineContactNetwork = *contactNetwork.level(fineLevel);
std::vector<size_t> maxL(nBodies), cL(nBodies), fL(nBodies);
@@ -71,13 +71,13 @@ void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwor
*/
- if (fineLevel == contactNetwork.size()-1) {
+ if (fineLevel == contactNetwork.nLevels()-1) {
const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
const auto& contactCouplings = nBodyAssembler.getContactCouplings();
- const std::vector<const MatrixType*> mortarTransferOperators(nBodyAssembler.nCouplings());
- const std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nBodyAssembler.nCouplings());
- const std::vector<std::array<int,2> > gridIdx(nBodyAssembler.nCouplings());
+ std::vector<const MatrixType*> mortarTransferOperators(nBodyAssembler.nCouplings());
+ std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nBodyAssembler.nCouplings());
+ std::vector<std::array<int,2> > gridIdx(nBodyAssembler.nCouplings());
for (size_t i=0; i<nBodyAssembler.nCouplings(); i++) {
mortarTransferOperators[i] = &contactCouplings[i]->mortarLagrangeMatrix();
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 7b6875ce..77c31c1c 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -1,13 +1,11 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
+//#ifdef HAVE_CONFIG_H
+//#include "config.h"
+//#endif
#include <dune/solvers/common/wrapownshare.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include <dune/solvers/solvers/umfpacksolver.hh>
-#include "solverfactory.hh"
-
#include "../utils/debugutils.hh"
template <class Functional, class BitVector>
@@ -21,8 +19,8 @@ SolverFactory<Functional, BitVector>::SolverFactory(
nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
- auto linearSolver = Dune::Solvers::wrap_own_share<LinearSolver>(std::forward<LinearSolver>(linearSolver));
- tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver, DefectProjection(), LineSearchSolver());
+ auto linearSolver_ptr = Dune::Solvers::wrap_own_share<std::decay_t<LinearSolver>>(std::forward<LinearSolver>(linearSolver));
+ tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, DefectProjection(), LineSearchSolver());
tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
tnnmgStep_->setIgnore(ignoreNodes);
}
@@ -40,4 +38,4 @@ auto SolverFactory<Functional, BitVector>::step()
return tnnmgStep_;
}
-#include "solverfactory_tmpl.cc"
+//#include "solverfactory_tmpl.cc"
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index 05235798..48e4a63c 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -55,4 +55,7 @@ class SolverFactory {
// TNNMGStep
std::shared_ptr<Step> tnnmgStep_;
};
+
+#include "solverfactory.cc"
+
#endif
diff --git a/src/spatial-solving/solverfactory_tmpl.cc b/src/spatial-solving/solverfactory_tmpl.cc
index 8e3b2a15..488c5191 100644
--- a/src/spatial-solving/solverfactory_tmpl.cc
+++ b/src/spatial-solving/solverfactory_tmpl.cc
@@ -19,14 +19,18 @@ using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
using MyLinearSolver = Dune::Solvers::LoopSolver<Vector, BitVector>;
+using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
template class SolverFactory<MyFunctional, BitVector>;
-template SolverFactory<MyFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
+
+/*template<> SolverFactory<MyFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
MyFunctional&,
MyLinearSolver&&,
- const BitVector&);
+ const BitVector&);*/
template class SolverFactory<MyZeroFunctional, BitVector>;
-template SolverFactory<MyZeroFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
+
+/*template<>
+template<> SolverFactory<MyZeroFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
MyZeroFunctional&,
MyLinearSolver&&,
- const BitVector&);
+ const BitVector&);*/
--
GitLab