Thanks, Barry, On Sat, Jan 27, 2018 at 3:10 PM, Smith, Barry F. <[email protected]> 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 <[email protected]> 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. <[email protected]> > 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 < > [email protected]> 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 <[email protected]> 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. <[email protected]> > 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 <[email protected]> > wrote: > > >> > > > >> > > > >> > > > >> > On Fri, Jan 26, 2018 at 3:10 PM, Smith, Barry F. < > [email protected]> wrote: > > >> > > > >> > > > >> > > On Jan 26, 2018, at 2:15 PM, Kong, Fande <[email protected]> > wrote: > > >> > > > > >> > > > > >> > > > > >> > > On Mon, Jan 8, 2018 at 2:15 PM, Smith, Barry F. < > [email protected]> wrote: > > >> > > > > >> > > > > >> > > > On Jan 8, 2018, at 2:59 PM, Alexander Lindsay < > [email protected]> 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> > > > > > >
