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 <alexlindsay...@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.k...@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. <bsm...@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.k...@inl.gov> wrote:
>> >
>> >
>> >
>> > On Fri, Jan 26, 2018 at 3:10 PM, Smith, Barry F. <bsm...@mcs.anl.gov> 
>> > wrote:
>> >
>> >
>> > > On Jan 26, 2018, at 2:15 PM, Kong, Fande <fande.k...@inl.gov> wrote:
>> > >
>> > >
>> > >
>> > > On Mon, Jan 8, 2018 at 2:15 PM, Smith, Barry F. <bsm...@mcs.anl.gov> 
>> > > wrote:
>> > >
>> > >
>> > > > On Jan 8, 2018, at 2:59 PM, Alexander Lindsay 
>> > > > <alexlindsay...@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>

Attachment: ex3.c
Description: ex3.c

Reply via email to