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

[Steffen Müthing]

Hello Eike,

> 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:

that's pretty neat! The problem is that by doing this, you have ventured into an area of PDELab that most users
(and developers) rarely, if ever, touch... the LA backend interface is not very fleshed out yet, mostly because we
lacked good examples of more advanced setups to figure out a good interface design.

I'm afraid there are no obvious, "clean" solutions to your questions, but here are some thoughts:

> 
> 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 take it that's for the mapping "global flat index" -> "blocked matrix index" in the backend. The problem here is that the backend only exposes
some static methods and cannot carry any dynamic state (like the block size). Until that changes, I'm afraid the only thing you might be able to
do is to store the block size on the vector / matrix wrappers (those exposed by the backends) and then read them back in from those objects in
the access() methods of the backend. Of course, that would require a non-standard constructor that gets passed the block size in addition to
the GFS / GOS. As I said, it's an ugly hack... ;-)

> * 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.

Well, what you want is a custom reordering of the IndexSet which belongs to the GridView of the GFS. PDELab just arranges the DOFs in the
same way the IndexSet does (and sorts different geometry types according to their natural ordering, given by operator<() ). The precise layout
of the IndexSet actually depends on the grid you used, UG would for example yield a different pattern than what you observe (I take it you
currently use YaspGrid, which iterates (from inner to outer loop) over x,y,z, giving you the blocking in x-direction). If you want to change that,
the easiest way to do so will be to determine the correct re-ordering of the IndexSet and permute the entity indices wherever they are used. If I
remember correctly, you probably only need to touch the implementation of the following methods of GridFunctionSpace in gridfunctionspace.hh:

- globalIndices()
- entityOffset()
- dataHandleGlobalIndices()
- update()

There are three different implementations of GridFunctionSpace in that file, you only need to modify the first one in your case.

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

As I laid out above, I'm afraid you need both...

Steffen

> 
> Thanks,
> 
> Eike
> 
> _______________________________________________
> dune-pdelab mailing list
> dune-pdelab at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-pdelab

Steffen Müthing
Universität Stuttgart
Institut für Parallele und Verteilte Systeme
Universitätsstr. 38
70569 Stuttgart
Tel: +49 711 685 88429
Fax: +49 711 685 88340
Email: steffen.muething at ipvs.uni-stuttgart.de