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>
