Actual source code: snesmfj.c


  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:     MatSNESMFGetSNES - returns the SNES associated with a matrix created with MatCreateSNESMF()

 56:     Not collective

 58:     Input Parameter:
 59: .   J - the matrix

 61:     Output Parameter:
 62: .   snes - the SNES object

 64:     Level: advanced

 66: .seealso: MatCreateSNESMF()
 67: @*/
 68: PetscErrorCode MatSNESMFGetSNES(Mat J,SNES *snes)
 69: {
 70:   MatMFFD        j;

 74:   MatShellGetContext(J,&j);
 75:   *snes = (SNES)j->ctx;
 76:   return(0);
 77: }

 79: /*
 80:    MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the
 81:     base from the SNES context

 83: */
 84: static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt)
 85: {
 87:   MatMFFD        j;
 88:   SNES           snes;
 89:   Vec            u,f;
 90:   DM             dm;
 91:   DMSNES         dms;

 94:   MatShellGetContext(J,&j);
 95:   snes = (SNES)j->ctx;
 96:   MatAssemblyEnd_MFFD(J,mt);

 98:   SNESGetSolution(snes,&u);
 99:   SNESGetDM(snes,&dm);
100:   DMGetDMSNES(dm,&dms);
101:   if ((j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction) && !dms->ops->computemffunction) {
102:     SNESGetFunction(snes,&f,NULL,NULL);
103:     MatMFFDSetBase_MFFD(J,u,f);
104:   } else {
105:     /* f value known by SNES is not correct for other differencing function */
106:     MatMFFDSetBase_MFFD(J,u,NULL);
107:   }
108:   return(0);
109: }

111: /*
112:    MatAssemblyEnd_SNESMF_UseBase - Calls MatAssemblyEnd_MFFD() and then sets the
113:     base from the SNES context. This version will cause the base to be used for differencing
114:     even if the func is not SNESComputeFunction. See: MatSNESMFUseBase()

116: */
117: static PetscErrorCode MatAssemblyEnd_SNESMF_UseBase(Mat J,MatAssemblyType mt)
118: {
120:   MatMFFD        j;
121:   SNES           snes;
122:   Vec            u,f;

125:   MatAssemblyEnd_MFFD(J,mt);
126:   MatShellGetContext(J,&j);
127:   snes = (SNES)j->ctx;
128:   SNESGetSolution(snes,&u);
129:   SNESGetFunction(snes,&f,NULL,NULL);
130:   MatMFFDSetBase_MFFD(J,u,f);
131:   return(0);
132: }

134: /*
135:     This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer
136:   uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF().
137: */
138: static PetscErrorCode  MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F)
139: {

143:   MatMFFDSetBase_MFFD(J,U,F);
144:   J->ops->assemblyend = MatAssemblyEnd_MFFD;
145:   return(0);
146: }

148: static PetscErrorCode  MatSNESMFSetReuseBase_SNESMF(Mat J,PetscBool use)
149: {
151:   if (use) {
152:     J->ops->assemblyend = MatAssemblyEnd_SNESMF_UseBase;
153:   } else {
154:     J->ops->assemblyend = MatAssemblyEnd_SNESMF;
155:   }
156:   return(0);
157: }

159: /*@
160:     MatSNESMFSetReuseBase - Causes the base vector to be used for differencing even if the function provided to SNESSetFunction() is not the
161:                        same as that provided to MatMFFDSetFunction().

163:     Logically Collective on Mat

165:     Input Parameters:
166: +   J - the MatMFFD matrix
167: -   use - if true always reuse the base vector instead of recomputing f(u) even if the function in the MatSNESMF is
168:           not SNESComputeFunction()

170:     Notes:
171:     Care must be taken when using this routine to insure that the function provided to MatMFFDSetFunction(), call it F_MF() is compatible with
172:            with that provided to SNESSetFunction(), call it F_SNES(). That is, (F_MF(u + h*d) - F_SNES(u))/h has to approximate the derivative

174:     Developer Notes:
175:     This was provided for the MOOSE team who desired to have a SNESSetFunction() function that could change configurations due
176:                      to contacts while the function provided to MatMFFDSetFunction() cannot. Except for the possibility of changing the configuration
177:                      both functions compute the same mathematical function so the differencing makes sense.

179:     Level: advanced

181: .seealso: MatCreateSNESMF(), MatSNESMFGetReuseBase()
182: @*/
183: PetscErrorCode  MatSNESMFSetReuseBase(Mat J,PetscBool use)
184: {

189:   PetscTryMethod(J,"MatSNESMFSetReuseBase_C",(Mat,PetscBool),(J,use));
190:   return(0);
191: }

