Permalink
Browse files

fix MATSNESMF so that MatMFFDSetFunction() works correctly with it

(cherry picked from commit 5eb111a)

Conflicts:
	src/snes/mf/snesmfj.c
  • Loading branch information...
1 parent 4ded9de commit 18d1b208785a8fad2f6fed974bf524d1bc08f966 @BarrySmith BarrySmith committed with jedbrown Feb 13, 2014
Showing with 46 additions and 8 deletions.
  1. +1 −0 src/mat/impls/mffd/mffdimpl.h
  2. +20 −0 src/snes/examples/tests/ex7.c
  3. +25 −8 src/snes/mf/snesmfj.c
@@ -52,6 +52,7 @@ struct _p_MatMFFD { /* context for default matrix-free SNES */
PetscScalar vscale,vshift; /* diagonal scale and shift by scalars */
Vec dlscale,drscale,dshift; /* diagonal scale and shift by vectors */
+ void *ctx; /* this is used by MatCreateSNESMF() to store the SNES object */
};
PETSC_EXTERN PetscFunctionList MatMFFDList;
@@ -6,6 +6,7 @@ static char help[] = "Solves u`` + u^{2} = f with Newton-like methods. Using\n\
extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
+extern PetscErrorCode OtherFunctionForDifferencing(void*,Vec,Vec);
extern PetscErrorCode FormInitialGuess(SNES,Vec);
extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
@@ -66,11 +67,15 @@ int main(int argc,char **argv)
/* Set various routines and options */
ierr = SNESSetFunction(snes,r,FormFunction,F);CHKERRQ(ierr);
if (user.variant) {
+ /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */
ierr = MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J);CHKERRQ(ierr);
ierr = MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes);CHKERRQ(ierr);
} else {
/* create matrix free matrix for Jacobian */
ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr);
+ /* demonstrates differencing a different function than FormFunction() to apply a matrix operator */
+ /* note we use the same context for this function as FormFunction, the F vector */
+ ierr = MatMFFDSetFunction(J,OtherFunctionForDifferencing,F);CHKERRQ(ierr);
}
/* Set various routines and options */
@@ -115,6 +120,21 @@ PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
ierr = VecRestoreArray((Vec)dummy,&FF);CHKERRQ(ierr);
return 0;
}
+
+/*
+
+ Example function that when differenced produces the same matrix free Jacobian as FormFunction()
+ this is provided to show how a user can provide a different function
+*/
+PetscErrorCode OtherFunctionForDifferencing(void *dummy,Vec x,Vec f)
+{
+ PetscErrorCode ierr;
+
+ ierr = FormFunction(NULL,x,f,dummy);CHKERRQ(ierr);
+ ierr = VecShift(f,1.0);CHKERRQ(ierr);
+ return 0;
+}
+
/* -------------------- Form initial approximation ----------------- */
#undef __FUNCT__
View
@@ -63,15 +63,20 @@ PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt)
{
PetscErrorCode ierr;
MatMFFD j = (MatMFFD)J->data;
- SNES snes = (SNES)j->funcctx;
+ SNES snes = (SNES)j->ctx;
Vec u,f;
PetscFunctionBegin;
ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr);
ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr);
- ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr);
- ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr);
+ 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);
+ }
PetscFunctionReturn(0);
}
@@ -110,15 +115,23 @@ PetscErrorCode MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F)
Level: advanced
- Warning:
- If MatMFFDSetBase() is ever called on jac then this routine will NO longer get
- the x from the SNES object and MatMFFDSetBase() must from that point on be used to
- change the base vector x.
- Notes: The difference between this routine and MatCreateMFFD() is that this matrix
+ Notes:
+ You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
+ preconditioner matrix
+
+ If you wish to provide a different function to do differencing on to compute the matrix-free operator than
+ that provided to SNESSetFunction() then call MatMFFDSetFunction() with your function after this call.
+
+ The difference between this routine and MatCreateMFFD() is that this matrix
automatically gets the current base vector from the SNES object and not from an
explicit call to MatMFFDSetBase().
+ Warning:
+ If MatMFFDSetBase() is ever called on jac then this routine will NO longer get
+ the x from the SNES object and MatMFFDSetBase() must from that point on be used to
+ change the base vector x.
+
.seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin()
MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateMFFD(),
MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
@@ -128,6 +141,7 @@ PetscErrorCode MatCreateSNESMF(SNES snes,Mat *J)
{
PetscErrorCode ierr;
PetscInt n,N;
+ MatMFFD mf;
PetscFunctionBegin;
if (snes->vec_func) {
@@ -141,6 +155,9 @@ PetscErrorCode MatCreateSNESMF(SNES snes,Mat *J)
ierr = DMRestoreGlobalVector(snes->dm,&tmp);CHKERRQ(ierr);
} else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
ierr = MatCreateMFFD(PetscObjectComm((PetscObject)snes),n,n,N,N,J);CHKERRQ(ierr);
+ mf = (MatMFFD)(*J)->data;
+ mf->ctx = snes;
+
ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);CHKERRQ(ierr);
(*J)->ops->assemblyend = MatAssemblyEnd_SNESMF;

0 comments on commit 18d1b20

Please sign in to comment.