.

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;