diff --git a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index 1b75e51e3f30ef7546010ef630f23eea543850a0..9d6d64e6eaf8391c6fe6e88290c3a76140130ee1 100644
--- a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -24,14 +24,12 @@
template <class ContactNetwork, class MatrixType, class VectorType>
class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
-public:
- //struct Mode { enum ModeType {additive, multiplicative}; };
-
private:
using Base = LinearIterationStep<MatrixType, VectorType>;
- using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
- //using PatchSolver = Dune::Solvers::UMFPackSolver<MatrixType, VectorType>;
+ using PatchSolver = LoopSolver<Vector, BitVector>;
+ using PatchSolverStep = TruncatedBlockGSStep<MatrixType, VectorType>;
+ using CoarseSolver = Dune::Solvers::UMFPackSolver<MatrixType, VectorType>;
using LevelContactNetwork = typename ContactNetwork::LevelContactNetwork;
using LevelPatchPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
@@ -57,10 +55,10 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
std::vector<std::shared_ptr<EnergyNorm<MatrixType, VectorType>>> levelErrorNorms_;
std::vector<std::shared_ptr<LinearIterationStep<MatrixType, VectorType>>> levelItSteps_;
std::vector<std::shared_ptr<PatchSolver>> levelSolvers_;
+ CoarseSolver coarseSolver_;
- std::vector<BitVector> recompute_;
- //std::vector<BitVector> truncation_;
- std::vector<std::shared_ptr<MGTransfer>> mgTransfer_; //TODO
+ //std::vector<BitVector> recompute_;
+ std::vector<std::shared_ptr<MGTransfer>> mgTransfer_;
public:
MultilevelPatchPreconditioner(const Dune::ParameterTree& parset,
@@ -84,8 +82,6 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
}
}
- const int coarsestLevel = maxLevel;
-
for (size_t i=maxLevel+1; i<activeLevels_.size(); i++) {
if (activeLevels_[i][0]) {
// init local patch preconditioners on each level
@@ -104,16 +100,6 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
levelSolvers_.resize(size());
levelItSteps_.resize(size());
levelErrorNorms_.resize(size());
- //auto coarseStep =
- levelItSteps_[0] = std::make_shared<TruncatedBlockGSStep<MatrixType, VectorType>>();//Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
- levelErrorNorms_[0] = std::make_shared<EnergyNorm<MatrixType, VectorType>>(*levelItSteps_[0].get());
- levelSolvers_[0] = std::make_shared<PatchSolver>(*levelItSteps_[0].get(),
- parset.get<size_t>("maximumIterations"),
- parset.get<double>("tolerance"),
- *levelErrorNorms_[0].get(),
- parset.get<Solver::VerbosityMode>("verbosity"),
- false); // absolute error
- //std::make_shared<PatchSolver>();
for (size_t i=1; i<levelSolvers_.size(); i++) {
//auto gsStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
@@ -123,13 +109,10 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
parset.get<size_t>("maximumIterations"),
parset.get<double>("tolerance"),
*levelErrorNorms_[i].get(),
- parset.get<Solver::VerbosityMode>("verbosity"),
- false); // absolute error
- //std::make_shared<PatchSolver>();
+ parset.get<Solver::VerbosityMode>("verbosity"));
levelPatchPreconditioners_[i]->setPatchSolver(levelSolvers_[i]);
}
-
// init multigrid transfer
mgTransfer_.resize(levelPatchPreconditioners_.size()-1);
for (size_t i=0; i<mgTransfer_.size(); i++) {
@@ -137,19 +120,23 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
mgTransfer_[i]->setup(contactNetwork_, levelPatchPreconditioners_[i+1]->level(), levelPatchPreconditioners_[i+1]->fineLevel());
}
- // Remove any recompute filed so that initially the full transferoperator is assembled
+ // from dune/contact/solvers/nsnewtonmgstep.cc
+ /*// Remove any recompute filed so that initially the full transferoperator is assembled
for (size_t i=0; i<mgTransfer_.size(); i++)
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(nullptr);
// specify the subset of entries to be reassembled at each iteration
- /* recompute_.resize(size());
+ recompute_.resize(size());
recompute_[recompute_.size()-1] = contactNetwork_.nBodyAssembler().totalHasObstacle_;
- for (int i=mgTransfer_.size()-1; i>=0; i--) {
- std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
+ for (int i=mgTransfer_.size(); i>=1; i--) {
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i-1])->restrict(recompute_[i], recompute_[i-1]);
// std::dynamic_pointer_cast<DenseMultigridTransfer<VectorType, BitVector, MatrixType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
- std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(&recompute_[i]);
- print(recompute_[i], "recompute: ");
- }*/
+
+ //print(recompute_[i], "recompute: ");
+ }
+
+ for (size_t i=0; i<this->mgTransfer_.size(); i++)
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(&recompute_[i]); */
}
void setMatrix(const MatrixType& mat) override {
@@ -161,32 +148,16 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrictToFathers(ignoreHierarchy_[i+1], ignoreHierarchy_[i]);
}
- /* truncation_.resize(size());
- truncation_[size()-1] = Base::ignore();
- print(truncation_.back(), "truncation: ");
- for (int i=mgTransfer_.size()-1; i>=0; i--) {
- std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(truncation_[i+1], truncation_[i]);
- std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setCriticalBitField(&truncation_[i+1]);
- print(truncation_[i], "truncation: ");
- }*/
-
size_t maxLevel = levelPatchPreconditioners_.size()-1;
levelPatchPreconditioners_[maxLevel]->setMatrix(mat);
- //print(mat, "levelMat_" + std::to_string(maxLevel) + ":");
-
for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
const auto& transfer = *mgTransfer_[i];
- //print(*levelPatchPreconditioners_[i+1]->getMatrix(), "levelMat:");
- //print(levelMat_[i], "coarseMat:");
-
transfer.galerkinRestrictSetOccupation(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
transfer.galerkinRestrict(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
levelPatchPreconditioners_[i]->setMatrix(levelMat_[i]);
- //print(levelMat_[i], "levelMat_" + std::to_string(i) + ":");
-
// Set solution vector sizes for the lower levels
Dune::MatrixVector::resize(levelX_[i], levelMat_[i]);
//Dune::MatrixVector::resize(levelRhs_[i], levelMat_[i]);
@@ -198,58 +169,9 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
coarseTransfer.galerkinRestrictSetOccupation(fineMat, levelMat_[0]);
coarseTransfer.galerkinRestrict(fineMat, levelMat_[0]);
- //print(levelMat_[0], "levelMat_0:");
-
- /*dynamic_cast<Base&>(levelSolvers_[0]->getIterationStep()).setMatrix(levelMat_[0]); */
-
Dune::MatrixVector::resize(levelX_[0], levelMat_[0]);
- //Dune::MatrixVector::resize(levelRhs_[0], levelMat_[0]);
}
- void setTruncation(const BitVector& truncation) {
- std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_.back())->setCriticalBitField(&truncation);
- }
-
- /*void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override {
- LinearIterationStep<MatrixType, VectorType>::setProblem(mat, x, rhs);
-
- truncation_.resize(size());
- truncation_[size()-1] = this->ignore();
- print(truncation_.back(), "truncation: ");
- for (int i=mgTransfer_.size()-1; i>=0; i--) {
- std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(truncation_[i+1], truncation_[i]);
- std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setCriticalBitField(&truncation_[i]);
- print(truncation_[i], "truncation: ");
- }
-
- size_t maxLevel = levelPatchPreconditioners_.size()-1;
- levelX_[maxLevel] = x;
- levelRhs_[maxLevel] = rhs;
- levelPatchPreconditioners_[maxLevel]->setProblem(mat, levelX_[maxLevel], levelRhs_[maxLevel]);
-
- for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
- const auto& transfer = *mgTransfer_[i];
-
- transfer.restrict(levelX_[i+1], levelX_[i]);
- transfer.restrict(levelRhs_[i+1], levelRhs_[i]);
- transfer.galerkinRestrict(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
-
- levelPatchPreconditioners_[i]->setProblem(levelMat_[i], levelX_[i], levelRhs_[i]);
- }
-
- // set coarse global problem
- const auto& coarseTransfer = *mgTransfer_[0];
-
- coarseTransfer.restrict(levelX_[1], levelX_[0]);
- coarseTransfer.restrict(levelRhs_[1], levelRhs_[0]);
- coarseTransfer.galerkinRestrict(levelMat_[1], levelMat_[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);
- } */
-
void iterate() {
size_t maxLevel = levelPatchPreconditioners_.size()-1;
levelX_[maxLevel] = *this->getIterate();
@@ -280,33 +202,15 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
VectorType x;
// solve coarse global problem
- auto& coarseSolver = levelSolvers_[0];
- coarseSolver->getIterationStep().setIgnore(ignoreHierarchy_[0]);
-
- dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(levelSolvers_[0]->getIterationStep()).setProblem(levelMat_[0], levelX_[0], levelRhs_[0]);
-
- coarseSolver->check();
- coarseSolver->preprocess();
- coarseSolver->solve();
-
- mgTransfer_[0]->prolong(levelX_[0], x);
-
-
- /* BitVector emptyIgnore0(levelX_[0].size(), false);
LocalProblem<MatrixType, VectorType> localProblem(levelMat_[0], levelRhs_[0], ignoreHierarchy_[0]);
Vector newR;
localProblem.getLocalRhs(levelX_[0], newR);
- auto& coarseSolver = levelSolvers_[0];
- coarseSolver->setProblem(localProblem.getMat(), levelX_[0], newR);
- coarseSolver->preprocess();
- coarseSolver->solve();
-
- mgTransfer_[0]->prolong(levelX_[0], x);*/
-
- //print(levelX_[0], "MultilevelPatchPreconditioner: levelX_[0]");
- //print(x, "MultilevelPatchPreconditioner: x0");
+ coarseSolver_.setProblem(localProblem.getMat(), levelX_[0], newR);
+ coarseSolver_.preprocess();
+ coarseSolver_.solve();
+ mgTransfer_[0]->prolong(levelX_[0], x);
for (size_t i=1; i<levelPatchPreconditioners_.size()-1; i++) {
const auto& transfer = *mgTransfer_[i];
@@ -315,18 +219,12 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
preconditioner.setIgnore(ignoreHierarchy_[i]);
preconditioner.apply(levelX_[i], levelRhs_[i]);
- //print(levelX_[i], "MultilevelPatchPreconditioner: levelX_[i] " + std::to_string(i));
-
x += preconditioner.getSol();
VectorType newX;
transfer.prolong(x, newX);
- //print(newX, "MultilevelPatchPreconditioner: newX " + std::to_string(i));
-
x = newX;
-
- //print(x, "MultilevelPatchPreconditioner: x " + std::to_string(i));
}
auto& preconditioner = *levelPatchPreconditioners_.back();
@@ -336,29 +234,15 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
x += preconditioner.getSol();
- //print(x, "MultilevelPatchPreconditioner: xmax");
-
-
- /* BitVector emptyIgnore(this->x_->size(), false);
- LocalProblem<MatrixType, VectorType> globalProblem(*this->mat_.get(), *this->rhs_, emptyIgnore);
- globalProblem.getLocalRhs(*this->x_, newR);
- VectorType totalX = *(this->x_);
- auto& globalSolver = levelSolvers_.back();
- globalSolver->setProblem(globalProblem.getMat(), totalX, newR);
- globalSolver->preprocess();
- globalSolver->solve();
-
- print(totalX, "MultilevelPatchPreconditioner: totalX"); */
-
*(this->x_) = x;
-
- //print(*(this->x_), "MultilevelPatchPreconditioner: sol");
}
void iterateMult() {
*(this->x_) = 0;
+ DUNE_THROW(Dune::Exception, "Not implemented!");
+
/*
VectorType x;