[petsc-dev] Preconditioner for Schur complement (pcfieldsplit)

[Thomas Witkowski]

Matthew Knepley wrote:
> On Tue, Apr 5, 2011 at 8:38 AM, Thomas Witkowski 
> <thomas.witkowski at tu-dresden.de 
> <mailto:thomas.witkowski at tu-dresden.de>> wrote:
>
>     Matthew Knepley wrote:
>
>         On Tue, Apr 5, 2011 at 4:28 AM, Thomas Witkowski
>         <thomas.witkowski at tu-dresden.de
>         <mailto:thomas.witkowski at tu-dresden.de>
>         <mailto:thomas.witkowski at tu-dresden.de
>         <mailto:thomas.witkowski at tu-dresden.de>>> wrote:
>
>            Jed Brown wrote:
>
>                On Tue, Mar 29, 2011 at 15:29, Thomas Witkowski
>                <thomas.witkowski at tu-dresden.de
>         <mailto:thomas.witkowski at tu-dresden.de>
>                <mailto:thomas.witkowski at tu-dresden.de
>         <mailto:thomas.witkowski at tu-dresden.de>>
>                <mailto:thomas.witkowski at tu-dresden.de
>         <mailto:thomas.witkowski at tu-dresden.de>
>                <mailto:thomas.witkowski at tu-dresden.de
>         <mailto:thomas.witkowski at tu-dresden.de>>>> wrote:
>
>                   Is there any example which makes use of it?
>
>
>                src/snes/examples/tutorials/ex28.c is a (contrived)
>                multi-physics prototype that uses MatNest tranparently
>         if you
>                run with -pack_dm_mat_type nest (currently, you should add
>                -pc_type fieldsplit -pc_fieldsplit_type additive because
>                off-diagonal blocks are not automatically preallocated).
>
>                There are also two tests that use the low-level API:
>
>                $ grep -l MatNest src/**/tests/*.c
>                src/ksp/ksp/examples/tests/ex22.c
>                src/snes/examples/tests/ex17.c
>
>            So, I reimplemented my code to assemble the matrices to a
>         MatNest
>            matrix. This seems to work fine. But how to define the
>            FieldSplits? I did it in the following way:
>
>              KSP solver;
>              PC pc;
>              KSPCreate(PETSC_COMM_WORLD, &solver);
>              KSPGetPC(solver, &pc);
>              const PetscInt interiorField[1] = {0};
>              const PetscInt boundaryField[1] = {1};
>              PCFieldSplitSetFields(pc, "interior", 1, interiorField);
>              PCFieldSplitSetFields(pc, "boundary", 1, boundaryField);
>
>            When I run my code with the options "-ksp_type preonly -pc_type
>            fieldsplit -pc_fieldsplit_type schur", I get the following
>         errors
>            in KSPSolve:
>
>            [0]PETSC ERROR: --------------------- Error Message
>            ------------------------------------
>            [0]PETSC ERROR: Petsc has generated inconsistent data!
>            [0]PETSC ERROR: Unhandled case, must have at least two fields!
>
>            Is the creation of the split wrong? If not, what is a good
>         way to
>            figure out the problem?
>
>
>         This is almost certainly a problem with the setup order. In
>         this code, the PC type is not yet set
>         so the FieldSplit calls are ignored. When it finally is set,
>         there are no fields defined. I would call
>         KSPSetFromOptions() before retrieving the PC.
>
>     This generates exactly the same error message.
>
>
> Then there is no choice but to use the debugger to see what is going on.
I did so and the problem is as you have described before. But, when I 
understood the implementation of fieldsplit right, than it should not be 
necessary to set the splits using PCFieldSplitSetFields. This could be 
done in PCFieldSplitSetDefaults, when I make use of the parameter 
"-pc_fieldsplit_defaults". The problem is that MatGetBlockSize on my 
MatNest (2x2 blocks) returns 1. Therefore PCFieldSplitSetDefaults 
creates only one split and exits with the error due to

  if (jac->nsplits < 2) 
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unhandled case, must have at 
least two fields");

I have not a very good overview of PETSc's internal matters, but 
wouldn't be it correct that MatGetBlockSize returns the number of nested 
matrices of a MatNest?

Thomas
>
>    Matt
>  
>
>
>     Thomas
>
>
>           Matt
>          
>
>            Thomas
>
>
>
>
>         -- 
>         What most experimenters take for granted before they begin
>         their experiments is infinitely more interesting than any
>         results to which their experiments lead.
>         -- Norbert Wiener
>
>
>
>
>
> -- 
> What most experimenters take for granted before they begin their 
> experiments is infinitely more interesting than any results to which 
> their experiments lead.
> -- Norbert Wiener