Actual source code: dmsnes.c

petsc-3.3-p7 2013-05-11
  1: #include <petsc-private/snesimpl.h>   /*I "petscsnes.h" I*/
  2: #include <petscdm.h>            /*I "petscdm.h" I*/

  6: /*@C
  7:   SNESDMComputeFunction - This is a universal function evaluation routine that
  8:   may be used with SNESSetFunction() as long as the user context has a DM
  9:   as its first record and the user has called DMSetLocalFunction().

 11:   Collective on SNES

 13:   Input Parameters:
 14: + snes - the SNES context
 15: . X - input vector
 16: . F - function vector
 17: - ptr - pointer to a structure that must have a DM as its first entry.
 18:         This ptr must have been passed into SNESDMComputeFunction() as the context.

 20:   Level: intermediate

 22: .seealso: DMSetLocalFunction(), DMSetLocalJacobian(), SNESSetFunction(), SNESSetJacobian()
 23: @*/
 24: PetscErrorCode SNESDMComputeFunction(SNES snes, Vec X, Vec F, void *ptr)
 25: {
 26:   DM               dm = *(DM*) ptr;
 27:   PetscErrorCode (*lf)(DM, Vec, Vec, void *);
 28:   Vec              localX, localF;
 29:   PetscInt         N, n;
 30:   PetscErrorCode   ierr;

 36:   if (!dm) SETERRQ(((PetscObject)snes)->comm, PETSC_ERR_ARG_WRONGSTATE, "Looks like you called SNESSetFromFuntion(snes,SNESDMComputeFunction,) without the DM context");

 39:   /* determine whether X = localX */
 40:   DMGetLocalVector(dm, &localX);
 41:   DMGetLocalVector(dm, &localF);
 42:   VecGetSize(X, &N);
 43:   VecGetSize(localX, &n);

 45:   if (n != N){ /* X != localX */
 46:     /* Scatter ghost points to local vector, using the 2-step process
 47:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
 48:     */
 49:     DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
 50:     DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
 51:   } else {
 52:     DMRestoreLocalVector(dm, &localX);
 53:     localX = X;
 54:   }
 55:   DMGetLocalFunction(dm, &lf);
 56:   (*lf)(dm, localX, localF, ptr);
 57:   if (n != N){
 58:     DMRestoreLocalVector(dm, &localX);
 59:   }
 60:   DMLocalToGlobalBegin(dm, localF, ADD_VALUES, F);
 61:   DMLocalToGlobalEnd(dm, localF, ADD_VALUES, F);
 62:   DMRestoreLocalVector(dm, &localF);
 63:   return(0);
 64: }

 68: /*
 69:   SNESDMComputeJacobian - This is a universal Jacobian evaluation routine that
 70:   may be used with SNESSetJacobian() as long as the user context has a DM
 71:   as its first record and the user has called DMSetLocalJacobian().

 73:   Collective on SNES

 75:   Input Parameters:
 76: + snes - the SNES context
 77: . X - input vector
 78: . J - Jacobian
 79: . B - Jacobian used in preconditioner (usally same as J)
 80: . flag - indicates if the matrix changed its structure
 81: - ptr - pointer to a structure that must have a DM as its first entry.
 82:         This ptr must have been passed into SNESDMComputeFunction() as the context.

 84:   Level: intermediate

 86: .seealso: DMSetLocalFunction(), DMSetLocalJacobian(), SNESSetFunction(), SNESSetJacobian()
 87: */
 88: PetscErrorCode SNESDMComputeJacobian(SNES snes, Vec X, Mat *J, Mat *B, MatStructure *flag, void *ptr)
 89: {
 90:   DM               dm = *(DM*) ptr;
 91:   PetscErrorCode (*lj)(DM, Vec, Mat, Mat, void *);
 92:   Vec              localX;
 93:   PetscErrorCode   ierr;

 96:   DMGetLocalVector(dm, &localX);
 97:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
 98:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
 99:   DMGetLocalJacobian(dm, &lj);
100:   (*lj)(dm, localX, *J, *B, ptr);
101:   DMRestoreLocalVector(dm, &localX);
102:   /* Assemble true Jacobian; if it is different */
103:   if (*J != *B) {
104:     MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY);
105:     MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY);
106:   }
107:   MatSetOption(*B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
108:   *flag = SAME_NONZERO_PATTERN;
109:   return(0);
110: }

114: static PetscErrorCode PetscContainerDestroy_SNESDM(void *ctx)
115: {
117:   SNESDM sdm = (SNESDM)ctx;

120:   if (sdm->destroy) {(*sdm->destroy)(sdm);}
121:   PetscFree(sdm);
122:   return(0);
123: }

