Use TNNMG solver

dune/tectonic/myblockproblem.hh

@@ -23,10 +23,28 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
   typedef typename MyConvexProblemType::LocalMatrixType LocalMatrixType;
 
   int static const block_size = MyConvexProblemType::block_size;
+  int static const coarse_block_size = block_size;
 
   /** \brief Solves one local system using a modified gradient method */
   class IterateObject;
 
+  struct Linearization {
+    static const int block_size = coarse_block_size;
+
+    typedef typename MyBlockProblem<
+        MyConvexProblemTypeTEMPLATE>::LocalMatrixType LocalMatrixType;
+    typedef Dune::BCRSMatrix<typename Linearization::LocalMatrixType>
+    MatrixType;
+    typedef Dune::BlockVector<
+        Dune::FieldVector<double, Linearization::block_size>> VectorType;
+    typedef Dune::BitSetVector<Linearization::block_size> BitVectorType;
+    typename Linearization::MatrixType A;
+    typename Linearization::VectorType b;
+    typename Linearization::BitVectorType ignore;
+
+    Dune::BitSetVector<Linearization::block_size> truncation;
+  };
+
   MyBlockProblem(Dune::ParameterTree const &parset,
                  MyConvexProblemType &problem)
       : parset(parset), problem(problem) {
@@ -38,6 +56,110 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
         0); // safety: acceptance factor for inexact minimization
   }
 
