[petsc-users] Change Amat in FormJacobian

[Anton Popov]

Hi,

I want a simple(?) thing. I want to decide and be able to assign the 
Jacobian matrix (Amat) on the fly within the FormJacobian function (i.e. 
during Newton iteration) to one of the following options:

1) built-in MFFD approximation
2) assembled preconditioner matrix (Pmat)

I have not found this scenario demonstrated in any example, therefore 
I’m asking how to do that.

Currently I do the following:

1) setup Amat as a shell matrix with a MATOP_MULT operation that simply 
retrieves a matrix object form its context and calls MatMult on it.

2) if I need MFFD, I put a matrix generated with MatCreateSNESMF in the 
Amat context (of course I also call MatMFFDComputeJacobian before that).

3) if I need Pmat, I simply put Pmat in the Amat context.

4) call MatAssemblyBegin/End on Amat

So far so good.

However, shell Amat and assembled Pmat generate a problem if Galerkin 
multigrid is requested as a preconditioner (I just test on 1 CPU):

[0]PETSC ERROR: MatPtAP requires A, shell, to be compatible with P, 
seqaij (Misses composed function MatPtAP_shell_seqaij_C)
[0]PETSC ERROR: #1 MatPtAP()
[0]PETSC ERROR: #2 MatGalerkin()
[0]PETSC ERROR: #3 PCSetUp_MG()
[0]PETSC ERROR: #4 PCSetUp()
[0]PETSC ERROR: #5 KSPSetUp()
[0]PETSC ERROR: #6 KSPSolve()
[0]PETSC ERROR: #7 SNESSolve_NEWTONLS()
[0]PETSC ERROR: #8 SNESSolve()

It seems like PETSc tries to apply Galerkin coarsening to the shell Amat 
matrix instead of the assembled Pmat. Why?

Pmat is internally generated by SNES using a DMDA attached to SNES, so 
it has everything to define Galerkin coarsening. And it actually works 
if I just get rid of the shell Amat.

I can get around this problem by wrapping a PC object in a shell matrix 
and pass it as Pmat. This is a rather ugly approach, so I wonder how to 
resolve the situation in a better way, if it is possible.

Thanks.
Anton