[dune-pdelab] Implementing tri-diagonal matrices in DUNE/PDELab

[Eike Mueller]

Dear dune-pdelab,

I already posted this on the main dune mailing list, but as it is  
really a pdelab-question I'm reposting it here. Apologies for any  
confusion!

I'm trying to exploit the strong vertical coupling of the problem I'm  
solving (essentially a 3d Poisson equation, discretised by a 7 point  
FV stencil) by using a line- instead of a point smoother and have  
implemented this using the BTDMatrix class (i.e. my matrix is a block  
matrix, where each block is tridiagonal and the vertical degrees of  
freedom for each horizontal position (x,y) are grouped together in a  
block of the solution vector).

I tried to build this into the PDELab framework by writing a new  
version of ISTLVectorBackend, based on

Dune::BlockVector<Dune::BlockVector<Dune::BlockFieldVector<E,1>>>

instead of

Dune::BlockVector<BlockFieldVector<E,BLOCKSIZE>>

and a modified version of ISTLBCRSMatrixBackend, based on

Dune::BCRSMatrix<Dune::BTDMatrix<Dune::FieldMatrix<E,1,1>>>

instead of

Dune::BCRSMatrix<Dune::FieldMatrix<E,BLOCKROWSIZE,BLOCKCOLSIZE>>

and this all works fine, i.e. I really just need to replace the  
backends in my code and everything else stays the same:

[...]
// (2) Construct finite element space
    typedef Dune::PDELab::P0LocalFiniteElementMap<Coord, double, dim>  
FEM;
    FEM fem(Dune::GeometryType(Dune::GeometryType::cube,dim));
    typedef Dune::PDELab::OverlappingConformingDirichletConstraints CON;

    // (3) Set up grid function space
    typedef ISTLBlockVectorBackend<32> VBE;
//    typedef Dune::PDELab::ISTLVectorBackend<1> VBE;
    typedef Dune::PDELab::GridFunctionSpace<GV, FEM, CON, VBE> GFS;
    GFS gfs(gv, fem);

    // (4) Set up grid operator space
    LocalOperator lop(epsilon2);
    typedef ISTLBlockTridiagMatrixBackend<32> MBE;
//    typedef Dune::PDELab::ISTLBCRSMatrixBackend<1,1> MBE;
    typedef Dune::PDELab::EmptyTransformation CC;
    typedef Dune::PDELab::GridOperatorSpace<GFS, GFS, LocalOperator,  
CC, CC, MBE> GOS;
    GOS gos(gfs,gfs,lop);
[...]

If I solve a Poisson problem on a unit cube with a strong anisotropy  
in the x-direction and use SSOR as a preconditioner (i.e. an x-line  
smoother), the tridiagonal solve is very efficient and the solver  
converges much faster than just by using the default Vector/Matrix  
backends (i.e. a point smoother). The only issues I still have are:

* I have to specify the size of the x-direction, i.e. the block size  
as a template parameter to my backends (I'd rather specify it at  
runtime) and
* I don't seem to have any control over how the degrees of freedom are  
mapped from the grid onto the vector container, for example in my case  
all points on an x-line are stored in one block of the block vector,  
but I actually want the vertical degrees at a given (x,y) to be  
grouped together. I guess I could modify the access() methods in the  
backends, but the backends do not know anything about the grid  
function space.

Maybe rewriting the backends is also not the right way forward and the  
grid function space class needs to be adapted as well?

Thanks,

Eike