<div dir="ltr"><div dir="ltr">On Fri, Mar 12, 2021 at 9:55 AM feng wang <<a href="mailto:snailsoar@hotmail.com">snailsoar@hotmail.com</a>> wrote:<br></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">
<div dir="ltr">
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Hi Mat,</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Thanks for your reply. I will try the parallel implementation.</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
I've got a serial matrix-free GMRES working, but I would like to know why my initial version of matrix-free implementation does not work and there is still something I don't understand. I did some debugging and find that the callback function to compute the
RHS for the matrix-free matrix is called twice by Petsc when it computes the finite difference Jacobian, but it should only be called once. I don't know why, could you please give some advice?</div></div></blockquote><div><br></div><div>F is called once to calculate the base point and once to get the perturbation. The base point is not recalculated, so if you do many iterates, it is amortized.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Thanks,</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Feng<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<div>
<div id="gmail-m_-6641396561494424030appendonsend"></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<hr style="display:inline-block;width:98%">
<div id="gmail-m_-6641396561494424030divRplyFwdMsg" dir="ltr"><font style="font-size:11pt" face="Calibri, sans-serif" color="#000000"><b>From:</b> Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>
<b>Sent:</b> 12 March 2021 12:05<br>
<b>To:</b> feng wang <<a href="mailto:snailsoar@hotmail.com" target="_blank">snailsoar@hotmail.com</a>><br>
<b>Cc:</b> Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>>; <a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b> Re: [petsc-users] Questions on matrix-free GMRES implementation</font>
<div> </div>
</div>
<div>
<div dir="ltr">
<div dir="ltr">On Fri, Mar 12, 2021 at 6:02 AM feng wang <<a href="mailto:snailsoar@hotmail.com" target="_blank">snailsoar@hotmail.com</a>> wrote:<br>
</div>
<div>
<blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr">
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Hi Barry,</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Thanks for your advice.</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
You are right on this. somehow there is some inconsistency when I compute the right hand side (true RHS + time-stepping contribution to the diagonal matrix) to compute the finite difference Jacobian. If I just use the call back function to recompute my RHS
before I call <b>MatMFFDSetBase</b>, then it works like a charm. But now I end up with computing my RHS three times. 1st time is to compute the true RHS, the rest two is for computing finite difference Jacobian.
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
In my previous buggy version, I only compute RHS twice. If possible, could you elaborate on your comments "<font face="Calibri, Helvetica, sans-serif">Also be careful about </font><span style="font-family:Calibri,Helvetica,sans-serif"> </span><span style="font-family:Calibri,Helvetica,sans-serif">petsc_baserhs</span>",
so I may possibly understand what was going on with my buggy version.</div>
</div>
</blockquote>
<div><br>
</div>
<div>Our FD implementation is simple. It approximates the action of the Jacobian as</div>
<div><br>
</div>
<div> J(b) v = (F(b + h v) - F(b)) / h ||v||</div>
<div><br>
</div>
<div>where h is some small parameter and b is the base vector, namely the one that you are linearizing around. In a Newton step, b is the previous solution</div>
<div>and v is the proposed solution update.</div>
<div> </div>
<blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr">
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Besides, for a parallel implementation, my code already has its own partition method, is it possible to allow petsc read in a user-defined partition? if not what is a better way to do this?</div>
</div>
</blockquote>
<div><br>
</div>
<div>Sure</div>
<div><br>
</div>
<div> <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetSizes.html" target="_blank">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetSizes.html</a></div>
<div><br>
</div>
<div> Thanks,</div>
<div><br>
</div>
<div> Matt</div>
<div> </div>
<blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr">
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Many thanks,</div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
Feng<br>
</div>
<div>
<div id="gmail-m_-6641396561494424030x_gmail-m_-4210113834389506889appendonsend"></div>
<div style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt;color:rgb(0,0,0)">
<br>
</div>
<hr style="display:inline-block;width:98%">
<div id="gmail-m_-6641396561494424030x_gmail-m_-4210113834389506889divRplyFwdMsg" dir="ltr"><font style="font-size:11pt" face="Calibri, sans-serif" color="#000000"><b>From:</b> Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>><br>
<b>Sent:</b> 11 March 2021 22:15<br>
<b>To:</b> feng wang <<a href="mailto:snailsoar@hotmail.com" target="_blank">snailsoar@hotmail.com</a>><br>
<b>Cc:</b> <a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b> Re: [petsc-users] Questions on matrix-free GMRES implementation</font>
<div> </div>
</div>
<div>
<div><br>
</div>
Feng,
<div><br>
</div>
<div> The first thing to check is that for each linear solve that involves a new operator (values in the base vector) the MFFD matrix knows it is using a new operator. </div>
<div><br>
</div>
<div>The easiest way is to call <font face="Calibri, Helvetica, sans-serif">MatMFFDSetBase() before each solve that involves a new operator (new values in the base vector). Also be careful about </font><span style="font-family:Calibri,Helvetica,sans-serif"> </span><span style="font-family:Calibri,Helvetica,sans-serif">petsc_baserhs,
when you change the base vector's values you also need to change the </span><span style="font-family:Calibri,Helvetica,sans-serif">petsc_baserhs values to the function evaluation at that point.</span></div>
<div><span style="font-family:Calibri,Helvetica,sans-serif"><br>
</span></div>
<div><span style="font-family:Calibri,Helvetica,sans-serif">If that is correct I would check with a trivial function evaluator to make sure the infrastructure is all set up correctly. For examples use for the matrix free a 1 4 1 operator applied matrix free. </span></div>
<div><br>
</div>
<div>Barry</div>
<div><br>
</div>
<div>
<div><br>
<blockquote type="cite">
<div>On Mar 11, 2021, at 7:35 AM, feng wang <<a href="mailto:snailsoar@hotmail.com" target="_blank">snailsoar@hotmail.com</a>> wrote:</div>
<br>
<div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
Dear All,</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
I am new to petsc and trying to implement a matrix-free GMRES. I have assembled an approximate Jacobian matrix just for preconditioning. After reading some previous questions on this topic, my approach is:</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
the matrix-free matrix is created as:</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
ierr = MatCreateMFFD(*A_COMM_WORLD, iqe*blocksize, iqe*blocksize, PETSC_DETERMINE, PETSC_DETERMINE, &petsc_A_mf); CHKERRQ(ierr);
<div> ierr = MatMFFDSetFunction(petsc_A_mf, FormFunction_mf, this); CHKERRQ(ierr);</div>
<div></div>
<div><br>
</div>
<div>KSP linear operator is set up as:</div>
<div><br>
</div>
<div> ierr = KSPSetOperators(petsc_ksp, petsc_A_mf, petsc_A_pre); CHKERRQ(ierr); //petsc_A_pre is my assembled pre-conditioning matrix<br>
</div>
<div><br>
</div>
<div>Before calling KSPSolve, I do:</div>
<div><br>
</div>
ierr = MatMFFDSetBase(petsc_A_mf, petsc_csv, petsc_baserhs); CHKERRQ(ierr); //petsc_csv is the flow states, petsc_baserhs is the pre-computed right hand side</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
The call back function is defined as:</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
PetscErrorCode cFdDomain::FormFunction_mf(void *ctx, Vec in_vec, Vec out_vec)
<div> {</div>
<div> PetscErrorCode ierr;</div>
<div> cFdDomain *user_ctx;</div>
<div><br>
</div>
<div> cout << "FormFunction_mf called\n";</div>
<div><br>
</div>
<div> //in_vec: flow states</div>
<div> //out_vec: right hand side + diagonal contributions from CFL number</div>
<div><br>
</div>
<div> user_ctx = (cFdDomain*)ctx;</div>
<div><br>
</div>
<div> //get perturbed conservative variables from petsc</div>
<div> user_ctx->petsc_getcsv(in_vec);</div>
<div><br>
</div>
<div> //get new right side</div>
<div> user_ctx->petsc_fd_rhs();</div>
<div><br>
</div>
<div> //set new right hand side to the output vector</div>
<div> user_ctx->petsc_setrhs(out_vec);</div>
<div><br>
</div>
<div> ierr = 0;</div>
<div> return ierr;</div>
<div> }</div>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
The linear system I am solving is (J+D)x=RHS. J is the Jacobian matrix. D is a diagonal matrix and it is used to stabilise the solution at the start but reduced gradually when the solution moves on to recover Newton's method. I add D*x to the true right side
when non-linear function is computed to work out finite difference Jacobian, so when finite difference is used, it actually computes (J+D)*dx.<span> </span><br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
The code runs but diverges in the end. If I don't do matrix-free and use my approximate Jacobian matrix, GMRES works. So something is wrong with my matrix-free implementation. Have I missed something in my implementation? Besides, is there a way to check
if the finite difference Jacobian matrix is computed correctly in a matrix-free implementation?</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
<br>
</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
Thanks for your help in advance.</div>
<div style="font-style:normal;font-variant-caps:normal;font-weight:normal;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;font-family:Calibri,Helvetica,sans-serif;font-size:12pt">
Feng <span> </span></div>
</div>
</blockquote>
</div>
<br>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
-- <br>
<div dir="ltr">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener</div>
<div><br>
</div>
<div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>