<div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Wed, Sep 5, 2018 at 9:34 PM Oleksandr Koshkarov <<a href="mailto:olk548@mail.usask.ca">olk548@mail.usask.ca</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello all,<br>
<br>
I got lost a little with petsc and want to clarify something.<br>
<br>
I want to do my nonlinear solves with matrix free method, but I want a <br>
real matrix preconditioner and invert it with LU.<br>
<br>
without preconditioner it is simple:<br>
<br>
Mat Jmf;<br>
<br>
MatCreateSNESMF(snes,&Jmf);<br>
<br>
SNESSetJacobian(snes,Jmf,Jmf,MatMFFDComputeJacobian,0);<br>
<br>
<br>
But with preconditioner, I am lost:<br>
<br>
Mat Jmf,J;<br>
<br>
MatCreateSNESMF(snes,&Jmf);<br>
<br>
DMCreateMatrix(da,&J);<br>
<br>
SNESSetJacobian(snes,Jmf,J,FromJacobian,0);<br>
<br>
//And here I have problems:<br>
<br>
PetscErrorCode FormJacobian(SNES snes,Vec X,Mat Jmf,Mat J,void *ptr){<br>
<br>
// I am constructing it only once in a while, not every nonlinear solve <br>
(this snes is inside time stepper).<br>
<br>
if (sometimes) {  MatSetValuesStencil(J,...); ...    }<br>
<br>
   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);<br>
   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);<br>
<br>
// assembling Jmf<br>
   MatAssemblyBegin(Jmf,MAT_FINAL_ASSEMBLY);<br>
   MatAssemblyEnd(Jmf,MAT_FINAL_ASSEMBLY);<br>
<br>
   return(0);<br>
}<br>
<br>
Here I have two concerns:<br>
<br>
1) I am changing J only once in a while, because it is expansive to do <br>
otherwise. Am I doing it correctly? During iteration when I am not <br>
changing J, will petsc know this?</blockquote><div><br></div><div>Yes, every matrix has a state variable that is updated when values are changed.</div><div>The LU looks at this state variable to decide whether to refactor.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> I do not want petsc to perform LU <br>
every time when I am not changing J. And when I will change J, will it <br>
work?</blockquote><div><br></div><div>Yes, but it is easy to check in the log. Use -log_view and look at the number of calls</div><div>to PCLUFactorNumeric.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> (My snes is inside time stepper, and I want to construct <br>
preconditioner only say every 100 time steps)<br>
<br>
2) What about Jmf? how can I let petsc know that it should use matrix <br>
free finite difference? Will it work if I let it as it is? (I think so, <br>
but I am not sure - my information source is <br>
ts/examples/tutorials/ex15.c, but I do not understand it fully)<br></blockquote><div><br></div><div>It appears you have it correct. I always use MMS to make sure my solver is doing what I expect.</div><div><br></div><div>  Thanks,</div><div><br></div><div>    Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Thank you and best regards,<br>
<br>
Alex.<br>
<br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature" data-smartmail="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>