[petsc-users] Context behind SNESComputeFunction call
Smith, Barry F.
bsmith at mcs.anl.gov
Sat Jan 27 16:10:06 CST 2018
Fande,
I have done something similar in the branch https://bitbucket.org/petsc/petsc/pull-requests/845/added-matsnesmfsetreusebase-at-request-of/diff attached is your test case.
I don't think a command line option makes sense for this functionality.
Barry
I didn't do it with a SNESSetReuseBaseMFFD() because I didn't want the functionality of the MatSNESMFCreate() function to "bleed over" into the SNES object. That is, it is not the business of the SNES object to even know about details of specific matrix types it might use.
> On Jan 27, 2018, at 11:18 AM, Fande Kong <fdkong.jd at gmail.com> wrote:
>
> Hi Barry,
>
> Thanks so much!
>
> Does this requires us to attach a shell function "MatAssemblyEnd_SNESMF"? We need to maintain "MatAssemblyEnd_SNESMF" in the moose side?
>
>
>
> Could we just add a flag into the original routine "MatAssemblyEnd_SNESMF" in PETSc?
>
>
> static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt)
> {
> PetscErrorCode ierr;
> MatMFFD j = (MatMFFD)J->data;
> SNES snes = (SNES)j->ctx;
> PetscBool refuse;
> Vec u,f;
>
> PetscFunctionBegin;
> ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr);
>
> ierr = SNESGetReuseBaseMFFD(SNES snes, & refuse);CHKERRQ(ierr);
>
> ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr);
> if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction || reuse) {
> 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);
> }
> PetscFunctionReturn(0);
> }
>
>
>
> We may need new APIs SNESGetReuseBaseMFFD and SNESSetReuseBaseMFFD, and a command line option "-snes_reuse_base_mffd". "reuse" is false by default, and it is consistent with current setting.
>
>
> Please let me know what do you think. If you agree, I can make a PR on this.
>
> Fande,
>
>
> On Fri, Jan 26, 2018 at 10:47 PM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
>
> 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/31638c09/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex3.c
Type: application/octet-stream
Size: 13541 bytes
Desc: ex3.c
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180127/31638c09/attachment-0001.obj>
More information about the petsc-users
mailing list