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>

// 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/multigridstep.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.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>
struct MultigridTestSuite
{
    template <class GridType>
    bool check(const GridType& grid)
    {
        double tol = 1e-12;
        int maxIterations = 200;

        bool passed = true;

        typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
        typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
        typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
        typedef typename Dune::BlockVector<VectorBlock> Vector;
        typedef typename Dune::BitSetVector<blocksize> BitVector;

        typedef typename GridType::LevelGridView GridView;


        const GridView gridView = grid.levelGridView(grid.maxLevel());

        Matrix A;
        constructPQ1Pattern(gridView, A);
        A=0.0;
        assemblePQ1Stiffness(gridView, A);

        Vector u(A.N()), u_ex(A.N());
        u = 0;
        for (size_t i=0; i<u.size(); ++i)
            for (size_t j=0; j<blocksize; ++j)
            {
                u_ex[i][j] = (1.0*rand())/RAND_MAX;
//                A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]); // in case of no dirichlet values you can make the matrix pd hereby
            }

        Vector rhs(A.N());
        A.mv(u_ex, rhs);

        BitVector ignore(A.N());
        ignore.unsetAll();
        markBoundaryDOFs(gridView, ignore);

        // set dirichlet values
        for (size_t i=0; i<u.size(); ++i)
            for (size_t j=0; j<blocksize; ++j)
                if (ignore[i][j])
                    u[i][j] = u_ex[i][j];

        typedef EnergyNorm<Matrix, Vector> Norm;
        typedef ::LoopSolver<Vector> Solver;
        typedef MultigridStep<Matrix, Vector, BitVector> MGStep;
        typedef MultigridTransfer<Vector, BitVector, Matrix> Transfer;
        typedef CompressedMultigridTransfer<Vector, BitVector, Matrix> TransferImplementation;
        typedef TruncatedBlockGSStep<Matrix, Vector> Smoother;

        // we need a vector of pointers to the transfer operator base class
        std::vector<Transfer*> transfer(grid.maxLevel());
        for (int 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;
        }

        // we want to control errors with respect to the energy norm induced by the matrix
        Norm energynorm(A);

        // set up smoothers and basesolver
        Smoother smoother, basesolver_step;

        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(A,u,rhs);
        mgStep.setSmoother(&smoother);
        mgStep.setMGType(1,3,3);
        mgStep.ignoreNodes_ = &ignore;
        mgStep.basesolver_  = &basesolver;

        // create loop solver
        Solver solver(&mgStep, maxIterations, tol, &energynorm, Solver::FULL);

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

        u = mgStep.getSol();

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


        // now solve with the same operator but different rhs (and dirichlet values)

        // create new exact solution
        for (size_t i=0; i<u.size(); ++i)
            for (size_t j=0; j<blocksize; ++j)
                u_ex[i][j] = (1.0*rand())/RAND_MAX;

        // ... and its corresponding rhs
        A.mv(u_ex, rhs);

        // set initial iterate (remember that the iterationStep member points to u)
        u = 0.0;
        // set dirichlet values
        for (size_t i=0; i<u.size(); ++i)
            for (size_t j=0; j<blocksize; ++j)
                if (ignore[i][j])
                    u[i][j] = u_ex[i][j];

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

        solver.solve();

        u = mgStep.getSol();

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

        return passed;
    }
};



int main(int argc, char** argv) try
{
    bool passed(true);

    MultigridTestSuite<1> testsuite1;
//    MultigridTestSuite<2> testsuite2;
//    MultigridTestSuite<8> testsuite8;
    passed = checkWithStandardGrids(testsuite1);
//    passed = checkWithStandardGrids(testsuite2);
//    passed = checkWithStandardGrids(testsuite8);

    return passed ? 0 : 1;
}
catch (Dune::Exception e) {

    std::cout << e << std::endl;

}