+  std::string virtual getOutput(bool header = false) const {
+    // TODO: implement (not urgent)
+    return std::string();
+  }
+
+  void projectCoarseCorrection(VectorType const &u,
+                               typename Linearization::VectorType const &v,
+                               VectorType &projected_v,
+                               Linearization const &linearization) const {
+    // TODO: implement (not urgent)
+    projected_v = v;
+  };
+
+  double computeDampingParameter(VectorType const &u,
+                                 VectorType const &projected_v) const {
+    // TODO: implement
+    return 1.0;
+  };
+
+  void assembleTruncate(VectorType const &u, Linearization &linearization,
+                        Dune::BitSetVector<block_size> const &ignore) const {
+    // we can just copy the ignore information
+    linearization.ignore = ignore;
+
+    // determine truncation pattern
+    linearization.truncation.resize(u.size());
+    linearization.truncation.unsetAll();
+    Interval<double> ab;
+    for (int i = 0; i < u.size(); ++i) {
+      for (int j = 0; j < block_size; ++j) {
+        // problem.phi.smoothDomain(i, u[i][j], regularity_tol, ab, j);
+
+        // // if (phi is not regular at u) or (u is not in smooth domain) or
+        // (component should be ignored)
+        // // truncate this component
+        // if ((problem.phi.regularity(i, u[i][j], j) > regularity_tol)
+        //     or (u[i][j] < ab[0])
+        //     or (u[i][j] > ab[1])
+        //     or ignore[i][j])
+        //   linearization.truncation[i][j] = true;
+      }
+    }
+
+    // construct sparsity pattern for linearization
+    Dune::MatrixIndexSet indices(problem.A.N(), problem.A.M());
+    indices.import(problem.A);
+    // problem.phi.addHessianIndices(indices);
+
+    // construct matrix from pattern and initialize it
+    indices.exportIdx(linearization.A);
+    linearization.A = 0.0;
+
+    // compute quadratic part of hessian (linearization.A += problem.A)
+    for (int i = 0; i < problem.A.N(); ++i) {
+      typename MatrixType::row_type::ConstIterator it = problem.A[i].begin();
+      typename MatrixType::row_type::ConstIterator end = problem.A[i].end();
+      for (; it != end; ++it)
+        StaticMatrix::axpy(linearization.A[i][it.index()], 1.0, *it);
+    }
+
+    // compute nonlinearity part of hessian
+    problem.phi.addHessian(u, linearization.A);
+
+    // compute quadratic part of gradient
+    linearization.b.resize(u.size());
+    linearization.b = 0.0;
+    problem.A.template umv<VectorType, VectorType>(u, linearization.b);
+    linearization.b -= problem.f;
+
+    // compute nonlinearity part of gradient
+    problem.phi.addGradient(u, linearization.b);
+
+    // -grad is needed for Newton step
+    linearization.b *= -1.0;
+
+    // apply truncation to system
+    typename Linearization::MatrixType::row_type::Iterator it;
+    typename Linearization::MatrixType::row_type::Iterator end;
+    for (int row = 0; row < linearization.A.N(); ++row) {
+      it = linearization.A[row].begin();
+      end = linearization.A[row].end();
+      for (; it != end; ++it) {
+        int col = it.index();
+        for (int i = 0; i < (*it).N(); ++i) {
+          typename Linearization::MatrixType::block_type::row_type::Iterator
+          blockIt = (*it)[i].begin();
+          typename Linearization::MatrixType::block_type::row_type::Iterator
+          blockEnd = (*it)[i].end();
+          for (; blockIt != blockEnd; ++blockIt)
+            if (linearization.truncation[row][i] or linearization
+                    .truncation[col][blockIt.index()])
+              (*blockIt) = 0.0;
+        }
+      }
+
+      for (int j = 0; j < block_size; ++j)
+        if (linearization.truncation[row][j])
+          linearization.b[row][j] = 0.0;
+    }
+
+    // for(int j=0; j<block_size; ++j)
+    //   outStream << std::setw(9) << linearization.truncation.countmasked(j);
+  }
+
   /** \brief Constructs and returns an iterate object */
   IterateObject getIterateObject() {
     return IterateObject(parset, bisection, problem);

src/one-body-sample.cc

@@ -49,6 +49,11 @@
 #include <dune/tectonic/myconvexproblem.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
 int const dim = 2;
 
 template <class GridView>
@@ -242,17 +247,20 @@ int main(int argc, char *argv[]) {
     u1 = 0;
     VectorType u2 = u1;
     VectorType u3 = u1;
+    VectorType u4 = u1;
 
     VectorType u1_diff_new(grid.size(grid.maxLevel(), dim));
     u1_diff_new = 0;
     VectorType u2_diff_new = u1_diff_new;
     VectorType u3_diff_new = u1_diff_new;
+    VectorType u4_diff_new = u1_diff_new;
 
     CellVectorType vonMisesStress;
 
     VectorType b1;
     VectorType b2;
     VectorType b3;
+    VectorType b4;
 
     auto nodalIntegrals =
         Dune::make_shared<Dune::BlockVector<Dune::FieldVector<double, 1>>>();
@@ -272,17 +280,19 @@ int main(int argc, char *argv[]) {
           leafView, p1Basis, neumannNodes, b1, run);
       b2 = b1;
       b3 = b1;
+      b4 = b1;
 
       // b -= linear update
       stiffnessMatrix.mmv(u1, b1);
       stiffnessMatrix.mmv(u2, b2);
       stiffnessMatrix.mmv(u3, b3);
+      stiffnessMatrix.mmv(u4, b4);
 
       if (parset.get<bool>("useNonlinearGS")) {
         MyConvexProblemType myConvexProblem(
             stiffnessMatrix, *myGlobalNonlinearity, b1, u1_diff_new);
-        MyBlockProblemType myBlockProblem(parset, myConvexProblem);
-        nonlinearGSStep.setProblem(u1_diff_new, myBlockProblem);
+        auto myBlockProblem = new MyBlockProblemType(parset, myConvexProblem);
+        nonlinearGSStep.setProblem(u1_diff_new, *myBlockProblem);
 
         LoopSolver<VectorType> solver(&nonlinearGSStep, solver_maxIterations,
                                       solver_tolerance, &energyNorm, verbosity);
@@ -291,6 +301,72 @@ int main(int argc, char *argv[]) {
 
       u1 += u1_diff_new;
 
+      if (parset.get<bool>("useTNNMG")) {
+        MyConvexProblemType myConvexProblem(
+            stiffnessMatrix, *myGlobalNonlinearity, b4, u4_diff_new);
+        auto myBlockProblem = new MyBlockProblemType(parset, myConvexProblem);
+
+        // {{{ Linear Solver;
+        TruncatedBlockGSStep<OperatorType, VectorType> *linearBaseSolverStep =
+            new TruncatedBlockGSStep<OperatorType, VectorType>;
+
+        EnergyNorm<OperatorType, VectorType> *baseEnergyNorm =
+            new EnergyNorm<OperatorType, VectorType>(*linearBaseSolverStep);
+        LoopSolver<VectorType> *linearBaseSolver = new LoopSolver<VectorType>(
+            linearBaseSolverStep, solver_maxIterations, solver_tolerance,
+            baseEnergyNorm, Solver::QUIET);
+        // }}}
+
+        // {{{ Smoothers
+        TruncatedBlockGSStep<OperatorType, VectorType> *linearPresmoother =
+            new TruncatedBlockGSStep<OperatorType, VectorType>;
+        TruncatedBlockGSStep<OperatorType, VectorType> *linearPostsmoother =
+            new TruncatedBlockGSStep<OperatorType, VectorType>;
+        // }}}
+
+        MultigridStep<OperatorType, VectorType> *linearIterationStep =
+            new MultigridStep<OperatorType, VectorType>;
+        linearIterationStep->setNumberOfLevels(levels);
+        linearIterationStep->setMGType(1, 3, 3); // 1/3/3 (mu/nu1/nu2)
+        linearIterationStep->basesolver_ = linearBaseSolver;
+        linearIterationStep->setSmoother(linearPresmoother, linearPostsmoother);
+
+        // {{{ Transfer operators
+        linearIterationStep->mgTransfer_.resize(levels - 1);
+        for (size_t i = 0; i < linearIterationStep->mgTransfer_.size(); i++) {
+          CompressedMultigridTransfer<VectorType> *newTransferOp =
+              new CompressedMultigridTransfer<VectorType>;
+          newTransferOp->setup(grid, i, i + 1);
+          linearIterationStep->mgTransfer_[i] = newTransferOp;
+        }
+        // }}}
+
+        // // Now the actual nonlinear solver
+        typedef MyBlockProblemType TNNMGProblemType;
+        typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType;
+        typedef TruncatedNonsmoothNewtonMultigrid<
+            TNNMGProblemType, NonlinearSmootherType> TNNMGStepType;
+
+        TNNMGProblemType *tnnmgProblem = myBlockProblem;
+        NonlinearSmootherType *nonlinearSmoother = new NonlinearSmootherType;
+
+        auto multigridStep =
+            new TNNMGStepType(*linearIterationStep, *nonlinearSmoother);
+        multigridStep->setProblem(u4_diff_new, *tnnmgProblem);
+        multigridStep->setSmoothingSteps(3, 1, 3); // 3/1/3 (pre/linear/post)
+        multigridStep->ignoreNodes_ = &ignoreNodes;
+
+        auto energyNorm =
+            new EnergyNorm<OperatorType, VectorType>(stiffnessMatrix);
+
+        LoopSolver<VectorType> overallSolver(
+            multigridStep, solver_maxIterations, solver_tolerance, energyNorm,
+            verbosity);
+        overallSolver.solve();
+      }
+
+      u4 += u4_diff_new;
+
       { // Compute von Mises stress and write everything to a file
         auto *displacement =
             new BasisGridFunction<P1Basis, VectorType>(p1Basis, u1);
@@ -351,12 +427,19 @@ int main(int argc, char *argv[]) {
     std::cout << "sup |u1 - u3| = " << diff3.infinity_norm() << ", "
               << "|u1 - u3| = " << diff3.two_norm() << std::endl;
 
+    VectorType diff4 = u4;
+    diff4 -= u1;
+    std::cout << "sup |u1 - u4| = " << diff4.infinity_norm() << ", "
+              << "|u1 - u4| = " << diff4.two_norm() << std::endl;
+
     // Print displacement on frictional boundary
     for (size_t i = 0; i < frictionalNodes.size(); ++i)
       if (frictionalNodes[i][0])
         std::cout << boost::format("u1[%02d] = %+3e, %|40t|u2[%02d] = %+3e "
-                                   "%|80t|u3[%02d] = %+3e") %
-                         i % u1[i] % i % u2[i] % i % u3[i] << std::endl;
+                                   "%|80t|u3[%02d] = %+3e %|120t|u4[%02d] = "
+                                   "%+3e") %
+                         i % u1[i] % i % u2[i] % i % u3[i] % i %
+                         u4[i] << std::endl;
   }
   catch (Dune::Exception &e) {
     Dune::derr << "Dune reported error: " << e << std::endl;

src/one-body-sample.parset

@@ -6,6 +6,7 @@ verbose = false
 useNonlinearGS = true
 useGS = false
 useTruncatedGS = false
+useTNNMG = true
 
 [grid]
 refinements = 5