<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class=""><div class=""><br class=""></div><div class=""> Anton,</div><div class=""><br class=""></div> -pc_mg_galerkin pmat<div class=""><br class=""></div><div class=""> Though it seems simple, there is some subtly in swapping out matrices with SNES.</div><div class=""><br class=""></div><div class=""> When using multigrid with SNES there are at least five distinct uses of the Jacobian operator. </div><div class=""><ol class="MailOutline"><ul class=""><li class="">Perform matrix-vector product in line search to check Wolf sufficient decrease convergence criteria</li><li class="">Perform the matrix-vector product for the Krylov accelerator of the system</li><li class="">Perform smoothing on the finest level of MG</li><li class="">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</li><li class="">When using Galerkin products to compute the coarser grid operator performing the Galerkin matrix triple product</li></ul></ol><div class=""><br class=""></div><div class="">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</div><div class=""><br class=""></div><div class="">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.</div><div class=""><br class=""></div><div class="">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.</div><div class=""><br class=""></div><div class="">Barry</div><div class=""><br class=""></div><div class=""><br class=""></div></div><div class=""><div><br class=""><blockquote type="cite" class=""><div class="">On Jun 13, 2021, at 11:21 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class=""><div dir="ltr" class="">On Sun, Jun 13, 2021 at 10:55 AM Anton Popov <<a href="mailto:popov@uni-mainz.de" class="">popov@uni-mainz.de</a>> wrote:<br class=""></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi,<br class="">
<br class="">
I want a simple(?) thing. I want to decide and be able to assign the <br class="">
Jacobian matrix (Amat) on the fly within the FormJacobian function (i.e. <br class="">
during Newton iteration) to one of the following options:<br class="">
<br class="">
1) built-in MFFD approximation<br class="">
2) assembled preconditioner matrix (Pmat)<br class="">
<br class="">
I have not found this scenario demonstrated in any example, therefore <br class="">
I’m asking how to do that.<br class="">
<br class="">
Currently I do the following:<br class="">
<br class="">
1) setup Amat as a shell matrix with a MATOP_MULT operation that simply <br class="">
retrieves a matrix object form its context and calls MatMult on it.<br class="">
<br class="">
2) if I need MFFD, I put a matrix generated with MatCreateSNESMF in the <br class="">
Amat context (of course I also call MatMFFDComputeJacobian before that).<br class="">
<br class="">
3) if I need Pmat, I simply put Pmat in the Amat context.<br class="">
<br class="">
4) call MatAssemblyBegin/End on Amat<br class="">
<br class="">
So far so good.<br class="">
<br class="">
However, shell Amat and assembled Pmat generate a problem if Galerkin <br class="">
multigrid is requested as a preconditioner (I just test on 1 CPU):<br class="">
<br class="">
[0]PETSC ERROR: MatPtAP requires A, shell, to be compatible with P, <br class="">
seqaij (Misses composed function MatPtAP_shell_seqaij_C)<br class="">
[0]PETSC ERROR: #1 MatPtAP()<br class="">
[0]PETSC ERROR: #2 MatGalerkin()<br class="">
[0]PETSC ERROR: #3 PCSetUp_MG()<br class="">
[0]PETSC ERROR: #4 PCSetUp()<br class="">
[0]PETSC ERROR: #5 KSPSetUp()<br class="">
[0]PETSC ERROR: #6 KSPSolve()<br class="">
[0]PETSC ERROR: #7 SNESSolve_NEWTONLS()<br class="">
[0]PETSC ERROR: #8 SNESSolve()<br class="">
<br class="">
It seems like PETSc tries to apply Galerkin coarsening to the shell Amat <br class="">
matrix instead of the assembled Pmat. Why?<br class="">
<br class="">
Pmat is internally generated by SNES using a DMDA attached to SNES, so <br class="">
it has everything to define Galerkin coarsening. And it actually works <br class="">
if I just get rid of the shell Amat.<br class="">
<br class="">
I can get around this problem by wrapping a PC object in a shell matrix <br class="">
and pass it as Pmat. This is a rather ugly approach, so I wonder how to <br class="">
resolve the situation in a better way, if it is possible.<br class=""></blockquote><div class=""><br class=""></div><div class="">Hi Anton,</div><div class=""><br class=""></div><div class="">You are correct that there is no specific test, so I can believe that a bug might be lurking here.</div><div class="">I believe the way it is supposed to work is that you set the type of Galerkin coarsening that you</div><div class="">want</div><div class=""><br class=""></div><div class=""> <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCMGSetGalerkin.html" class="">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCMGSetGalerkin.html</a></div><div class=""><br class=""></div><div class="">So I am thinking you want 'pmat' only, and you would be in charge of making the coarse MFFD operator</div><div class="">and telling PCMG about it. I could imagine that you might want us to automatically do that, but we would</div><div class="">need some way to know that it is MFFD, and with the scheme above we do not have that.</div><div class=""><br class=""></div><div class="">What is the advantage of switching representations during the Newton iteration?</div><div class=""><br class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> Matt</div><div class=""> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Thanks.<br class="">
Anton<br class=""></blockquote></div>-- <br class=""><div dir="ltr" class="gmail_signature"><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a><br class=""></div></div></div></div></div></div></div></div>
</div></blockquote></div><br class=""></div></body></html>