193: static PetscErrorCode  MatSNESMFGetReuseBase_SNESMF(Mat J,PetscBool *use)
194: {
196:   if (J->ops->assemblyend == MatAssemblyEnd_SNESMF_UseBase) *use = PETSC_TRUE;
197:   else *use = PETSC_FALSE;
198:   return(0);
199: }

201: /*@
202:     MatSNESMFGetReuseBase - Determines if the base vector is to be used for differencing even if the function provided to SNESSetFunction() is not the
203:                        same as that provided to MatMFFDSetFunction().

205:     Logically Collective on Mat

207:     Input Parameter:
208: .   J - the MatMFFD matrix

210:     Output Parameter:
211: .   use - if true always reuse the base vector instead of recomputing f(u) even if the function in the MatSNESMF is
212:           not SNESComputeFunction()

214:     Notes:
215:     Care must be taken when using this routine to insure that the function provided to MatMFFDSetFunction(), call it F_MF() is compatible with
216:            with that provided to SNESSetFunction(), call it F_SNES(). That is, (F_MF(u + h*d) - F_SNES(u))/h has to approximate the derivative

218:     Developer Notes:
219:     This was provided for the MOOSE team who desired to have a SNESSetFunction() function that could change configurations due
220:                      to contacts while the function provided to MatMFFDSetFunction() cannot. Except for the possibility of changing the configuration
221:                      both functions compute the same mathematical function so the differencing makes sense.

223:     Level: advanced

225: .seealso: MatCreateSNESMF(), MatSNESMFSetReuseBase()
226: @*/
227: PetscErrorCode  MatSNESMFGetReuseBase(Mat J,PetscBool *use)
228: {

233:   PetscUseMethod(J,"MatSNESMFGetReuseBase_C",(Mat,PetscBool*),(J,use));
234:   return(0);
235: }

237: /*@
238:    MatCreateSNESMF - Creates a matrix-free matrix context for use with
239:    a SNES solver.  This matrix can be used as the Jacobian argument for
240:    the routine SNESSetJacobian(). See MatCreateMFFD() for details on how
241:    the finite difference computation is done.

243:    Collective on SNES

245:    Input Parameters:
246: .  snes - the SNES context

248:    Output Parameter:
249: .  J - the matrix-free matrix

251:    Level: advanced

253:    Notes:
254:      You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
255:      preconditioner matrix

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

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

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

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

272: .seealso: MatDestroy(), MatMFFDSetFunction(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin()
273:           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateMFFD(),
274:           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian(), MatSNESMFSetReuseBase(), MatSNESMFGetReuseBase()

276: @*/
277: PetscErrorCode  MatCreateSNESMF(SNES snes,Mat *J)
278: {
280:   PetscInt       n,N;
281:   MatMFFD        mf;

284:   if (snes->vec_func) {
285:     VecGetLocalSize(snes->vec_func,&n);
286:     VecGetSize(snes->vec_func,&N);
287:   } else if (snes->dm) {
288:     Vec tmp;
289:     DMGetGlobalVector(snes->dm,&tmp);
290:     VecGetLocalSize(tmp,&n);
291:     VecGetSize(tmp,&N);
292:     DMRestoreGlobalVector(snes->dm,&tmp);
293:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
294:   MatCreateMFFD(PetscObjectComm((PetscObject)snes),n,n,N,N,J);
295:   MatShellGetContext(*J,&mf);
296:   mf->ctx = snes;

298:   if (snes->npc && snes->npcside== PC_LEFT) {
299:     MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunctionDefaultNPC,snes);
300:   } else {
301:     DM     dm;
302:     DMSNES dms;

304:     SNESGetDM(snes,&dm);
305:     DMGetDMSNES(dm,&dms);
306:     MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))(dms->ops->computemffunction ? SNESComputeMFFunction : SNESComputeFunction),snes);
307:   }
308:   (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF;

310:   PetscObjectComposeFunction((PetscObject)*J,"MatMFFDSetBase_C",MatMFFDSetBase_SNESMF);
311:   PetscObjectComposeFunction((PetscObject)*J,"MatSNESMFSetReuseBase_C",MatSNESMFSetReuseBase_SNESMF);
312:   PetscObjectComposeFunction((PetscObject)*J,"MatSNESMFGetReuseBase_C",MatSNESMFGetReuseBase_SNESMF);
313:   return(0);
314: }