127: /* Attaches the SNESDM to the coarse level.
128:  * Under what conditions should we copy versus duplicate?
129:  */
130: static PetscErrorCode DMCoarsenHook_SNESDM(DM dm,DM dmc,void *ctx)
131: {

135:   DMSNESCopyContext(dm,dmc);
136:   return(0);
137: }

141: /* This could restrict auxiliary information to the coarse level.
142:  */
143: static PetscErrorCode DMRestrictHook_SNESDM(DM dm,Mat Restrict,Vec rscale,Mat Inject,DM dmc,void *ctx)
144: {

147:   return(0);
148: }

152: /*@C
153:    DMSNESGetContext - get read-only private SNESDM context from a DM

155:    Not Collective

157:    Input Argument:
158: .  dm - DM to be used with SNES

160:    Output Argument:
161: .  snesdm - private SNESDM context

163:    Level: developer

165:    Notes:
166:    Use DMSNESGetContextWrite() if write access is needed. The DMSNESSetXXX API should be used wherever possible.

168: .seealso: DMSNESGetContextWrite()
169: @*/
170: PetscErrorCode DMSNESGetContext(DM dm,SNESDM *snesdm)
171: {
173:   PetscContainer container;
174:   SNESDM         sdm;

178:   PetscObjectQuery((PetscObject)dm,"SNESDM",(PetscObject*)&container);
179:   if (container) {
180:     PetscContainerGetPointer(container,(void**)snesdm);
181:   } else {
182:     PetscInfo(dm,"Creating new SNESDM\n");
183:     PetscContainerCreate(((PetscObject)dm)->comm,&container);
184:     PetscNewLog(dm,struct _n_SNESDM,&sdm);
185:     PetscContainerSetPointer(container,sdm);
186:     PetscContainerSetUserDestroy(container,PetscContainerDestroy_SNESDM);
187:     PetscObjectCompose((PetscObject)dm,"SNESDM",(PetscObject)container);
188:     DMCoarsenHookAdd(dm,DMCoarsenHook_SNESDM,DMRestrictHook_SNESDM,PETSC_NULL);
189:     PetscContainerGetPointer(container,(void**)snesdm);
190:     PetscContainerDestroy(&container);
191:   }
192:   return(0);
193: }

197: /*@C
198:    DMSNESGetContextWrite - get write access to private SNESDM context from a DM

200:    Not Collective

202:    Input Argument:
203: .  dm - DM to be used with SNES

205:    Output Argument:
206: .  snesdm - private SNESDM context

208:    Level: developer

210: .seealso: DMSNESGetContext()
211: @*/
212: PetscErrorCode DMSNESGetContextWrite(DM dm,SNESDM *snesdm)
213: {
215:   SNESDM         sdm;

219:   DMSNESGetContext(dm,&sdm);
220:   if (!sdm->originaldm) sdm->originaldm = dm;
221:   if (sdm->originaldm != dm) {  /* Copy on write */
222:     PetscContainer container;
223:     SNESDM         oldsdm = sdm;
224:     PetscInfo(dm,"Copying SNESDM due to write\n");
225:     PetscContainerCreate(((PetscObject)dm)->comm,&container);
226:     PetscNewLog(dm,struct _n_SNESDM,&sdm);
227:     PetscMemcpy(sdm,oldsdm,sizeof *sdm);
228:     PetscContainerSetPointer(container,sdm);
229:     PetscContainerSetUserDestroy(container,PetscContainerDestroy_SNESDM);
230:     PetscObjectCompose((PetscObject)dm,"SNESDM",(PetscObject)container);
231:     PetscContainerDestroy(&container);
232:   }
233:   *snesdm = sdm;
234:   return(0);
235: }

239: /*@C
240:    DMSNESCopyContext - copies a DM context to a new DM

242:    Logically Collective

244:    Input Arguments:
245: +  dmsrc - DM to obtain context from
246: -  dmdest - DM to add context to

248:    Level: developer

250:    Note:
251:    The context is copied by reference. This function does not ensure that a context exists.

253: .seealso: DMSNESGetContext(), SNESSetDM()
254: @*/
255: PetscErrorCode DMSNESCopyContext(DM dmsrc,DM dmdest)
256: {
258:   PetscContainer container;

263:   PetscObjectQuery((PetscObject)dmsrc,"SNESDM",(PetscObject*)&container);
264:   if (container) {
265:     PetscObjectCompose((PetscObject)dmdest,"SNESDM",(PetscObject)container);
266:     DMCoarsenHookAdd(dmdest,DMCoarsenHook_SNESDM,DMRestrictHook_SNESDM,PETSC_NULL);
267:   }
268:   return(0);
269: }

