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

[Eike Mueller]

Dear dune-pdelab,

just to say that, with a lot of help from the dune developers in  
Heidelberg, this has now been solved. Actually it was enough to  
reorder the (local) IndexSet in the gridview class, which is then used  
to construct the grid function space. This can be done quite easily by  
deriving a class from the appropriate gridview which contains a new,  
reordered index set (basically realised by a lookup table, which is  
constructed by traversing the grid in the constructor).
As I always have the same number of cells in the vertical direction,  
the vector/matrix stored by the ISTL backends then has the desired  
block structure, i.e. it can be split into blocks of the length n_z,  
and I don't even need a nested date structure such as  
BCRSMatrix<BTDMatrix<...>>. I then use the SeqOverlappingSchwarz  
preconditioner (with ILU0, which is exact on the blocks), which needs  
to know the dofs that are blocked together, but that is easy as the  
layout of the vector is known because of the regular ordering. It  
works both for YaspGrid and ALUGrid.

For my problem with strong coupling in the vertical direction the line  
smoother as a preconditioner alone is very efficient.

Thank you very much for all your help.

Eike

On 2 Dec 2011, at 09:43, Steffen Müthing wrote:

> 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
>
>
>
>
>
>
> _______________________________________________
> dune-pdelab mailing list
> dune-pdelab at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-pdelab
>