Refactor test.

dune/solvers/test/gssteptest.cc

@@ -16,34 +16,6 @@
 
 #include "common.hh"
 
-template <class Vector>
-class Analyser {
-public:
-  Analyser(Vector const &reference,
-           Norm<Vector> const &norm,
-           double tolerance)
-    : reference_(reference), norm_(norm), tolerance_(tolerance)
-  {}
-
-  bool
-  analyse(Vector const &candidate, std::string testCase) {
-    const auto normDiff = norm_.diff(reference_, candidate);
-    const std::string s = Dune::formatString("%*s", -40, testCase.c_str());
-    std::cout << "error for solver \"" << s
-              << "\": " << normDiff << std::endl;
-    if (normDiff > tolerance_) {
-      std::cerr << "Error too large for solver \"" << s << "\"."
-                << std::endl;
-      return false;
-    }
-    return true;
-  }
-private:
-  Vector const &reference_;
-  Norm<Vector> const &norm_;
-  double const tolerance_;
-};
-
 /**  \brief test for the GSStep classes
   *
   *  This test tests if the Dirichlet problem for a Laplace operator
@@ -55,34 +27,24 @@ struct GSTestSuite
   template <class GridType>
   bool check(const GridType& grid)
   {
-    double stepTol    = 1e-12; // termination criterion
-    double solveTol   = 1.3e-5;  // error in the solution
-    int maxIterations = 10000;
-
     bool passed = true;
 
-    using Problem =
-    SymmetricSampleProblem<1, typename GridType::LevelGridView>;
+    using Problem = SymmetricSampleProblem<1, typename GridType::LevelGridView>;
     Problem p(grid.levelGridView(grid.maxLevel()));
 
-    const auto verbosity = Solver::REDUCED;
+    const auto verbosity = Solver::QUIET;
     const bool relativeErrors = false;
 
     using Matrix = typename Problem::Matrix;
     using Vector = typename Problem::Vector;
 
-    Analyser<Vector> analyser(p.u_ex, p.energyNorm, solveTol);
-
     using LoopSolver = ::LoopSolver<Vector>;
     using Step = LinearIterationStep<Matrix, Vector>;
 
-    auto test = [&](Step* step, std::string name) {
+    auto solve = [&](Step* step, double stepTol, size_t maxIterations) {
       Vector u_copy = p.u;
-      Vector rhs_copy = p.rhs;
-
-      step->setProblem(p.A, u_copy, rhs_copy);
-      typename Step::BitVector ignore(u_copy.size(), false);
-      step->ignoreNodes_ = &ignore;
+      step->setProblem(p.A, u_copy, p.rhs);
+      step->ignoreNodes_ = &p.ignore;
 
       LoopSolver solver(step, maxIterations, stepTol,
                         &p.energyNorm, verbosity, relativeErrors);
@@ -90,7 +52,25 @@ struct GSTestSuite
       solver.preprocess();
       solver.solve();
 
-      passed = passed and analyser.analyse(u_copy, name);
+      return u_copy;
+    };
+
+    auto analyse = [&](const Vector& v, std::string testCase, double tolerance) {
+      const auto normDiff = p.energyNorm.diff(p.u_ex, v);
+      const std::string s = Dune::formatString("%*s", -60, testCase.c_str());
+      std::cout << "error for solver \"" << s
+                << "\": " << normDiff << std::endl;
+      if (normDiff > tolerance) {
+        std::cerr << "Error too large for solver \"" << testCase << "\"."
+                  << std::endl;
+        return false;
+      }
+      return true;
+    };
+
+    auto test = [&](Step* step, std::string name) {
+      auto result = solve(step, 1e-12, 2000);
+      passed = passed and analyse(result, name, 1e-7);
     };
 
     {
@@ -98,14 +78,34 @@ struct GSTestSuite
       test(&gsStep, "BlockGS");
     }
 
+    // test regularizable block GS
+    for(auto reg : {false, true}) {
+      for(auto type : {"Direct", "LDLt", "CG"}) {
+        Dune::ParameterTree config;
+        config["type"] = type;
+        config["regularize_diagonal"] = reg ? "1" : "0"; // for SemiDefiniteBlockGS
+        config["regularize_diag"] = reg ? "1e-10" : "0.0"; // for Dune::Solvers::*BlockGS
+        config["maxIter"] = "100"; // for Dune::Solvers::CGBlockGS
+        config["tol"] = "1e-10"; // for Dune::Solvers::CGBlockGS
+
+        {
+          SemiDefiniteBlockGSStep<Matrix, Vector> gsStep(config);
+          test(&gsStep, Dune::formatString("SemiDefiniteBlockGS %s %s",
+                                           type, reg ? "regularized" : ""));
+        }
+      }
+    }
+
+    // test projected block GS
     {
       size_t size = p.u.size();
       using GSStep = ProjectedBlockGSStep<Matrix, Vector>;
-      ProjectedBlockGSStep<Matrix, Vector> gsStep;
       typename GSStep::HasObstacle hasObstacle(size, false);
-      gsStep.hasObstacle_ = &hasObstacle;
-      test(&gsStep, "ProjectedBlockGS w/o obstacle");
-
+      {
+        ProjectedBlockGSStep<Matrix, Vector> gsStep;
+        gsStep.hasObstacle_ = &hasObstacle;
+        test(&gsStep, "ProjectedBlockGS free");
+      }
       hasObstacle.setAll();
       std::vector<typename GSStep::Obstacle> obstacles(size);
       for(size_t i=0; i<size; ++i) {
@@ -115,19 +115,11 @@ struct GSTestSuite
           obstacles[i].upper(j) = j+1;
         }
       }
-      gsStep.obstacles_ = &obstacles;
-      test(&gsStep, "ProjectedBlockGS w obstacle");
-    }
-
-    for(auto reg : {false, true}) {
-      for(auto type : {"direct", "ldlt", "cg"}) {
-        Dune::ParameterTree config;
-        config["type"] = type;
-        config["regularize_diagonal"] = reg ? "1" : "0";
-
-        SemiDefiniteBlockGSStep<Matrix, Vector> gsStep(config);
-        test(&gsStep, Dune::formatString("SemiDefiniteBlockGS %s %s",
-                                         type, reg ? "regularized" : ""));
+      {
+        ProjectedBlockGSStep<Matrix, Vector> gsStep;
+        gsStep.hasObstacle_ = &hasObstacle;
+        gsStep.obstacles_ = &obstacles;
+        test(&gsStep, "ProjectedBlockGS obstacle");
       }
     }