[petsc-users] Context behind SNESComputeFunction call
Fande Kong
fdkong.jd at gmail.com
Sat Jan 27 17:10:00 CST 2018
Thanks, Barry,
On Sat, Jan 27, 2018 at 3:10 PM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
>
> 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 like this design. I will try out at the moose side.
Fande,
>
>
> 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/033a079f/attachment.html>
More information about the petsc-users
mailing list