[petsc-users] Context behind SNESComputeFunction call
Smith, Barry F.
bsmith at mcs.anl.gov
Fri Jan 26 23:47:48 CST 2018
Thanks, I now understand the situation.
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.
Please let me know if it works for you and what you think,
Barry
> On Jan 26, 2018, at 5:54 PM, Alexander Lindsay <alexlindsay239 at gmail.com> wrote:
>
> 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.
>
> 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.
>
> On Jan 26, 2018, at 5:34 PM, Kong, Fande <fande.kong at inl.gov> wrote:
>
>> Hi Barry,
>>
>> I made minor changes on src/snes/examples/tutorials/ex2.c to demonstrate this issue. Please see the attachment.
>>
>> ./ex2 -snes_monitor -ksp_monitor -snes_mf_operator 1
>>
>>
>> atol=1e-50, rtol=1e-08, stol=1e-08, maxit=50, maxf=10000
>> FormFunction is called
>> 0 SNES Function norm 5.414682427127e+00
>> 0 KSP Residual norm 9.559082033938e-01
>> FormFunction is called
>> FormFunction is called
>> 1 KSP Residual norm 1.703870633386e-09
>> FormFunction is called
>> FormFunction is called
>> 1 SNES Function norm 2.952582481151e-01
>> 0 KSP Residual norm 2.672054855433e-02
>> FormFunction is called
>> FormFunction is called
>> 1 KSP Residual norm 1.519298012177e-10
>> FormFunction is called
>> FormFunction is called
>> 2 SNES Function norm 4.502289047587e-04
>> 0 KSP Residual norm 4.722075651268e-05
>> FormFunction is called
>> FormFunction is called
>> 1 KSP Residual norm 3.834927363659e-14
>> FormFunction is called
>> FormFunction is called
>> 3 SNES Function norm 1.390073376131e-09
>> number of SNES iterations = 3
>>
>> Norm of error 1.49795e-10, Iterations 3
>>
>> "FormFunction" is called TWICE at "0 KSP".
>>
>> If we comment out MatMFFDSetFunction:
>>
>> /* ierr = MatMFFDSetFunction(Jacobian,FormFunction_MFFD,(void*)snes);CHKERRQ(ierr); */
>>
>>
>> ./ex2 -snes_monitor -ksp_monitor -snes_mf_operator 1
>>
>> atol=1e-50, rtol=1e-08, stol=1e-08, maxit=50, maxf=10000
>> FormFunction is called
>> 0 SNES Function norm 5.414682427127e+00
>> 0 KSP Residual norm 9.559082033938e-01
>> FormFunction is called
>> 1 KSP Residual norm 1.703870633386e-09
>> FormFunction is called
>> FormFunction is called
>> 1 SNES Function norm 2.952582481151e-01
>> 0 KSP Residual norm 2.672054855433e-02
>> FormFunction is called
>> 1 KSP Residual norm 1.519298012177e-10
>> FormFunction is called
>> FormFunction is called
>> 2 SNES Function norm 4.502289047587e-04
>> 0 KSP Residual norm 4.722075651268e-05
>> FormFunction is called
>> 1 KSP Residual norm 3.834927363659e-14
>> FormFunction is called
>> FormFunction is called
>> 3 SNES Function norm 1.390073376131e-09
>> number of SNES iterations = 3
>>
>> Norm of error 1.49795e-10, Iterations 3
>>
>> "FormFunction" is called ONCE at "0 KSP".
>>
>> Hopefully, this example makes the point clear.
>>
>>
>> Fande,
>>
>> On Fri, Jan 26, 2018 at 3:50 PM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
>>
>>
>> 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.
>>
>> Barry
>>
>>
>> > On Jan 26, 2018, at 4:32 PM, Kong, Fande <fande.kong at inl.gov> wrote:
>> >
>> >
>> >
>> > On Fri, Jan 26, 2018 at 3:10 PM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
>> >
>> >
>> > > On Jan 26, 2018, at 2:15 PM, Kong, Fande <fande.kong at inl.gov> wrote:
>> > >
>> > >
>> > >
>> > > On Mon, Jan 8, 2018 at 2:15 PM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
>> > >
>> > >
>> > > > On Jan 8, 2018, at 2:59 PM, Alexander Lindsay <alexlindsay239 at gmail.com> wrote:
>> > > >
>> > > > Is there any elegant way to tell whether SNESComputeFunction is being called under different conceptual contexts?
>> > > >
>> > > > E.g. non-linear residual evaluation vs. Jacobian formation from finite differencing vs. Jacobian-vector products from finite differencing?
>> > >
>> > > Under normal usage with the options database no.
>> > >
>> > > Hi Barry,
>> > >
>> > > How difficult to provide an API? Is it possible?
>> > >
>> > >
>> > >
>> > > 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.
>> > >
>> > > 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?
>> >
>> > 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
>> >
>> > /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
>> > if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
>> > ierr = (*ctx->func)(ctx->funcctx,U,F);CHKERRQ(ierr);
>> > }
>> >
>> > since the if is satisfied it does not compute the function at the base location. To double check I ran src/snes/examples/tutorials/ex19 with -snes_mf in the debugger and verified that the "extra" function evaluations are not done.
>> >
>> > In MatAssemblyEnd_SNESMF,
>> >
>> > if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction) {
>> > ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr);
>> > ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr);
>> > } else {
>> > /* f value known by SNES is not correct for other differencing function */
>> > ierr = MatMFFDSetBase_MFFD(J,u,NULL);CHKERRQ(ierr);
>> > }
>> >
>> >
>> > Will hit ierr = MatMFFDSetBase_MFFD(J,u,NULL);CHKERRQ(ierr), because SNES and MAT have different function pointers.
>> >
>> > In MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F),
>> >
>> > if (F) {
>> > if (ctx->current_f_allocated) {ierr = VecDestroy(&ctx->current_f);CHKERRQ(ierr);}
>> > ctx->current_f = F;
>> > ctx->current_f_allocated = PETSC_FALSE;
>> > } else if (!ctx->current_f_allocated) {
>> > ierr = MatCreateVecs(J,NULL,&ctx->current_f);CHKERRQ(ierr);
>> >
>> > ctx->current_f_allocated = PETSC_TRUE;
>> > }
>> >
>> > Because F=NULL, we then have ctx->current_f_allocated = PETSC_TRUE.
>> >
>> > Then, the following if statement is true:
>> >
>> > if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
>> > ierr = (*ctx->func)(ctx->funcctx,U,F);CHKERRQ(ierr);
>> > }
>> >
>> >
>> > Fande,
>> >
>> >
>> >
>> > Barry
>> >
>> >
>> > >
>> > > Fande,
>> > >
>> > >
>> > >
>> > > Barry
>> > >
>> > >
>> > >
>> > > >
>> > > > Alex
>> > >
>> > >
>>
>>
>> <ex2.c>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180127/148c925f/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex3.c
Type: application/octet-stream
Size: 14125 bytes
Desc: ex3.c
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180127/148c925f/attachment-0001.obj>
More information about the petsc-users
mailing list