[petsc-users] Change Amat in FormJacobian

[Anton Popov]

Hi Barry & Matt,

thanks for your quick response. These options were exactly what I needed 
and expected:

-pc_mg_galerkin pmat
-pc_use_amat false

I just assumed that it’s a default behavior of the PC object.

So to clarify my case, I don't use nonlinear multigrid. Galerkin is 
expected to deal with Pmat only, and it's enough if Amat implements a 
matrix-vector product for the Krylov accelerator.

Matt, the reason for switching Amat during the iteration is a quite 
common Picard-Newton combination. Jacobian matrix gives accurate updates 
close to the solution, but is rather unstable far form the solution. 
Picard matrix (approximate Jacobian) is quite the opposite – it’s kind 
of stable, but slow. So the idea is to begin the iteration with Picard 
matrix, and switch to the Jacobian later.

If the assembled matrices are used, then the standard SNES interface is 
just perfect. I can decide how to fill the matrices. But I don’t bother 
with  Jacobian assembly and want to use a built-in MFFD approximation 
instead. I did quite a few tests previously and figured out that MFFD is 
practically the same as closed-from matrix-free Jacobian for the later 
stages of the iteration. The Picard matrix still does a good job as a 
preconditioner. But it is important to start the iteration with Picard 
and only change to MFFD later.

Is my workaround with the shell matrix acceptable, or there is a better 
solution?

Thanks,
Anton

On 13.06.21 20:52, Barry Smith wrote:
>
>   Anton,
>
>   -pc_mg_galerkin pmat
>
>   Though it seems simple, there is some subtly in swapping out 
> matrices with SNES.
>
>   When using multigrid with SNES there are at least five distinct uses 
> of the Jacobian operator.
>
>       * Perform matrix-vector product in line search to check Wolf
>         sufficient decrease convergence criteria
>       * Perform the matrix-vector product for the Krylov accelerator
>         of the system
>       * Perform smoothing on the finest level of MG
>       * Perform the matrix-vector product needed on the finest level
>         of MG to compute the residual that will be restricted to the
>         coarser level of multigrid
>       * When using Galerkin products to compute the coarser grid
>         operator performing the Galerkin matrix triple product
>
>
> when one swaps out the mat, which of these do they wish to change? The 
> first two seem to naturally go together as do the last three. In your 
> case I presume you want to swap for the first two, but always use pmat 
> for the last three? To achieve this you will also need -pc_use_amat  false
>
> If you run with -info and -snes_view it will print out some of the 
> information about which operator it is using for each case, but 
> probably not all of them.
>
> Note: if the pmat is actually an accurate computation of the Jacobian 
> then it is likely best not to ever use a matrix-free product. It is 
> only when pmat is approximated in some specific way that using the 
> matrix-free product would be useful to insure the "Newton" method 
> actually computes a Newton step.
>
> Barry
>
>
>
>> On Jun 13, 2021, at 11:21 AM, Matthew Knepley <knepley at gmail.com 
>> <mailto:knepley at gmail.com>> wrote:
>>
>> On Sun, Jun 13, 2021 at 10:55 AM Anton Popov <popov at uni-mainz.de 
>> <mailto:popov at uni-mainz.de>> wrote:
>>
>>     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.
>>
>>
>> Hi Anton,
>>
>> You are correct that there is no specific test, so I can believe that 
>> a bug might be lurking here.
>> I believe the way it is supposed to work is that you set the type of 
>> Galerkin coarsening that you
>> want
>>
>> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCMGSetGalerkin.html
>>
>> So I am thinking you want 'pmat' only, and you would be in charge of 
>> making the coarse MFFD operator
>> and telling PCMG about it. I could imagine that you might want us to 
>> automatically do that, but we would
>> need some way to know that it is MFFD, and with the scheme above we 
>> do not have that.
>>
>> What is the advantage of switching representations during the Newton 
>> iteration?
>>
>>   Thanks,
>>
>>      Matt
>>
>>     Thanks.
>>     Anton
>>
>> -- 
>> 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/20210614/74e2c88c/attachment-0001.html>