Add test for TNNMGStep

This is just a modified version of the NonlinearGSTest using the newer TNNMGStep.
Up to now, commits that lead TNNMGStep based implementations to fail would not be covered by the CI.

dune/tnnmg/test/CMakeLists.txt

@@ -6,3 +6,4 @@ dune_add_test(SOURCES nonlineargstest)
 dune_add_test(SOURCES nonlineargsperformancetest NAME nonlineargsperformancetest_corrections)
 dune_add_test(SOURCES nonlineargsperformancetest NAME nonlineargsperformancetest_iterates COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
 dune_add_test(SOURCES simplexsolvertest)
+dune_add_test(SOURCES tnnmgsteptest)

dune/tnnmg/test/tnnmgsteptest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include <config.h>
+
+#include <iostream>
+#include <sstream>
+
+// dune-common includes
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/stringutility.hh>
+#include <dune/common/test/testsuite.hh>
+
+// dune-istl includes
+#include <dune/istl/bcrsmatrix.hh>
+
+// dune-grid includes
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+// dune-solver includes
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/common/defaultbitvector.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+
+// dune-tnnmg includes
+#include <dune/tnnmg/functionals/quadraticfunctional.hh>
+#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
+#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
+
+#include <dune/tnnmg/projections/obstacledefectprojection.hh>
+
+#include <dune/solvers/test/common.hh>
+
+template <class DomainType, class RangeType, class F>
+class FunctionFromLambda :
+    public Dune::VirtualFunction<DomainType, RangeType>
+{
+    public:
+        FunctionFromLambda(F f):
+            f_(f)
+    {}
+
+        void evaluate(const DomainType& x, RangeType& y) const
+        {
+            y = f_(x);
+        }
+
+    private:
+        F f_;
+};
+
+template<class Domain, class Range, class F>
+auto makeFunction(F&& f)
+{
+    return FunctionFromLambda<Domain, Range, std::decay_t<F>>(std::forward<F>(f));
+}
+
+
+// TODO include those last two parameters
+template<class MatrixType, class VectorType, class BitVector, class TransferOperators>
+void solveProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs, const BitVector& ignore, const TransferOperators& transfer, int maxIterations=100, double tolerance=1.0e-8)
+{
+    using Vector = VectorType;
+    using Matrix = MatrixType;
+
+    auto lower = Vector(rhs.size());
+    auto upper = Vector(rhs.size());
+
+    lower = 0;
+    upper = 1;
+
+    using Functional = Dune::TNNMG::BoxConstrainedQuadraticFunctional<Matrix&, Vector&, Vector&, Vector&, double>;
+    auto J = Functional(mat, rhs, lower, upper);
+
+    auto localSolver = gaussSeidelLocalSolver(Dune::TNNMG::ScalarObstacleSolver());
+
+    using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, decltype(localSolver), BitVector>;
+    auto nonlinearSmoother = std::make_shared<NonlinearSmoother>(J, x, localSolver);
+
+    auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+
+    auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+    linearMultigridStep->setMGType(1, 3, 3);
+    linearMultigridStep->setSmoother(&smoother);
+    linearMultigridStep->setTransferOperators(transfer);
+
+    using Linearization = Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<Functional, BitVector>;
+    using DefectProjection = Dune::TNNMG::ObstacleDefectProjection;
+    using LineSearchSolver = Dune::TNNMG::ScalarObstacleSolver;
+    using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+
+    using Solver = LoopSolver<Vector>;
+    using Norm =  EnergyNorm<Matrix, Vector>;
+
+
+    using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+    int mu=1; // #multigrid steps in Newton step
+    auto step = Step(J, x, nonlinearSmoother, linearMultigridStep, mu, DefectProjection(), LineSearchSolver());
+
+    auto norm = Norm(mat);
+    auto solver = Solver(step, 1e9, 0, norm, Solver::FULL);
+
+    step.setIgnore(ignore);
+    step.setPreSmoothingSteps(3);
+
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12.5e", J(x));
+            },
+            "   energy      ");
+
+    double initialEnergy = J(x);
+    solver.addCriterion(
+            [&](){
+            static double oldEnergy=initialEnergy;
+            double currentEnergy = J(x);
+            double decrease = currentEnergy - oldEnergy;
+            oldEnergy = currentEnergy;
+            return Dune::formatString("   % 12.5e", decrease);
+            },
+            "   decrease    ");
+
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12.5e", step.lastDampingFactor());
+            },
+            "   damping     ");
+
+
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12d", step.linearization().truncated().count());
+            },
+            "   truncated   ");
+
+
+    std::vector<double> correctionNorms;
+    solver.addCriterion(Dune::Solvers::correctionNormCriterion(step, norm, 1e-10, correctionNorms));
+
+    solver.preprocess();
+    solver.solve();
+    std::cout << correctionNorms.size() << std::endl;
+}
+
+
+
+
+
+template <class GridType>
+Dune::TestSuite checkWithGrid(const GridType& grid)
+{
+    Dune::TestSuite suite;
+
+    static const int dim = GridType::dimension;
+
+    using MatrixBlock = typename Dune::FieldMatrix<double, 1, 1>;
+    using VectorBlock = typename Dune::FieldVector<double, 1>;
+    using Matrix = typename Dune::BCRSMatrix<MatrixBlock>;
+    using Vector = typename Dune::BlockVector<VectorBlock>;
+    using BitVector = typename Dune::Solvers::DefaultBitVector_t<Vector>;
+
+    using GridView = typename GridType::LeafGridView;
+    using DomainType = typename Dune::FieldVector<double, dim>;
+    using RangeType = typename Dune::FieldVector<double, 1>;
+
+
+    const GridView gridView = grid.leafGridView();
+
+    Matrix A;
+    constructPQ1Pattern(gridView, A);
+    A = 0.0;
+    assemblePQ1Stiffness(gridView, A);
+
+    Matrix M;
+    constructPQ1Pattern(gridView, M);
+    M = 0.0;
+    assemblePQ1Mass(gridView, M);
+    A += M;
+
+    Vector rhs(A.N());
+    rhs = 0;
+    auto f = makeFunction<DomainType, RangeType>([](const DomainType& x) { return (x.two_norm() < .7) ? 200 : 0;});
+    assemblePQ1RHS(gridView, rhs, f);
+
+    Vector u(A.N());
+    u = 0;
+
+    BitVector ignore(A.N());
+    ignore.unsetAll();
+    markBoundaryDOFs(gridView, ignore);
+
+
+    using TransferOperator = CompressedMultigridTransfer<Vector>;
+    using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+    TransferOperators transfer(grid.maxLevel());
+    for (size_t i = 0; i < transfer.size(); ++i)
+    {
+        // create transfer operator from level i to i+1
+        transfer[i] = std::make_shared<TransferOperator>();
+        transfer[i]->setup(grid, i, i+1);
+    }
+
+    solveProblem(A, u, rhs, ignore, transfer);
+
+    // TODO: let suite check something
+    return suite;
+}
+
+
+template <int dim>
+Dune::TestSuite checkWithYaspGrid(int refine)
+{
+    using GridType = Dune::YaspGrid<dim>;
+
+    typename Dune::FieldVector<typename GridType::ctype,dim> L(1.0);
+    typename std::array<int,dim> s;
+    std::fill(s.begin(), s.end(), 1);
+
+    GridType grid(L, s);
+
+    for (int i = 0; i < refine; ++i)
+        grid.globalRefine(1);
+
+    std::cout << "Testing with YaspGrid<" << dim << ">" << std::endl;
+    std::cout << "Number of levels: " << (grid.maxLevel()+1) << std::endl;
+    std::cout << "Number of leaf nodes: " << grid.leafGridView().size(dim) << std::endl;
+
+    auto suite = checkWithGrid(grid);
+
+    return suite;
+}
+
+
+
+
+int main(int argc, char** argv) try
+{
+    Dune::MPIHelper::instance(argc, argv);
+    Dune::TestSuite suite;
+
+
+    int refine1d = 1;
+    int refine2d = 6;
+    int refine3d = 1;
+
+    if (argc>1)
+        std::istringstream(argv[1]) >> refine1d;
+    if (argc>2)
+        std::istringstream(argv[2]) >> refine2d;
+    if (argc>3)
+        std::istringstream(argv[3]) >> refine3d;
+
+    suite.subTest(checkWithYaspGrid<1>(refine1d));
+    suite.subTest(checkWithYaspGrid<2>(refine2d));
+
+    return suite.exit();
+}
+catch (Dune::Exception e) {
+
+    std::cout << e << std::endl;
+
+}
+