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 = 3Norm 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 = 3Norm 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
> > >
> > >
>
>
static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
This example employs a user-defined monitoring routine.\n\n";
/*T
Concepts: SNES^basic uniprocessor example
Concepts: SNES^setting a user-defined monitoring routine
Processors: 1
T*/
/*
Include "petscdraw.h" so that we can use PETSc drawing routines.
Include "petscsnes.h" so that we can use SNES solvers. Note that this
file automatically includes:
petscsys.h - base PETSc routines petscvec.h - vectors
petscmat.h - matrices
petscis.h - index sets petscksp.h - Krylov subspace methods
petscviewer.h - viewers petscpc.h - preconditioners
petscksp.h - linear solvers
*/
#include <petscsnes.h>
/*
User-defined routines
*/
extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
extern PetscErrorCode FormInitialGuess(Vec);
extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
extern PetscErrorCode FormFunction_SNES(SNES snes,Vec x, Vec r,void* ctx){
FormFunction(snes,x,r,ctx);
return 0;
}
extern PetscErrorCode FormFunction_MFFD(void*snes_ctx,Vec x, Vec r){
SNES snes = (SNES) snes_ctx;
void *ctx;
SNESGetFunction(snes,NULL,NULL,&ctx);
FormFunction(snes,x,r,ctx);
return 0;
}
/*
User-defined context for monitoring
*/
typedef struct {
PetscViewer viewer;
} MonitorCtx;
int main(int argc,char **argv)
{
SNES snes; /* SNES context */
Vec x,r,F,U; /* vectors */
Mat J; /* Jacobian matrix */
MonitorCtx monP; /* monitoring context */
PetscErrorCode ierr;
PetscInt its,n = 5,i,maxit,maxf;
PetscMPIInt size;
PetscScalar h,xp,v,none = -1.0;
PetscReal abstol,rtol,stol,norm;
ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
h = 1.0/(n-1);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create nonlinear solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create vector data structures; set function evaluation routine
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Note that we form 1 vector from scratch and then duplicate as needed.
*/
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
ierr = VecDuplicate(x,&F);CHKERRQ(ierr);
ierr = VecDuplicate(x,&U);CHKERRQ(ierr);
/*
Set function evaluation routine and vector
*/
ierr = SNESSetFunction(snes,r,FormFunction_SNES,(void*)F);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create matrix data structure; set Jacobian evaluation routine
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(J);CHKERRQ(ierr);
ierr = MatSeqAIJSetPreallocation(J,3,NULL);CHKERRQ(ierr);
/*
Set Jacobian matrix data structure and default Jacobian evaluation
routine. User can override with:
-snes_fd : default finite differencing approximation of Jacobian
-snes_mf : matrix-free Newton-Krylov method with no preconditioning
(unless user explicitly sets preconditioner)
-snes_mf_operator : form preconditioning matrix as set by the user,
but use matrix-free approx for Jacobian-vector
products within Newton-Krylov method
*/
ierr = SNESSetJacobian(snes,NULL,J,FormJacobian,NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Customize nonlinear solver; set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Set an optional user-defined monitoring routine
*/
ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);CHKERRQ(ierr);
/*ierr = SNESMonitorSet(snes,Monitor,&monP,0);CHKERRQ(ierr);*/
/*
Set names for some vectors to facilitate monitoring (optional)
*/
ierr = PetscObjectSetName((PetscObject)x,"Approximate Solution");CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)U,"Exact Solution");CHKERRQ(ierr);
/*
Set SNES/KSP/KSP/PC runtime options, e.g.,
-snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
*/
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
ierr = SNESSetUp(snes);CHKERRQ(ierr);
Mat Jacobian;
ierr = SNESGetJacobian(snes,&Jacobian,NULL,NULL,NULL);CHKERRQ(ierr);
/* ierr = MatMFFDSetFunction(Jacobian,FormFunction_MFFD,(void*)snes);CHKERRQ(ierr); */
/*
Print parameters used for convergence testing (optional) ... just
to demonstrate this routine; this information is also printed with
the option -snes_view
*/
ierr = SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Initialize application:
Store right-hand-side of PDE and exact solution
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
xp = 0.0;
for (i=0; i<n; i++) {
v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
ierr = VecSetValues(F,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
v = xp*xp*xp;
ierr = VecSetValues(U,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
xp += h;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Evaluate initial guess; then solve nonlinear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Note: The user should initialize the vector, x, with the initial guess
for the nonlinear solver prior to calling SNESSolve(). In particular,
to employ an initial guess of zero, the user should explicitly set
this vector to zero by calling VecSet().
*/
ierr = FormInitialGuess(x);CHKERRQ(ierr);
ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n\n",its);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Check solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Check the error
*/
ierr = VecAXPY(x,none,U);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr);
/*
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
*/
ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = VecDestroy(&U);CHKERRQ(ierr); ierr = VecDestroy(&F);CHKERRQ(ierr);
ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&monP.viewer);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
/* ------------------------------------------------------------------- */
/*
FormInitialGuess - Computes initial guess.
Input/Output Parameter:
. x - the solution vector
*/
PetscErrorCode FormInitialGuess(Vec x)
{
PetscErrorCode ierr;
PetscScalar pfive = .50;
ierr = VecSet(x,pfive);CHKERRQ(ierr);
return 0;
}
/* ------------------------------------------------------------------- */
/*
FormFunction - Evaluates nonlinear function, F(x).
Input Parameters:
. snes - the SNES context
. x - input vector
. ctx - optional user-defined context, as set by SNESSetFunction()
Output Parameter:
. f - function vector
Note:
The user-defined context can contain any application-specific data
needed for the function evaluation (such as various parameters, work
vectors, and grid information). In this program the context is just
a vector containing the right-hand-side of the discretized PDE.
*/
PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
{
Vec g = (Vec)ctx;
const PetscScalar *xx,*gg;
PetscScalar *ff,d;
PetscErrorCode ierr;
PetscInt i,n;
ierr = PetscPrintf(PETSC_COMM_WORLD," FormFunction is called \n");CHKERRQ(ierr);
/*
Get pointers to vector data.
- For default PETSc vectors, VecGetArray() returns a pointer to
the data array. Otherwise, the routine is implementation dependent.
- You MUST call VecRestoreArray() when you no longer need access to
the array.
*/
ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
ierr = VecGetArray(f,&ff);CHKERRQ(ierr);
ierr = VecGetArrayRead(g,&gg);CHKERRQ(ierr);
/*
Compute function
*/
ierr = VecGetSize(x,&n);CHKERRQ(ierr);
d = (PetscReal)(n - 1); d = d*d;
ff[0] = xx[0];
for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
ff[n-1] = xx[n-1] - 1.0;
/*
Restore vectors
*/
ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(g,&gg);CHKERRQ(ierr);
return 0;
}
/* ------------------------------------------------------------------- */
/*
FormJacobian - Evaluates Jacobian matrix.
Input Parameters:
. snes - the SNES context
. x - input vector
. dummy - optional user-defined context (not used here)
Output Parameters:
. jac - Jacobian matrix
. B - optionally different preconditioning matrix
*/
PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
{
const PetscScalar *xx;
PetscScalar A[3],d;
PetscErrorCode ierr;
PetscInt i,n,j[3];
/*
Get pointer to vector data
*/
ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
/*
Compute Jacobian entries and insert into matrix.
- Note that in this case we set all elements for a particular
row at once.
*/
ierr = VecGetSize(x,&n);CHKERRQ(ierr);
d = (PetscReal)(n - 1); d = d*d;
/*
Interior grid points
*/
for (i=1; i<n-1; i++) {
j[0] = i - 1; j[1] = i; j[2] = i + 1;
A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
ierr = MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);CHKERRQ(ierr);
}
/*
Boundary points
*/
i = 0; A[0] = 1.0;
ierr = MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);CHKERRQ(ierr);
i = n-1; A[0] = 1.0;
ierr = MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);CHKERRQ(ierr);
/*
Restore vector
*/
ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
/*
Assemble matrix
*/
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (jac != B) {
ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
return 0;
}
/* ------------------------------------------------------------------- */
/*
Monitor - User-defined monitoring routine that views the
current iterate with an x-window plot.
Input Parameters:
snes - the SNES context
its - iteration number
norm - 2-norm function value (may be estimated)
ctx - optional user-defined context for private data for the
monitor routine, as set by SNESMonitorSet()
Note:
See the manpage for PetscViewerDrawOpen() for useful runtime options,
such as -nox to deactivate all x-window output.
*/
PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
{
PetscErrorCode ierr;
MonitorCtx *monP = (MonitorCtx*) ctx;
Vec x;
ierr = PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm);CHKERRQ(ierr);
ierr = SNESGetSolution(snes,&x);CHKERRQ(ierr);
/*ierr = VecView(x,monP->viewer);CHKERRQ(ierr);*/
return 0;
}
/*TEST
test:
args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
test:
suffix: 2
args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
requires: !single
test:
suffix: 3
args: -nox -snes_monitor_cancel -snes_monitor_short -malloc no -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
TEST*/
