[petsc-users] MatSchurComplementGetPmat voes

[Matthew Knepley]

On Fri, Jun 3, 2022 at 9:09 AM Arne Morten Kvarving via petsc-users <
petsc-users at mcs.anl.gov> wrote:

> Hi!
>
> I have a Chorin pressure correction solver with consistent pressure
> update, i.e.
> pressure solve is based on the Schur complement
>
> E = -A10*ainv(A00)*A01
>
> with A10 = divergence, A00 the mass matrix and A01 the gradient.
>
> I have had this implemented with petsc for a long time and it's working
> fine. However, I've done the schur-complement manually, ie using a MatShell.
>
> I now wanted to see if I can implement this using the petsc facilities for
> the schur-complement, but I get a confusing error when I call
> MatSchurComplementGetPmat().
>
> -----
>
> Code snippet:
>
>  MatCreateSchurComplement(m_blocks[0], m_blocks[0], m_blocks[1],
> m_blocks[2], nullptr, &E_operator);
> < ... setup the ksp for A00 >
>  MatSchurComplementSetAinvType(E_operator, MAT_SCHUR_COMPLEMENT_AINV_DIAG);
> MatView(E_operator);
> MatSchurComplementGetPmat(E_operator, MAT_INITIAL_MATRIX, &E_pc);
>
> -----
>
> This yields the output (I cut out the matrix elements):
> Mat Object: 1 MPI processes
>   type: schurcomplement
>   Schur complement A11 - A10 inv(A00) A01
>   A11 = 0
>   A10
>     Mat Object: 1 MPI processes
>       type: seqaij
>  KSP of A00
>     KSP Object: 1 MPI processes
>       type: preonly
>       maximum iterations=1000, initial guess is zero
>       tolerances:  relative=1e-06, absolute=1e-20, divergence=1e+06
>       left preconditioning
>       using DEFAULT norm type for convergence test
>     PC Object: 1 MPI processes
>       type: lu
>         out-of-place factorization
>         tolerance for zero pivot 2.22045e-14
>         matrix ordering: nd
>         factor fill ratio given 5., needed 1.02768
>           Factored matrix follows:
>             Mat Object: 1 MPI processes
>               type: seqaij
>               rows=72, cols=72
>               package used to perform factorization: petsc
>               total: nonzeros=4752, allocated nonzeros=4752
>                 using I-node routines: found 22 nodes, limit used is 5
>       linear system matrix = precond matrix:
>       Mat Object: 1 MPI processes
>         type: seqaij
>         rows=72, cols=72
>         total: nonzeros=4624, allocated nonzeros=5184
>         total number of mallocs used during MatSetValues calls=0
>           using I-node routines: found 24 nodes, limit used is 5
>   A01
>     Mat Object: 1 MPI processes
>       type: seqaij
>
> [0]PETSC ERROR: --------------------- Error Message
> --------------------------------------------------------------
> [0]PETSC ERROR: Invalid argument
> [0]PETSC ERROR: Wrong type of object: Parameter # 1
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.17.2, Jun 02, 2022
> [0]PETSC ERROR: ../d/bin/Stokes on a linux-gnu-cxx-opt named akvalung by
> akva Fri Jun  3 14:48:06 2022
> [0]PETSC ERROR: Configure options --with-mpi=0 --with-lapack-lib=-llapack
> --with-64-bit-indices=0 --with-shared-libraries=0 --COPTFLAGS=-O3
> --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-blas-lib=-lblas --CFLAGS=-fPIC
> --CXXFLAGS=-fPIC --FFLAGS=-fPIC
> [0]PETSC ERROR: #1 MatDestroy() at
> /home/akva/kode/petsc/petsc-3.17.2/src/mat/interface/matrix.c:1235
> [0]PETSC ERROR: #2 MatCreateSchurComplementPmat() at
> /home/akva/kode/petsc/petsc-3.17.2/src/ksp/ksp/utils/schurm/schurm.c:763
> [0]PETSC ERROR: #3 MatSchurComplementGetPmat_Basic() at
> /home/akva/kode/petsc/petsc-3.17.2/src/ksp/ksp/utils/schurm/schurm.c:785
> [0]PETSC ERROR: #4 MatSchurComplementGetPmat() at
> /home/akva/kode/petsc/petsc-3.17.2/src/ksp/ksp/utils/schurm/schurm.c:835
>
> where the errors come from the call call to obtain the preconditioner
> matrix.
> I don't see what I've done wrong, as far as I can see it's all following
> https://petsc.org/release/docs/manualpages/KSP/MatCreateSchurComplement.html#MatCreateSchurComplement
> MatCreateSchurComplement - Argonne National Laboratory
> <https://petsc.org/release/docs/manualpages/KSP/MatCreateSchurComplement.html#MatCreateSchurComplement>
> Notes The Schur complement is NOT explicitly formed! Rather, this function
> returns a virtual Schur complement that can compute the matrix-vector
> product by using formula S = A11 - A10 A^{-1} A01 for Schur complement S
> and a KSP solver to approximate the action of A^{-1}.. All four matrices
> must have the same MPI communicator.
> petsc.org
> **
>
> and
>
> https://petsc.org/release/docs/manualpages/KSP/MatSchurComplementGetPmat.html#MatSchurComplementGetPmat
>
> Looking into the code it seems to try to call MatDestroy() for the Sp
> matrix but as Sp has not been set up it fails (schurm.c:763)
> Removing that call as a test, it seems to succeed and I get the same
> solution as I do
> with my manual code.
>
> I'm sure I have done something stupid but I cannot see what, so any
> pointers would be appreciated.
>

This is not your fault. If the flag is MAT_INITIAL_MATRIX, we are expecting
the pointer to be initialized to NULL, but we never state this.
I think if you do this, the code will start working. I will fix GetPmat()
so that it does this automatically.

  Thanks,

     Matt


> cheers
>
> arnem
>
>

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

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220603/3c1f5141/attachment.html>