273: /*@C
274:    DMSNESSetFunction - set SNES residual evaluation function

276:    Not Collective

278:    Input Arguments:
279: +  dm - DM to be used with SNES
280: .  func - residual evaluation function, see SNESSetFunction() for calling sequence
281: -  ctx - context for residual evaluation

283:    Level: advanced

285:    Note:
286:    SNESSetFunction() is normally used, but it calls this function internally because the user context is actually
287:    associated with the DM.  This makes the interface consistent regardless of whether the user interacts with a DM or
288:    not. If DM took a more central role at some later date, this could become the primary method of setting the residual.

290: .seealso: DMSNESSetContext(), SNESSetFunction(), DMSNESSetJacobian()
291: @*/
292: PetscErrorCode DMSNESSetFunction(DM dm,PetscErrorCode (*func)(SNES,Vec,Vec,void*),void *ctx)
293: {
295:   SNESDM sdm;

299:   DMSNESGetContextWrite(dm,&sdm);
300:   if (func) sdm->computefunction = func;
301:   if (ctx)  sdm->functionctx = ctx;
302:   return(0);
303: }

307: /*@C
308:    DMSNESGetFunction - get SNES residual evaluation function

310:    Not Collective

312:    Input Argument:
313: .  dm - DM to be used with SNES

315:    Output Arguments:
316: +  func - residual evaluation function, see SNESSetFunction() for calling sequence
317: -  ctx - context for residual evaluation

319:    Level: advanced

321:    Note:
322:    SNESGetFunction() is normally used, but it calls this function internally because the user context is actually
323:    associated with the DM.

325: .seealso: DMSNESSetContext(), DMSNESSetFunction(), SNESSetFunction()
326: @*/
327: PetscErrorCode DMSNESGetFunction(DM dm,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx)
328: {
330:   SNESDM sdm;

334:   DMSNESGetContext(dm,&sdm);
335:   if (func) *func = sdm->computefunction;
336:   if (ctx)  *ctx = sdm->functionctx;
337:   return(0);
338: }

342: /*@C
343:    DMSNESSetGS - set SNES Gauss-Seidel relaxation function

345:    Not Collective

347:    Input Argument:
348: +  dm - DM to be used with SNES
349: .  func - relaxation function, see SNESSetGS() for calling sequence
350: -  ctx - context for residual evaluation

352:    Level: advanced

354:    Note:
355:    SNESSetGS() is normally used, but it calls this function internally because the user context is actually
356:    associated with the DM.  This makes the interface consistent regardless of whether the user interacts with a DM or
357:    not. If DM took a more central role at some later date, this could become the primary method of setting the residual.

359: .seealso: DMSNESSetContext(), SNESSetFunction(), DMSNESSetJacobian(), DMSNESSetFunction()
360: @*/
361: PetscErrorCode DMSNESSetGS(DM dm,PetscErrorCode (*func)(SNES,Vec,Vec,void*),void *ctx)
362: {
364:   SNESDM sdm;

368:   DMSNESGetContextWrite(dm,&sdm);
369:   if (func) sdm->computegs = func;
370:   if (ctx)  sdm->gsctx = ctx;
371:   return(0);
372: }

376: /*@C
377:    DMSNESGetGS - get SNES Gauss-Seidel relaxation function

379:    Not Collective

381:    Input Argument:
382: .  dm - DM to be used with SNES

384:    Output Arguments:
385: +  func - relaxation function, see SNESSetGS() for calling sequence
386: -  ctx - context for residual evaluation

388:    Level: advanced

390:    Note:
391:    SNESGetGS() is normally used, but it calls this function internally because the user context is actually
392:    associated with the DM.  This makes the interface consistent regardless of whether the user interacts with a DM or
393:    not. If DM took a more central role at some later date, this could become the primary method of setting the residual.

395: .seealso: DMSNESSetContext(), SNESGetGS(), DMSNESGetJacobian(), DMSNESGetFunction()
396: @*/
397: PetscErrorCode DMSNESGetGS(DM dm,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx)
398: {
400:   SNESDM sdm;

404:   DMSNESGetContext(dm,&sdm);
405:   if (func) *func = sdm->computegs;
406:   if (ctx)  *ctx = sdm->gsctx;
407:   return(0);
408: }

