multigridtest.cc

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:

#include <config.h>

#include <cmath>
#include <iostream>
#include <sstream>

// dune-common includes
#include <dune/common/bitsetvector.hh>
#include <dune/common/parallel/mpihelper.hh>

// dune-istl includes
#include <dune/istl/bcrsmatrix.hh>

// dune-grid includes

// dune-solver includes
#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>


#include "common.hh"

/**  \brief test for the MultigridStep class
  *
  *  This test tests if the Dirichlet problem for a laplace operator is solved correctly for a "random" rhs by a LoopSolver employing the MultigridStep.
  *  It furthermore tests the functionality to solve correctly for a different rhs with otherwise unaltered data without calling preprocess().
  *  Setting and using different smoothers in different levels is NOT a tested feature.
  */
template <size_t blocksize, bool trivialDirichletOnly = true>
struct MultigridTestSuite
{
    template <class GridType>
    bool check(const GridType& grid)
    {
        double tol = 1e-12;
        int maxIterations = 200;

        bool passed = true;

        using Problem =
            SymmetricSampleProblem<blocksize, typename GridType::LevelGridView, trivialDirichletOnly>;
        Problem p(grid.levelGridView(grid.maxLevel()));

        typedef typename Problem::Vector Vector;
        typedef typename Problem::Matrix Matrix;
        typedef typename Problem::BitVector BitVector;
        typedef ::LoopSolver<Vector> Solver;
        typedef MultigridStep<Matrix, Vector, BitVector> MGStep;
        typedef MultigridTransfer<Vector, BitVector, Matrix> Transfer;
        typedef CompressedMultigridTransfer<Vector, BitVector, Matrix> TransferImplementation;

        // we need a vector of pointers to the transfer operator base class
        std::vector<Transfer*> transfer(grid.maxLevel());
        for (size_t i = 0; i < transfer.size(); ++i)
        {
            // create transfer operator from level i to i+1 (note that this will only work for either uniformly refined grids or adaptive grids with RefinementType=COPY)
            TransferImplementation* t = new TransferImplementation;
            t->setup(grid, i, i+1);
            transfer[i] = t;
        }

        // set up smoothers and basesolver
        auto smoother =
            Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(
                Dune::Solvers::BlockGS::LocalSolvers::gs());
        auto basesolver_step =
            Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(
                Dune::Solvers::BlockGS::LocalSolvers::gs());

        typename Problem::Norm basenorm(basesolver_step);

        Solver basesolver(&basesolver_step, 1, 1e-10, &basenorm, Solver::QUIET);

        // create iteration step and set all kinds of data
        MGStep mgStep;
        mgStep.setTransferOperators(transfer);
        mgStep.setProblem(p.A,p.u,p.rhs);
        mgStep.setSmoother(&smoother);
        mgStep.setMGType(1,3,3);
        mgStep.ignoreNodes_ = &p.ignore;
        mgStep.basesolver_  = &basesolver;

        // create loop solver
        Solver solver(&mgStep, maxIterations, tol, &p.energyNorm, Solver::FULL);

        // solve problem
        solver.preprocess();
        solver.solve();

        p.u = mgStep.getSol();

        if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
        {
            std::cout << "The MultigridStep did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
            std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
            passed = false;
        }

        p.randomize();

        // reset rhs in mgStep without calling preprocess subsequently
        mgStep.setRhs(p.rhs);

        solver.solve();

        p.u = mgStep.getSol();

        if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
        {
            std::cout << "The MultigridStep did not produce a satisfactory result after calling setRhs(). ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
            std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
            passed = false;
        }

        return passed;
    }
};



int main(int argc, char** argv)
{
    Dune::MPIHelper::instance(argc, argv);
    bool passed(true);

    MultigridTestSuite<1> testsuite1;
    MultigridTestSuite<2> testsuite2;
    MultigridTestSuite<3, false> testsuite3f;
    MultigridTestSuite<8> testsuite8;
    passed = checkWithStandardGrids(testsuite1);
    passed = checkWithStandardGrids(testsuite2);
    passed = checkWithStandardGrids(testsuite3f);
    passed = checkWithStandardGrids(testsuite8);

    return passed ? 0 : 1;
}