<div dir="ltr">Hi Barry,<div><br></div><div>Thanks so much! </div><div><br></div><div>Does this requires us to attach a shell function "MatAssemblyEnd_SNESMF"? We need to maintain "MatAssemblyEnd_SNESMF" in the moose side?</div><div><br></div><div><br></div><div><br></div><div>Could we just add a flag into the original routine "MatAssemblyEnd_SNESMF" in PETSc?</div><div><br></div><div><br></div><div><div><i>static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt)</i></div><div><i>{</i></div><div><i> PetscErrorCode ierr;</i></div><div><i> MatMFFD j = (MatMFFD)J->data;</i></div><div><i> SNES snes = (SNES)j->ctx;</i></div><div><i> PetscBool refuse;</i></div><div><i> Vec u,f;</i></div><div><i><br></i></div><div><i> PetscFunctionBegin;</i></div><div><i> ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr);</i></div><div><i> </i></div><div><i> <font color="#ff0000">ierr = SNESGetReuseBaseMFFD(SNES snes, & refuse);CHKERRQ(ierr);</font></i></div><div><i><br></i></div><div><i> ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr);</i></div><div><i> if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction <font color="#ff0000">|| reuse</font>) {</i></div><div><i> ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr);</i></div><div><i> ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr);</i></div><div><i> } else {</i></div><div><i> /* f value known by SNES is not correct for other differencing function */</i></div><div><i> ierr = MatMFFDSetBase_MFFD(J,u,NULL);CHKERRQ(ierr);</i></div><div><i> }</i></div><div><i> PetscFunctionReturn(0);</i></div><div><i>}</i></div></div><div><br></div><div><br></div><div><br></div><div>We may need new APIs <span style="color:rgb(255,0,0)">SNESGetReuseBaseMFFD </span><font color="#000000">and</font><span style="color:rgb(255,0,0)"> </span><font color="#ff0000">SNESSetReuseBaseMFFD,</font><font color="#000000"> and a command line option "-</font><font color="#ff0000">snes_reuse_base_mffd</font><font color="#000000">". "reuse" is false by default, and it is consistent with current setting.</font></div><div><br></div><div><br></div><div>Please let me know what do you think. If you agree, I can make a PR on this.</div><div><br></div><div>Fande,</div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Jan 26, 2018 at 10:47 PM, Smith, Barry F. <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>
<div class="m_-3841611664697760969BodyFragment"><font size="2"><span style="font-size:11pt">
<div class="m_-3841611664697760969PlainText"><br>
Thanks, I now understand the situation.<br>
<br>
I have a tentative solution for you that does not require complex new APIs. I have added one function to the master branch MatSNESMFGetSNES() and attach a modified example that utilizes this to avoid the extra function evaluations.
<br>
<br>
</div>
</span></font></div>
<div class="m_-3841611664697760969BodyFragment"><font size="2"><span style="font-size:11pt">
<div class="m_-3841611664697760969PlainText"><br>
<br>
Please let me know if it works for you and what you think,<br>
<br>
Barry<div><div class="h5"><br>
<br>
<br>
> On Jan 26, 2018, at 5:54 PM, Alexander Lindsay <<a href="mailto:alexlindsay239@gmail.com" target="_blank">alexlindsay239@gmail.com</a>> wrote:<br>
> <br>
> We are doing something non-standard. We set different functions for snes and mffd. When a function is called through snes, we know we are doing a non-linear residual evaluation and we allow an update of our mechanical contact nodes. When a function is called
through mffd, we know we are within the linear solve, and we do not allow our contact state to change. If after a contact update our non-linear residual is not deemed to be converged, then we would hope to be able to re-use it for the zeroth-linear residual.
<br>
> <br>
> I know that you don't like this way of handling contact, but it integrates well for us into our multiphysics framework and works well in many cases.
<br>
> <br>
> On Jan 26, 2018, at 5:34 PM, Kong, Fande <<a href="mailto:fande.kong@inl.gov" target="_blank">fande.kong@inl.gov</a>> wrote:<br>
> <br>
>> Hi Barry,<br>
>> <br>
>> I made minor changes on src/snes/examples/tutorials/<wbr>ex2.c to demonstrate this issue. Please see the attachment.<br>
>> <br>
>> ./ex2 -snes_monitor -ksp_monitor -snes_mf_operator 1<br>
>> <br>
>> <br>
>> atol=1e-50, rtol=1e-08, stol=1e-08, maxit=50, maxf=10000<br>
>> FormFunction is called <br>
>> 0 SNES Function norm 5.414682427127e+00 <br>
>> 0 KSP Residual norm 9.559082033938e-01 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 1 KSP Residual norm 1.703870633386e-09 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 1 SNES Function norm 2.952582481151e-01 <br>
>> 0 KSP Residual norm 2.672054855433e-02 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 1 KSP Residual norm 1.519298012177e-10 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 2 SNES Function norm 4.502289047587e-04 <br>
>> 0 KSP Residual norm 4.722075651268e-05 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 1 KSP Residual norm 3.834927363659e-14 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 3 SNES Function norm 1.390073376131e-09 <br>
>> number of SNES iterations = 3<br>
>> <br>
>> Norm of error 1.49795e-10, Iterations 3<br>
>> <br>
>> "FormFunction" is called TWICE at "0 KSP". <br>
>> <br>
>> If we comment out MatMFFDSetFunction:<br>
>> <br>
>> /* ierr = MatMFFDSetFunction(Jacobian,<wbr>FormFunction_MFFD,(void*)snes)<wbr>;CHKERRQ(ierr); */<br>
>> <br>
>> <br>
>> ./ex2 -snes_monitor -ksp_monitor -snes_mf_operator 1<br>
>> <br>
>> atol=1e-50, rtol=1e-08, stol=1e-08, maxit=50, maxf=10000<br>
>> FormFunction is called <br>
>> 0 SNES Function norm 5.414682427127e+00 <br>
>> 0 KSP Residual norm 9.559082033938e-01 <br>
>> FormFunction is called <br>
>> 1 KSP Residual norm 1.703870633386e-09 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 1 SNES Function norm 2.952582481151e-01 <br>
>> 0 KSP Residual norm 2.672054855433e-02 <br>
>> FormFunction is called <br>
>> 1 KSP Residual norm 1.519298012177e-10 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 2 SNES Function norm 4.502289047587e-04 <br>
>> 0 KSP Residual norm 4.722075651268e-05 <br>
>> FormFunction is called <br>
>> 1 KSP Residual norm 3.834927363659e-14 <br>
>> FormFunction is called <br>
>> FormFunction is called <br>
>> 3 SNES Function norm 1.390073376131e-09 <br>
>> number of SNES iterations = 3<br>
>> <br>
>> Norm of error 1.49795e-10, Iterations 3<br>
>> <br>
>> "FormFunction" is called ONCE at "0 KSP". <br>
>> <br>
>> Hopefully, this example makes the point clear.<br>
>> <br>
>> <br>
>> Fande,<br>
>> <br>
>> On Fri, Jan 26, 2018 at 3:50 PM, Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
>> <br>
>> <br>
>> So you are doing something non-standard? Are you not just using -snes_mf or -snes_mf_operator? Can you send me a sample code that has the extra function evaluations? Because if you run through regular usage with the debugger you will see there is no extra
evaluation.<br>
>> <br>
>> Barry<br>
>> <br>
>> <br>
>> > On Jan 26, 2018, at 4:32 PM, Kong, Fande <<a href="mailto:fande.kong@inl.gov" target="_blank">fande.kong@inl.gov</a>> wrote:<br>
>> ><br>
>> ><br>
>> ><br>
>> > On Fri, Jan 26, 2018 at 3:10 PM, Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
>> ><br>
>> ><br>
>> > > On Jan 26, 2018, at 2:15 PM, Kong, Fande <<a href="mailto:fande.kong@inl.gov" target="_blank">fande.kong@inl.gov</a>> wrote:<br>
>> > ><br>
>> > ><br>
>> > ><br>
>> > > On Mon, Jan 8, 2018 at 2:15 PM, Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
>> > ><br>
>> > ><br>
>> > > > On Jan 8, 2018, at 2:59 PM, Alexander Lindsay <<a href="mailto:alexlindsay239@gmail.com" target="_blank">alexlindsay239@gmail.com</a>> wrote:<br>
>> > > ><br>
>> > > > Is there any elegant way to tell whether SNESComputeFunction is being called under different conceptual contexts?<br>
>> > > ><br>
>> > > > E.g. non-linear residual evaluation vs. Jacobian formation from finite differencing vs. Jacobian-vector products from finite differencing?<br>
>> > ><br>
>> > > Under normal usage with the options database no.<br>
>> > ><br>
>> > > Hi Barry,<br>
>> > ><br>
>> > > How difficult to provide an API? Is it possible?<br>
>> > ><br>
>> > ><br>
>> > ><br>
>> > > If you have some reason to know you could write three functions and provide them to SNESSetFunction(), MatMFFDSetFunction(), and MatFDColoringSetFunction(). Note that these functions have slightly different calling sequences but you can have all of
them call the same underlying common function if you are only interested in, for example, how many times the function is used for each purpose.<br>
>> > ><br>
>> > > If we use this way for the Jacobian-free Newton, the function evaluation will be called twice at the first linear iteration because the computed residual vector at the nonlinear step is not reused. Any way to reuse the function residual of the Newton
step instead of recomputing a new residual at the first linear iteration?<br>
>> ><br>
>> > It does reuse the function evaluation. Why do you think it does not? If you look at MatMult_MFFD() you will see the lines of code<br>
>> ><br>
>> > /* compute func(U) as base for differencing; only needed first time in and not when provided by user */<br>
>> > if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {<br>
>> > ierr = (*ctx->func)(ctx->funcctx,U,F)<wbr>;CHKERRQ(ierr);<br>
>> > }<br>
>> ><br>
>> > since the if is satisfied it does not compute the function at the base location. To double check I ran src/snes/examples/tutorials/<wbr>ex19 with -snes_mf in the debugger and verified that the "extra" function evaluations are not done.<br>
>> ><br>
>> > In MatAssemblyEnd_SNESMF,<br>
>> ><br>
>> > if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))<wbr>SNESComputeFunction) {<br>
>> > ierr = SNESGetFunction(snes,&f,NULL,<wbr>NULL);CHKERRQ(ierr);<br>
>> > ierr = MatMFFDSetBase_MFFD(J,u,f);<wbr>CHKERRQ(ierr);<br>
>> > } else {<br>
>> > /* f value known by SNES is not correct for other differencing function */<br>
>> > ierr = MatMFFDSetBase_MFFD(J,u,NULL);<wbr>CHKERRQ(ierr);<br>
>> > }<br>
>> ><br>
>> ><br>
>> > Will hit ierr = MatMFFDSetBase_MFFD(J,u,NULL);<wbr>CHKERRQ(ierr), because SNES and MAT have different function pointers.<br>
>> ><br>
>> > In MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F),<br>
>> ><br>
>> > if (F) {<br>
>> > if (ctx->current_f_allocated) {ierr = VecDestroy(&ctx->current_f);<wbr>CHKERRQ(ierr);}<br>
>> > ctx->current_f = F;<br>
>> > ctx->current_f_allocated = PETSC_FALSE;<br>
>> > } else if (!ctx->current_f_allocated) {<br>
>> > ierr = MatCreateVecs(J,NULL,&ctx-><wbr>current_f);CHKERRQ(ierr);<br>
>> ><br>
>> > ctx->current_f_allocated = PETSC_TRUE;<br>
>> > }<br>
>> ><br>
>> > Because F=NULL, we then have ctx->current_f_allocated = PETSC_TRUE.<br>
>> ><br>
>> > Then, the following if statement is true:<br>
>> ><br>
>> > if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {<br>
>> > ierr = (*ctx->func)(ctx->funcctx,U,F)<wbr>;CHKERRQ(ierr);<br>
>> > }<br>
>> ><br>
>> ><br>
>> > Fande,<br>
>> ><br>
>> ><br>
>> ><br>
>> > Barry<br>
>> ><br>
>> ><br>
>> > ><br>
>> > > Fande,<br>
>> > ><br>
>> > ><br>
>> > ><br>
>> > > Barry<br>
>> > ><br>
>> > ><br>
>> > ><br>
>> > > ><br>
>> > > > Alex<br>
>> > ><br>
>> > ><br>
>> <br>
>> <br>
>> <ex2.c><br>
<br>
</div></div></div>
</span></font></div>
</div>
</blockquote></div><br></div>