Improve the procedure of identifying blocks in dofmap

Issue #193 resolved
Johan Hake created an issue

The method:

DofMapBuilder::compute_blocksize

can only identify dofmap blocks of VectorFunctionSpace or from a MixedFunctionSpace where all subspaces are the same. Blocks are not identified for MixedFunctionSpace using a mixture of sub spaces like P1-P2.

With this in place we could:

  1. use MatNest for efficient block preconditioning
  2. block dofs within each subspace for efficiency
  3. more easily attach null spaces, etc (with are required for some preconditioners).
  4. make FunctionAssigner work for more than VectorFunctionSpaces

We should also remove all global (Real) dofs from the dof reordering graph algorithm and position these at a fixed place in the dofmap, for example the first dofs. Global dofs are connected to all other dofs, making the procedure of building the graph very time consuming and it does not add any information to the graph.

It would be nice to have available different strategies for how the dofs are ordered. For the FunctionAssigner to work we need all subspaces of the same family and degree to be ordered in the same way. For this to happen a subspace based dof block identification need to be possible.

Comments (17)

  1. Jan Blechta

    I believe that global DOFs should be placed in the end. This should help to sparse direct solvers as they typically need denser rows in the end of a matrix. This produces smaller fill-in during factorization.

  2. Jan Blechta

    @garth-wells Sure, but why to make it harder for these algorithms. Moreover, they are not bullet-proof.

  3. Johan Hake reporter

    @blechta putting it in the front was just a suggestion based on having a consistent rule.

  4. Nico Schlömer

    Johan correctly points out that it appears that many of the issues that we're having at the moment with block ("multiphysics") problems could be handled if we more closely adapt to PETSc's MatNest structure. For problems like Stokes, we'd probably have to climb down to the assembler to required it not to produce one large matrix, but four blocks subsequently stitched together by MatCreateNest. Block DOFs, subvector assignment, sub-functionspace access and everything else that is marked as a duplicate of this bug would be close to trivially solved. Something that I don't see clearly here is how (I)LU and AMG methods act on (nested) block matrices. (This is disregarding the fact that one wouldn't want to do that in the first place and rather exploit the block structure first.)

  5. Nico Schlömer

    @garth-wells The approach here goes a little further than #119 in that it suggests to adapt the assembler to produce matrix/vector blocks for MatCreateNest.

  6. Prof Garth Wells

    @nschloe The key ingredient for MatNest is to switch to local dof indexing. It is then possible to seamlessly create PETSc monolithic or nested matrices without changing the assembler.

  7. Jan Blechta

    @johanhake Is there another bug or will it be covered by #119?

    from dolfin import *
    
    mesh = UnitSquareMesh(20,20,"crossed")
    
    P = FunctionSpace(mesh, "Lagrange", 1)
    W = FunctionSpace(mesh, "Lagrange", 1)
    R = VectorFunctionSpace(mesh, "Lagrange", 2)
    M = VectorFunctionSpace(mesh, "Lagrange", 2)
    V = MixedFunctionSpace([P, W, R, M])
    
    assigner = FunctionAssigner(V.sub(0),P)
    

    returns in parallel

    *** Error:   Unable to create function assigner.
    *** Reason:  The receiving and assigning spaces do not have the same number of dofs per space.
    *** Where:   This error was encountered inside FunctionAssigner.cpp.
    
  8. Log in to comment