Construct TNNMG solver just once

src/one-body-sample.cc

@@ -269,6 +269,54 @@ int main(int argc, char *argv[]) {
     assemble_nonlinearity(grid.size(grid.maxLevel(), dim), parset,
                           myGlobalNonlinearity, nodalIntegrals);
 
+    // {{{ Set up TNNMG solver
+    auto linearBaseSolverStep =
+        new TruncatedBlockGSStep<OperatorType, VectorType>;
+    auto const baseEnergyNorm =
+        new EnergyNorm<OperatorType, VectorType>(*linearBaseSolverStep);
+    auto linearBaseSolver = new LoopSolver<VectorType>(
+        linearBaseSolverStep, solver_maxIterations, solver_tolerance,
+        baseEnergyNorm, Solver::QUIET);
+
+    // {{{ Smoothers
+    auto linearPresmoother = new TruncatedBlockGSStep<OperatorType, VectorType>;
+    auto linearPostsmoother =
+        new TruncatedBlockGSStep<OperatorType, VectorType>;
+    // }}}
+
+    auto linearIterationStep = new MultigridStep<OperatorType, VectorType>;
+    linearIterationStep->setNumberOfLevels(levels);
+    linearIterationStep->setMGType(parset.get<int>("solver.tnnmg.linear.mu"),
+                                   parset.get<int>("solver.tnnmg.linear.nu1"),
+                                   parset.get<int>("solver.tnnmg.linear.nu2"));
+    linearIterationStep->basesolver_ = linearBaseSolver;
+    linearIterationStep->setSmoother(linearPresmoother, linearPostsmoother);
+
+    // {{{ Transfer operators
+    std::vector<CompressedMultigridTransfer<VectorType> *> transferOperators;
+    transferOperators.resize(refinements);
+    for (auto &x : transferOperators)
+      x = new CompressedMultigridTransfer<VectorType>;
+
+    TransferOperatorAssembler<GridType> transferOperatorAssembler(grid);
+    transferOperatorAssembler.assembleOperatorPointerHierarchy(
+        transferOperators);
+    linearIterationStep->setTransferOperators(transferOperators);
+    // }}}
+
+    typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType;
+    typedef TruncatedNonsmoothNewtonMultigrid<
+        MyBlockProblemType, NonlinearSmootherType> TNNMGStepType;
+
+    auto nonlinearSmoother = new NonlinearSmootherType;
+    auto multigridStep =
+        new TNNMGStepType(*linearIterationStep, *nonlinearSmoother);
+    multigridStep->setSmoothingSteps(parset.get<int>("solver.tnnmg.main.nu1"),
+                                     parset.get<int>("solver.tnnmg.main.mu"),
+                                     parset.get<int>("solver.tnnmg.main.nu2"));
+    multigridStep->ignoreNodes_ = &ignoreNodes;
+    // }}}
+
     for (size_t run = 1; run <= parset.get<size_t>("timesteps"); ++run) {
       if (parset.get<bool>("printProgress")) {
         std::cout << ".";
@@ -303,81 +351,16 @@ int main(int argc, char *argv[]) {
       u1 += u1_diff;
 
       if (parset.get<bool>("useTNNMG")) {
-        // {{{ Linear Solver;
-        auto linearBaseSolverStep =
-            new TruncatedBlockGSStep<OperatorType, VectorType>;
-        auto const baseEnergyNorm =
-            new EnergyNorm<OperatorType, VectorType>(*linearBaseSolverStep);
-        auto linearBaseSolver = new LoopSolver<VectorType>(
-            linearBaseSolverStep, solver_maxIterations, solver_tolerance,
-            baseEnergyNorm, Solver::QUIET);
-        // }}}
-
-        // {{{ Smoothers
-        auto linearPresmoother =
-            new TruncatedBlockGSStep<OperatorType, VectorType>;
-        auto linearPostsmoother =
-            new TruncatedBlockGSStep<OperatorType, VectorType>;
-        // }}}
-
-        auto linearIterationStep = new MultigridStep<OperatorType, VectorType>;
-        linearIterationStep->setNumberOfLevels(levels);
-        linearIterationStep->setMGType(
-            parset.get<int>("solver.tnnmg.linear.mu"),
-            parset.get<int>("solver.tnnmg.linear.nu1"),
-            parset.get<int>("solver.tnnmg.linear.nu2"));
-        linearIterationStep->basesolver_ = linearBaseSolver;
-        linearIterationStep->setSmoother(linearPresmoother, linearPostsmoother);
-
-        // {{{ Transfer operators
-        std::vector<CompressedMultigridTransfer<VectorType> *>
-        transferOperators;
-        transferOperators.resize(refinements);
-        for (auto &x : transferOperators)
-          x = new CompressedMultigridTransfer<VectorType>;
-
-        TransferOperatorAssembler<GridType> transferOperatorAssembler(grid);
-        transferOperatorAssembler.assembleOperatorPointerHierarchy(
-            transferOperators);
-        linearIterationStep->setTransferOperators(transferOperators);
-        // }}}
-
-        typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType;
-        typedef TruncatedNonsmoothNewtonMultigrid<
-            MyBlockProblemType, NonlinearSmootherType> TNNMGStepType;
-
         MyConvexProblemType myConvexProblem(stiffnessMatrix,
                                             *myGlobalNonlinearity, b4, u4_diff);
         auto myBlockProblem = new MyBlockProblemType(parset, myConvexProblem);
-        auto nonlinearSmoother = new NonlinearSmootherType;
-        auto multigridStep =
-            new TNNMGStepType(*linearIterationStep, *nonlinearSmoother);
         multigridStep->setProblem(u4_diff, *myBlockProblem);
-        multigridStep->setSmoothingSteps(
-            parset.get<int>("solver.tnnmg.main.nu1"),
-            parset.get<int>("solver.tnnmg.main.mu"),
-            parset.get<int>("solver.tnnmg.main.nu2"));
-        multigridStep->ignoreNodes_ = &ignoreNodes;
-
-        auto const energyNorm =
-            new EnergyNorm<OperatorType, VectorType>(stiffnessMatrix);
 
         LoopSolver<VectorType> overallSolver(
-            multigridStep, solver_maxIterations, solver_tolerance, energyNorm,
+            multigridStep, solver_maxIterations, solver_tolerance, &energyNorm,
             verbosity);
         overallSolver.solve();
-        delete linearBaseSolverStep;
-        delete baseEnergyNorm;
-        delete linearBaseSolver;
-        delete linearPresmoother;
-        delete linearPostsmoother;
-        for (auto &x : transferOperators)
-          delete x;
-        delete linearIterationStep;
         delete myBlockProblem;
-        delete nonlinearSmoother;
-        delete multigridStep;
-        delete energyNorm;
       }
 
       u4 += u4_diff;
@@ -434,6 +417,19 @@ int main(int argc, char *argv[]) {
     }
     std::cout << std::endl;
 
+    // {{{ Clean up after the TNNMG solver
+    delete linearBaseSolverStep;
+    delete baseEnergyNorm;
+    delete linearBaseSolver;
+    delete linearPresmoother;
+    delete linearPostsmoother;
+    for (auto &x : transferOperators)
+      delete x;
+    delete linearIterationStep;
+    delete nonlinearSmoother;
+    delete multigridStep;
+    // }}}
+
     if (parset.get<bool>("printDifference")) {
       VectorType diff2 = u2;
       diff2 -= u1;