Actual source code: snesmfj.c

petsc-master 2017-07-25
Report Typos and Errors

  2:  #include <petsc/private/snesimpl.h>
  3:  #include <petscdm.h>
  4:  #include <../src/mat/impls/mffd/mffdimpl.h>
  5:  #include <petsc/private/matimpl.h>

  7: /*@C
  8:    MatMFFDComputeJacobian - Tells the matrix-free Jacobian object the new location at which
  9:        Jacobian matrix vector products will be computed at, i.e. J(x) * a. The x is obtained
 10:        from the SNES object (using SNESGetSolution()).

 12:    Logically Collective on SNES

 14:    Input Parameters:
 15: +   snes - the nonlinear solver context
 16: .   x - the point at which the Jacobian vector products will be performed
 17: .   jac - the matrix-free Jacobian object
 18: .   B - either the same as jac or another matrix type (ignored)
 19: .   flag - not relevent for matrix-free form
 20: -   dummy - the user context (ignored)

 22:    Level: developer

 24:    Warning:
 25:       If MatMFFDSetBase() is ever called on jac then this routine will NO longer get
 26:     the x from the SNES object and MatMFFDSetBase() must from that point on be used to
 27:     change the base vector x.

 29:    Notes:
 30:      This can be passed into SNESSetJacobian() as the Jacobian evaluation function argument
 31:      when using a completely matrix-free solver,
 32:      that is the B matrix is also the same matrix operator. This is used when you select
 33:      -snes_mf but rarely used directly by users. (All this routine does is call MatAssemblyBegin/End() on
 34:      the Mat jac.)

 36: .seealso: MatMFFDGetH(), MatCreateSNESMF(), MatCreateMFFD(), MATMFFD,
 37:           MatMFFDSetHHistory(), MatMFFDSetFunctionError(), MatCreateMFFD(), SNESSetJacobian()

 39: @*/
 40: PetscErrorCode  MatMFFDComputeJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
 41: {

 45:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
 46:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
 47:   return(0);
 48: }

 50: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat,MatAssemblyType);
 51: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat,Vec,Vec);

 53: /*
 54:    MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the
 55:     base from the SNES context

 57: */
 58: static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt)
 59: {
 61:   MatMFFD        j    = (MatMFFD)J->data;
 62:   SNES           snes = (SNES)j->ctx;
 63:   Vec            u,f;

 66:   MatAssemblyEnd_MFFD(J,mt);

 68:   SNESGetSolution(snes,&u);
 69:   if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction) {
 70:     SNESGetFunction(snes,&f,NULL,NULL);
 71:     MatMFFDSetBase_MFFD(J,u,f);
 72:   } else {
 73:     /* f value known by SNES is not correct for other differencing function */
 74:     MatMFFDSetBase_MFFD(J,u,NULL);
 75:   }
 76:   return(0);
 77: }

 79: /*
 80:     This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer
 81:   uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF().
 82: */
 83: static PetscErrorCode  MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F)
 84: {

 88:   MatMFFDSetBase_MFFD(J,U,F);

 90:   J->ops->assemblyend = MatAssemblyEnd_MFFD;
 91:   return(0);
 92: }

 94: /*@
 95:    MatCreateSNESMF - Creates a matrix-free matrix context for use with
 96:    a SNES solver.  This matrix can be used as the Jacobian argument for
 97:    the routine SNESSetJacobian(). See MatCreateMFFD() for details on how
 98:    the finite difference computation is done.

100:    Collective on SNES and Vec

102:    Input Parameters:
103: .  snes - the SNES context

105:    Output Parameter:
106: .  J - the matrix-free matrix

108:    Level: advanced


111:    Notes:
112:      You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
113:      preconditioner matrix

115:      If you wish to provide a different function to do differencing on to compute the matrix-free operator than
116:      that provided to SNESSetFunction() then call MatMFFDSetFunction() with your function after this call.

118:      The difference between this routine and MatCreateMFFD() is that this matrix
119:      automatically gets the current base vector from the SNES object and not from an
120:      explicit call to MatMFFDSetBase().

122:    Warning:
123:      If MatMFFDSetBase() is ever called on jac then this routine will NO longer get
124:      the x from the SNES object and MatMFFDSetBase() must from that point on be used to
125:      change the base vector x.

127:    Warning:
128:      Using a different function for the differencing will not work if you are using non-linear left preconditioning.


131: .seealso: MatDestroy(), MatMFFDSetFunction(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin()
132:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateMFFD(),
133:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()

135: @*/
136: PetscErrorCode  MatCreateSNESMF(SNES snes,Mat *J)
137: {
139:   PetscInt       n,N;
140:   MatMFFD        mf;

143:   if (snes->vec_func) {
144:     VecGetLocalSize(snes->vec_func,&n);
145:     VecGetSize(snes->vec_func,&N);
146:   } else if (snes->dm) {
147:     Vec tmp;
148:     DMGetGlobalVector(snes->dm,&tmp);
149:     VecGetLocalSize(tmp,&n);
150:     VecGetSize(tmp,&N);
151:     DMRestoreGlobalVector(snes->dm,&tmp);
152:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
153:   MatCreateMFFD(PetscObjectComm((PetscObject)snes),n,n,N,N,J);
154:   mf      = (MatMFFD)(*J)->data;
155:   mf->ctx = snes;

157:   if (snes->npc && snes->npcside== PC_LEFT) {
158:     MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunctionDefaultNPC,snes);
159:   } else {
160:     MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);
161:   }

163:   (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF;

165:   PetscObjectComposeFunction((PetscObject)*J,"MatMFFDSetBase_C",MatMFFDSetBase_SNESMF);
166:   return(0);
167: }