<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hi Barry & Matt,<br>
<br>
thanks for your quick response. These options were exactly what I
needed and expected:<br>
<br>
-pc_mg_galerkin pmat<br>
-pc_use_amat false<br>
<br>
I just assumed that it’s a default behavior of the PC object.<br>
<br>
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.<br>
<br>
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.<br>
<br>
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.<br>
<br>
Is my workaround with the shell matrix acceptable, or there is a
better solution? <br>
<br>
Thanks,<br>
Anton <br>
</p>
<div class="moz-cite-prefix">On 13.06.21 20:52, Barry Smith wrote:<br>
</div>
<blockquote type="cite"
cite="mid:D216B541-F44D-49E1-A6B0-3FFF7618C7C9@petsc.dev">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<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=""
moz-do-not-send="true">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="" moz-do-not-send="true">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="" moz-do-not-send="true">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=""
moz-do-not-send="true">https://www.cse.buffalo.edu/~knepley/</a><br
class="">
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br class="">
</div>
</blockquote>
</body>
</html>