<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none;"> P {margin-top:0;margin-bottom:0;} </style>
</head>
<body 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>Thanks for your help, I really appreciate it. <br>
</div>
<div><br>
</div>
<div>In the end I used a shell matrix to compute the matrix-vector product, it is clearer to me and there are more things under my control. I am now trying to do a parallel implementation, I have some questions on setting up parallel matrices and vectors for
a user-defined partition, could you please provide some advice? Suppose I have already got a partition for 2 CPUs. Each cpu is assigned a list of elements and also their halo elements.<br>
</div>
<div>
<ol>
<li>The global element index for each partition is not necessarily continuous, do I have to I re-order them to make them continuous?<br>
<span></span></li><li><span>When I set up the size of the matrix and vectors for each cpu, should I take into account the halo elements?
<br>
</span></li><li><span>In my serial version, when I initialize my RHS vector, I am not using VecSetValues, Instead I use VecGetArray/VecRestoreArray to assign the values. VecAssemblyBegin()/VecAssemblyEnd() is never used. would this still work for a parallel version?</span></li></ol>
<div>Thanks,</div>
<div>Feng<br>
</div>
</div>
<div><br>
<hr tabindex="-1" style="display:inline-block; width:98%">
<div id="divRplyFwdMsg" dir="ltr"><font style="font-size:11pt" face="Calibri, sans-serif" color="#000000"><b>From:</b> Barry Smith <bsmith@petsc.dev><br>
<b>Sent:</b> 12 March 2021 23:40<br>
<b>To:</b> feng wang <snailsoar@hotmail.com><br>
<b>Cc:</b> petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov><br>
<b>Subject:</b> Re: [petsc-users] Questions on matrix-free GMRES implementation</font>
<div> </div>
</div>
<div class="" style="word-wrap:break-word; line-break:after-white-space"><br class="">
<div><br class="">
<blockquote type="cite" class="">
<div class="">On Mar 12, 2021, at 9:37 AM, feng wang <<a href="mailto:snailsoar@hotmail.com" class="">snailsoar@hotmail.com</a>> wrote:</div>
<br class="x_Apple-interchange-newline">
<div class="">
<div dir="ltr" class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Hi Matt,</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Thanks for your prompt response.
<br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Below are my two versions. one is buggy and the 2nd one is working. For the first one, I add the diagonal contribution to the true RHS (variable: rhs) and then set the base point,
the callback function is somehow called twice afterwards to compute Jacobian. </div>
</div>
</div>
</blockquote>
<div><br class="">
</div>
Do you mean "to compute the Jacobian matrix-vector product?"</div>
<div><br class="">
</div>
<div> Is it only in the first computation of the product (for the given base vector) that it calls it twice or every matrix-vector product?</div>
<div><br class="">
</div>
<div> It is possible there is a bug in our logic; run in the debugger with a break point in <font class="" face="Calibri, Helvetica, sans-serif">FormFunction_mf and each time the function is hit in the debugger type where or bt to get the stack frames from
the calls. Send this. From this we can all see if it is being called excessively and why.</font></div>
<div><br class="">
<blockquote type="cite" class="">
<div class="">
<div dir="ltr" class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">For the 2nd one, I just call the callback function manually to recompute everything, the callback function is then called once as expected to compute the Jacobian. For me, both versions
should do the same things. but I don't know why in the first one the callback function is called twice after I set the base point. what could possibly go wrong?<br class="">
</div>
</div>
</div>
</blockquote>
<div><br class="">
</div>
The logic of how it is suppose to work is shown below.<br class="">
<blockquote type="cite" class="">
<div class="">
<div dir="ltr" class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"></div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Thanks,</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Feng<br class="">
</div>
<div class="">
<div id="x_appendonsend" class=""></div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><b class="">//This does not work</b><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"> fld->cnsv( iqs,iqe, q, aux, csv );
<div class=""> //add contribution of time-stepping</div>
<div class=""> for(iv=0; iv<nv; iv++)</div>
<div class=""> {</div>
<div class=""> for(iq=0; iq<nq; iq++)</div>
<div class=""> {</div>
<div class=""> //use conservative variables here</div>
<div class=""> rhs[iv][iq] = -rhs[iv][iq] + csv[iv][iq]*lhsa[nlhs-1][iq]/cfl;</div>
<div class=""> }</div>
<div class=""> }</div>
<div class=""> ierr = petsc_setcsv(petsc_csv); CHKERRQ(ierr);</div>
<div class=""> ierr = petsc_setrhs(petsc_baserhs); CHKERRQ(ierr);</div>
<div class=""> ierr = MatMFFDSetBase(petsc_A_mf, petsc_csv, petsc_baserhs); CHKERRQ(ierr);</div>
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><b class="">//This works</b><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"> fld->cnsv( iqs,iqe, q, aux, csv );
<div class=""> ierr = petsc_setcsv(petsc_csv); CHKERRQ(ierr);</div>
<div class=""> ierr = FormFunction_mf(this, petsc_csv, petsc_baserhs); //this is my callback function, now call it manually<br class="">
</div>
<div class=""> ierr = MatMFFDSetBase(petsc_A_mf, petsc_csv, petsc_baserhs); CHKERRQ(ierr);</div>
<div class=""><br class="">
</div>
<br class="">
</div>
</div>
</div>
</div>
</blockquote>
<div> Since you provide <span class="" style="font-family:Calibri,Helvetica,sans-serif">petsc_baserhs MatMFFD assumes (naturally) that you will keep the correct values in it. Hence for each new base value YOU need to compute the new values in </span><span class="" style="font-family:Calibri,Helvetica,sans-serif">petsc_baserhs.
This approach gives you a bit more control over reusing the information in </span><span class="" style="font-family:Calibri,Helvetica,sans-serif">petsc_baserhs. </span></div>
<div><span class="" style="font-family:Calibri,Helvetica,sans-serif"><br class="">
</span></div>
<div><span class="" style="font-family:Calibri,Helvetica,sans-serif"> If you would prefer that </span><span class="" style="font-family:Calibri,Helvetica,sans-serif">MatMFFD recomputes the base values, as needed, then you call </span><span class="" style="font-family:Calibri,Helvetica,sans-serif">FormFunction_mf(this,
petsc_csv, NULL); and PETSc will allocate a vector and fill it up as needed by calling your </span><span class="" style="font-family:Calibri,Helvetica,sans-serif">FormFunction_mf() But you need to call MatAssemblyBegin/End each time you the base input vector
</span><span class="" style="font-family:Calibri,Helvetica,sans-serif">this, petsc_csv values change. For </span><font class="" face="Calibri, Helvetica, sans-serif">example</font></div>
<div><font class="" face="Calibri, Helvetica, sans-serif"><br class="">
</font></div>
<div><font class="" face="Calibri, Helvetica, sans-serif"> MatAssemblyBegin(</font><span class="" style="font-family:Calibri,Helvetica,sans-serif">petsc_A_mf,...)</span></div>
<div><span class="" style="font-family:Calibri,Helvetica,sans-serif"> </span><font class="" face="Calibri, Helvetica, sans-serif">MatAssemblyEnd(</font><span class="" style="font-family:Calibri,Helvetica,sans-serif">petsc_A_mf,...)</span></div>
<div><span class="" style="font-family:Calibri,Helvetica,sans-serif"> KSPSolve()</span></div>
<div><span class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</span></div>
<div><span class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</span></div>
<div><br class="">
</div>
<br class="">
<blockquote type="cite" class="">
<div class="">
<div dir="ltr" class="">
<div class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"></div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<hr tabindex="-1" class="" style="display:inline-block; width:98%">
<div id="x_divRplyFwdMsg" dir="ltr" class=""><font class="" style="font-size:11pt" face="Calibri, sans-serif"><b class="">From:</b> Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>><br class="">
<b class="">Sent:</b> 12 March 2021 15:08<br class="">
<b class="">To:</b> feng wang <<a href="mailto:snailsoar@hotmail.com" class="">snailsoar@hotmail.com</a>><br class="">
<b class="">Cc:</b> Barry Smith <<a href="mailto:bsmith@petsc.dev" class="">bsmith@petsc.dev</a>>;
<a href="mailto:petsc-users@mcs.anl.gov" class="">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" class="">petsc-users@mcs.anl.gov</a>><br class="">
<b class="">Subject:</b> Re: [petsc-users] Questions on matrix-free GMRES implementation</font>
<div class=""> </div>
</div>
<div class="">
<div dir="ltr" class="">
<div dir="ltr" class="">On Fri, Mar 12, 2021 at 9:55 AM feng wang <<a href="mailto:snailsoar@hotmail.com" class="">snailsoar@hotmail.com</a>> wrote:<br class="">
</div>
<div class="x_x_gmail_quote">
<blockquote class="x_x_gmail_quote" style="margin:0px 0px 0px 0.8ex; border-left:1px solid rgb(204,204,204); padding-left:1ex">
<div dir="ltr" class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Hi Mat,</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Thanks for your reply. I will try the parallel implementation.</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">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 class=""><br class="">
</div>
<div class="">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 class=""><br class="">
</div>
<div class=""> Thanks,</div>
<div class=""><br class="">
</div>
<div class=""> Matt</div>
<div class=""> </div>
<blockquote class="x_x_gmail_quote" style="margin:0px 0px 0px 0.8ex; border-left:1px solid rgb(204,204,204); padding-left:1ex">
<div dir="ltr" class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Thanks,</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Feng<br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="">
<div id="x_x_gmail-m_-6641396561494424030appendonsend" class=""></div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<hr class="" style="display:inline-block; width:98%">
<div id="x_x_gmail-m_-6641396561494424030divRplyFwdMsg" dir="ltr" class=""><font class="" style="font-size:11pt" face="Calibri, sans-serif"><b class="">From:</b> Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank" class="">knepley@gmail.com</a>><br class="">
<b class="">Sent:</b> 12 March 2021 12:05<br class="">
<b class="">To:</b> feng wang <<a href="mailto:snailsoar@hotmail.com" target="_blank" class="">snailsoar@hotmail.com</a>><br class="">
<b class="">Cc:</b> Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank" class="">bsmith@petsc.dev</a>>;
<a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">petsc-users@mcs.anl.gov</a>><br class="">
<b class="">Subject:</b> Re: [petsc-users] Questions on matrix-free GMRES implementation</font>
<div class=""> </div>
</div>
<div class="">
<div dir="ltr" class="">
<div dir="ltr" class="">On Fri, Mar 12, 2021 at 6:02 AM feng wang <<a href="mailto:snailsoar@hotmail.com" target="_blank" class="">snailsoar@hotmail.com</a>> wrote:<br class="">
</div>
<div class="">
<blockquote class="" style="margin:0px 0px 0px 0.8ex; border-left:1px solid rgb(204,204,204); padding-left:1ex">
<div dir="ltr" class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Hi Barry,</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Thanks for your advice.</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">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 class="">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 class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">In my previous buggy version, I only compute RHS twice. If possible, could you elaborate on your comments "<font class="" face="Calibri, Helvetica, sans-serif">Also be careful about </font><span class="" style="font-family:Calibri,Helvetica,sans-serif"> </span><span class="" 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 class=""><br class="">
</div>
<div class="">Our FD implementation is simple. It approximates the action of the Jacobian as</div>
<div class=""><br class="">
</div>
<div class=""> J(b) v = (F(b + h v) - F(b)) / h ||v||</div>
<div class=""><br class="">
</div>
<div class="">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 class="">and v is the proposed solution update.</div>
<div class=""> </div>
<blockquote class="" style="margin:0px 0px 0px 0.8ex; border-left:1px solid rgb(204,204,204); padding-left:1ex">
<div dir="ltr" class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">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 class=""><br class="">
</div>
<div class="">Sure</div>
<div class=""><br class="">
</div>
<div class=""> <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetSizes.html" target="_blank" class="">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetSizes.html</a></div>
<div class=""><br class="">
</div>
<div class=""> Thanks,</div>
<div class=""><br class="">
</div>
<div class=""> Matt</div>
<div class=""> </div>
<blockquote class="" style="margin:0px 0px 0px 0.8ex; border-left:1px solid rgb(204,204,204); padding-left:1ex">
<div dir="ltr" class="">
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Many thanks,</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt">Feng<br class="">
</div>
<div class="">
<div id="x_x_gmail-m_-6641396561494424030x_gmail-m_-4210113834389506889appendonsend" class="">
</div>
<div class="" style="font-family:Calibri,Helvetica,sans-serif; font-size:12pt"><br class="">
</div>
<hr class="" style="display:inline-block; width:98%">
<div id="x_x_gmail-m_-6641396561494424030x_gmail-m_-4210113834389506889divRplyFwdMsg" dir="ltr" class="">
<font class="" style="font-size:11pt" face="Calibri, sans-serif"><b class="">From:</b> Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank" class="">bsmith@petsc.dev</a>><br class="">
<b class="">Sent:</b> 11 March 2021 22:15<br class="">
<b class="">To:</b> feng wang <<a href="mailto:snailsoar@hotmail.com" target="_blank" class="">snailsoar@hotmail.com</a>><br class="">
<b class="">Cc:</b> <a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">
petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">petsc-users@mcs.anl.gov</a>><br class="">
<b class="">Subject:</b> Re: [petsc-users] Questions on matrix-free GMRES implementation</font>
<div class=""> </div>
</div>
<div class="">
<div class=""><br class="">
</div>
Feng,
<div class=""><br class="">
</div>
<div class=""> 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 class=""><br class="">
</div>
<div class="">The easiest way is to call <font class="" 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 class="" style="font-family:Calibri,Helvetica,sans-serif"> </span><span class="" 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 class="" style="font-family:Calibri,Helvetica,sans-serif">petsc_baserhs values to the function evaluation at that point.</span></div>
<div class=""><span class="" style="font-family:Calibri,Helvetica,sans-serif"><br class="">
</span></div>
<div class=""><span class="" 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 class=""><br class="">
</div>
<div class="">Barry</div>
<div class=""><br class="">
</div>
<div class="">
<div class=""><br class="">
<blockquote type="cite" class="">
<div class="">On Mar 11, 2021, at 7:35 AM, feng wang <<a href="mailto:snailsoar@hotmail.com" target="_blank" class="">snailsoar@hotmail.com</a>> wrote:</div>
<br class="">
<div class="">
<div class="" 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 class="" 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 class="">
</div>
<div class="" 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 class="" 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 class="">
</div>
<div class="" 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 class="" 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 class="">
</div>
<div class="" 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 class=""> ierr = MatMFFDSetFunction(petsc_A_mf, FormFunction_mf, this); CHKERRQ(ierr);</div>
<div class=""></div>
<div class=""><br class="">
</div>
<div class="">KSP linear operator is set up as:</div>
<div class=""><br class="">
</div>
<div class=""> ierr = KSPSetOperators(petsc_ksp, petsc_A_mf, petsc_A_pre); CHKERRQ(ierr); //petsc_A_pre is my assembled pre-conditioning matrix<br class="">
</div>
<div class=""><br class="">
</div>
<div class="">Before calling KSPSolve, I do:</div>
<div class=""><br class="">
</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 class="" 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 class="">
</div>
<div class="" 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 class="" 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 class="">
</div>
<div class="" 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 class=""> {</div>
<div class=""> PetscErrorCode ierr;</div>
<div class=""> cFdDomain *user_ctx;</div>
<div class=""><br class="">
</div>
<div class=""> cout << "FormFunction_mf called\n";</div>
<div class=""><br class="">
</div>
<div class=""> //in_vec: flow states</div>
<div class=""> //out_vec: right hand side + diagonal contributions from CFL number</div>
<div class=""><br class="">
</div>
<div class=""> user_ctx = (cFdDomain*)ctx;</div>
<div class=""><br class="">
</div>
<div class=""> //get perturbed conservative variables from petsc</div>
<div class=""> user_ctx->petsc_getcsv(in_vec);</div>
<div class=""><br class="">
</div>
<div class=""> //get new right side</div>
<div class=""> user_ctx->petsc_fd_rhs();</div>
<div class=""><br class="">
</div>
<div class=""> //set new right hand side to the output vector</div>
<div class=""> user_ctx->petsc_setrhs(out_vec);</div>
<div class=""><br class="">
</div>
<div class=""> ierr = 0;</div>
<div class=""> return ierr;</div>
<div class=""> }</div>
</div>
<div class="" 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 class="">
</div>
<div class="" 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 class=""> </span><br class="">
</div>
<div class="" 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 class="">
</div>
<div class="" 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 class="" 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 class="">
</div>
<div class="" 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 class="" 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 class=""> </span></div>
</div>
</blockquote>
</div>
<br class="">
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br class="" clear="all">
<div class=""><br class="">
</div>
-- <br class="">
<div dir="ltr" class="">
<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>
</div>
</div>
</blockquote>
</div>
<br class="" clear="all">
<div class=""><br class="">
</div>
-- <br class="">
<div dir="ltr" class="x_x_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>
</div>
</div>
</div>
</blockquote>
</div>
<br class="">
</div>
</div>
</body>
</html>