412: /*@C
413:    DMSNESSetJacobian - set SNES Jacobian evaluation function

415:    Not Collective

417:    Input Argument:
418: +  dm - DM to be used with SNES
419: .  func - Jacobian evaluation function, see SNESSetJacobian() for calling sequence
420: -  ctx - context for residual evaluation

422:    Level: advanced

424:    Note:
425:    SNESSetJacobian() is normally used, but it calls this function internally because the user context is actually
426:    associated with the DM.  This makes the interface consistent regardless of whether the user interacts with a DM or
427:    not. If DM took a more central role at some later date, this could become the primary method of setting the Jacobian.

429: .seealso: DMSNESSetContext(), SNESSetFunction(), DMSNESGetJacobian(), SNESSetJacobian()
430: @*/
431: PetscErrorCode DMSNESSetJacobian(DM dm,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
432: {
434:   SNESDM sdm;

438:   DMSNESGetContextWrite(dm,&sdm);
439:   if (func) sdm->computejacobian = func;
440:   if (ctx)  sdm->jacobianctx = ctx;
441:   return(0);
442: }

446: /*@C
447:    DMSNESGetJacobian - get SNES Jacobian evaluation function

449:    Not Collective

451:    Input Argument:
452: .  dm - DM to be used with SNES

454:    Output Arguments:
455: +  func - Jacobian evaluation function, see SNESSetJacobian() for calling sequence
456: -  ctx - context for residual evaluation

458:    Level: advanced

460:    Note:
461:    SNESGetJacobian() is normally used, but it calls this function internally because the user context is actually
462:    associated with the DM.  This makes the interface consistent regardless of whether the user interacts with a DM or
463:    not. If DM took a more central role at some later date, this could become the primary method of setting the Jacobian.

465: .seealso: DMSNESSetContext(), SNESSetFunction(), DMSNESSetJacobian()
466: @*/
467: PetscErrorCode DMSNESGetJacobian(DM dm,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
468: {
470:   SNESDM sdm;

474:   DMSNESGetContext(dm,&sdm);
475:   if (func) *func = sdm->computejacobian;
476:   if (ctx)  *ctx = sdm->jacobianctx;
477:   return(0);
478: }

482: /*@C
483:    DMSNESSetPicard - set SNES Picard iteration matrix and RHS evaluation functions.

485:    Not Collective

487:    Input Argument:
488: +  dm - DM to be used with SNES
489: .  func - RHS evaluation function, see SNESSetFunction() for calling sequence
490: .  pjac - Picard matrix evaluation function, see SNESSetJacobian() for calling sequence
491: -  ctx - context for residual evaluation

493:    Level: advanced

495: .seealso: SNESSetPicard(), DMSNESSetFunction(), DMSNESSetJacobian()
496: @*/
497: PetscErrorCode DMSNESSetPicard(DM dm,PetscErrorCode (*pfunc)(SNES,Vec,Vec,void*),PetscErrorCode (*pjac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx)
498: {
500:   SNESDM sdm;

504:   DMSNESGetContext(dm,&sdm);
505:   if (pfunc) sdm->computepfunction = pfunc;
506:   if (pjac)  sdm->computepjacobian = pjac;
507:   if (ctx)   sdm->pctx             = ctx;
508:   return(0);
509: }


514: /*@C
515:    DMSNESGetPicard - get SNES Picard iteration evaluation functions

517:    Not Collective

519:    Input Argument:
520: .  dm - DM to be used with SNES

522:    Output Arguments:
523: +  pfunc - Jacobian evaluation function, see SNESSetJacobian() for calling sequence
524: .  pjac  - RHS evaluation function, see SNESSetFunction() for calling sequence
525: -  ctx - context for residual evaluation

527:    Level: advanced

529: .seealso: DMSNESSetContext(), SNESSetFunction(), DMSNESSetJacobian()
530: @*/
531: PetscErrorCode DMSNESGetPicard(DM dm,PetscErrorCode (**pfunc)(SNES,Vec,Vec,void*),PetscErrorCode (**pjac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx)
532: {
534:   SNESDM sdm;

538:   DMSNESGetContext(dm,&sdm);
539:   if (pfunc) *pfunc = sdm->computepfunction;
540:   if (pjac) *pjac   = sdm->computepjacobian;
541:   if (ctx)  *ctx    = sdm->pctx;
542:   return(0);
543: }


548: static PetscErrorCode SNESDefaultComputeFunction_DMLegacy(SNES snes,Vec X,Vec F,void *ctx)
549: {
551:   DM             dm;

554:   SNESGetDM(snes,&dm);
555:   DMComputeFunction(dm,X,F);
556:   return(0);
557: }

561: static PetscErrorCode SNESDefaultComputeJacobian_DMLegacy(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
562: {
564:   DM             dm;

567:   SNESGetDM(snes,&dm);
568:   DMComputeJacobian(dm,X,*A,*B,mstr);
569:   return(0);
570: }

574: /* Sets up calling of legacy DM routines */
575: PetscErrorCode DMSNESSetUpLegacy(DM dm)
576: {
578:   SNESDM         sdm;

581:   DMSNESGetContext(dm,&sdm);
582:   if (!sdm->computefunction) {DMSNESSetFunction(dm,SNESDefaultComputeFunction_DMLegacy,PETSC_NULL);}
583:   DMSNESGetContext(dm,&sdm);
584:   if (!sdm->computejacobian) {DMSNESSetJacobian(dm,SNESDefaultComputeJacobian_DMLegacy,PETSC_NULL);}
585:   return(0);
586: }