[petsc-users] MatSchurComplementGetPmat voes

[Arne Morten Kvarving]

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.

cheers

arnem

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220603/405da78e/attachment-0001.html>