[petsc-users] Change Amat in FormJacobian
Anton Popov
popov at uni-mainz.de
Mon Jun 14 10:27:42 CDT 2021
On 14.06.21 15:04, Dave May wrote:
>
> Hi Anton,
Hi Dave,
>
> On Mon, 14 Jun 2021 at 14:47, Anton Popov <popov at uni-mainz.de
> <mailto:popov at uni-mainz.de>> wrote:
>
> 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?
>
>
> Given what you write, it sounds like you already have a good heuristic
> for when to switch from Picard to Newton,
> Hence I think the simplest option is just to use two seperate SNES
> objects - one for performing Picard and one for Newton.
Yes, I considered this option initially. Sometimes it is necessary to
switch back and forth between the methods, so it becomes a bit messy to
organize this in the code.
But maybe if Newton fails after Picard, I should just reduce the time
step and restart, instead of switching back to Picard. Thanks, Dave.
Thanks,
Anton
> The stopping condition for the Picard object would encode your
> heuristic to abort earlier when the solution was deemed sufficiently
> accurate.
>
> Thanks,
> Dave
>
>
> 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/647220e7/attachment-0001.html>
More information about the petsc-users
mailing list