Actual source code: dmdasnes.c

petsc-3.6.3 2015-12-03
Report Typos and Errors
  1: #include <petscdmda.h>          /*I "petscdmda.h" I*/
  2: #include <petsc/private/dmimpl.h>
  3: #include <petsc/private/snesimpl.h>   /*I "petscsnes.h" I*/

  5: /* This structure holds the user-provided DMDA callbacks */
  6: typedef struct {
  7:   PetscErrorCode (*residuallocal)(DMDALocalInfo*,void*,void*,void*);
  8:   PetscErrorCode (*jacobianlocal)(DMDALocalInfo*,void*,Mat,Mat,void*);
  9:   PetscErrorCode (*objectivelocal)(DMDALocalInfo*,void*,PetscReal*,void*);
 10:   void       *residuallocalctx;
 11:   void       *jacobianlocalctx;
 12:   void       *objectivelocalctx;
 13:   InsertMode residuallocalimode;

 15:   /*   For Picard iteration defined locally */
 16:   PetscErrorCode (*rhsplocal)(DMDALocalInfo*,void*,void*,void*);
 17:   PetscErrorCode (*jacobianplocal)(DMDALocalInfo*,void*,Mat,Mat,void*);
 18:   void *picardlocalctx;
 19: } DMSNES_DA;

 23: static PetscErrorCode DMSNESDestroy_DMDA(DMSNES sdm)
 24: {

 28:   PetscFree(sdm->data);
 29:   return(0);
 30: }

 34: static PetscErrorCode DMSNESDuplicate_DMDA(DMSNES oldsdm,DMSNES sdm)
 35: {

 39:   PetscNewLog(sdm,(DMSNES_DA**)&sdm->data);
 40:   if (oldsdm->data) {
 41:     PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMSNES_DA));
 42:   }
 43:   return(0);
 44: }


 49: static PetscErrorCode DMDASNESGetContext(DM dm,DMSNES sdm,DMSNES_DA  **dmdasnes)
 50: {

 54:   *dmdasnes = NULL;
 55:   if (!sdm->data) {
 56:     PetscNewLog(dm,(DMSNES_DA**)&sdm->data);
 57:     sdm->ops->destroy   = DMSNESDestroy_DMDA;
 58:     sdm->ops->duplicate = DMSNESDuplicate_DMDA;
 59:   }
 60:   *dmdasnes = (DMSNES_DA*)sdm->data;
 61:   return(0);
 62: }

 66: static PetscErrorCode SNESComputeFunction_DMDA(SNES snes,Vec X,Vec F,void *ctx)
 67: {
 69:   DM             dm;
 70:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
 71:   DMDALocalInfo  info;
 72:   Vec            Xloc;
 73:   void           *x,*f;

 79:   if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
 80:   SNESGetDM(snes,&dm);
 81:   DMGetLocalVector(dm,&Xloc);
 82:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
 83:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
 84:   DMDAGetLocalInfo(dm,&info);
 85:   DMDAVecGetArray(dm,Xloc,&x);
 86:   switch (dmdasnes->residuallocalimode) {
 87:   case INSERT_VALUES: {
 88:     DMDAVecGetArray(dm,F,&f);
 89:     PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);
 90:     CHKMEMQ;
 91:     (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
 92:     CHKMEMQ;
 93:     PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);
 94:     DMDAVecRestoreArray(dm,F,&f);
 95:   } break;
 96:   case ADD_VALUES: {
 97:     Vec Floc;
 98:     DMGetLocalVector(dm,&Floc);
 99:     VecZeroEntries(Floc);
100:     DMDAVecGetArray(dm,Floc,&f);
101:     PetscLogEventBegin(SNES_FunctionEval,snes,X,F,0);
102:     CHKMEMQ;
103:     (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
104:     CHKMEMQ;
105:     PetscLogEventEnd(SNES_FunctionEval,snes,X,F,0);
106:     DMDAVecRestoreArray(dm,Floc,&f);
107:     VecZeroEntries(F);
108:     DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
109:     DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
110:     DMRestoreLocalVector(dm,&Floc);
111:   } break;
112:   default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
113:   }
114:   DMDAVecRestoreArray(dm,Xloc,&x);
115:   DMRestoreLocalVector(dm,&Xloc);
116:   if (snes->domainerror) {
117:     VecSetInf(F);
118:   }
119:   return(0);
120: }

124: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx)
125: {
127:   DM             dm;
128:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
129:   DMDALocalInfo  info;
130:   Vec            Xloc;
131:   void           *x;

137:   if (!dmdasnes->objectivelocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
138:   SNESGetDM(snes,&dm);
139:   DMGetLocalVector(dm,&Xloc);
140:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
141:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
142:   DMDAGetLocalInfo(dm,&info);
143:   DMDAVecGetArray(dm,Xloc,&x);
144:   CHKMEMQ;
145:   (*dmdasnes->objectivelocal)(&info,x,ob,dmdasnes->objectivelocalctx);
146:   CHKMEMQ;
147:   DMDAVecRestoreArray(dm,Xloc,&x);
148:   DMRestoreLocalVector(dm,&Xloc);
149:   return(0);
150: }


155: PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
156: {
158:   DM             dm;
159:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
160:   DMDALocalInfo  info;
161:   Vec            Xloc;
162:   void           *x;

165:   if (!dmdasnes->residuallocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
166:   SNESGetDM(snes,&dm);

168:   if (dmdasnes->jacobianlocal) {
169:     DMGetLocalVector(dm,&Xloc);
170:     DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
171:     DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
172:     DMDAGetLocalInfo(dm,&info);
173:     DMDAVecGetArray(dm,Xloc,&x);
174:     CHKMEMQ;
175:     (*dmdasnes->jacobianlocal)(&info,x,A,B,dmdasnes->jacobianlocalctx);
176:     CHKMEMQ;
177:     DMDAVecRestoreArray(dm,Xloc,&x);
178:     DMRestoreLocalVector(dm,&Xloc);
179:   } else {
180:     MatFDColoring fdcoloring;
181:     PetscObjectQuery((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject*)&fdcoloring);
182:     if (!fdcoloring) {
183:       ISColoring coloring;

185:       DMCreateColoring(dm,dm->coloringtype,&coloring);
186:       MatFDColoringCreate(B,coloring,&fdcoloring);
187:       switch (dm->coloringtype) {
188:       case IS_COLORING_GLOBAL:
189:         MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))SNESComputeFunction_DMDA,dmdasnes);
190:         break;
191:       default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"No support for coloring type '%s'",ISColoringTypes[dm->coloringtype]);
192:       }
193:       PetscObjectSetOptionsPrefix((PetscObject)fdcoloring,((PetscObject)dm)->prefix);
194:       MatFDColoringSetFromOptions(fdcoloring);
195:       MatFDColoringSetUp(B,coloring,fdcoloring);
196:       ISColoringDestroy(&coloring);
197:       PetscObjectCompose((PetscObject)dm,"DMDASNES_FDCOLORING",(PetscObject)fdcoloring);
198:       PetscObjectDereference((PetscObject)fdcoloring);

200:       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
201:        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
202:        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
203:        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
204:        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
205:        */
206:       PetscObjectDereference((PetscObject)dm);
207:     }
208:     MatFDColoringApply(B,fdcoloring,X,snes);
209:   }
210:   /* This will be redundant if the user called both, but it's too common to forget. */
211:   if (A != B) {
212:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
213:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
214:   }
215:   return(0);
216: }

220: /*@C
221:    DMDASNESSetFunctionLocal - set a local residual evaluation function

223:    Logically Collective

225:    Input Arguments:
226: +  dm - DM to associate callback with
227: .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
228: .  func - local residual evaluation
229: -  ctx - optional context for local residual evaluation

231:    Calling sequence:
232:    For PetscErrorCode (*func)(DMDALocalInfo *info,void *x, void *f, void *ctx),
233: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
234: .  x - dimensional pointer to state at which to evaluate residual (e.g. PetscScalar *x or **x or ***x)
235: .  f - dimensional pointer to residual, write the residual here (e.g. PetscScalar *f or **f or ***f)
236: -  ctx - optional context passed above

238:    Level: beginner

240: .seealso: DMDASNESSetJacobianLocal(), DMSNESSetFunction(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
241: @*/
242: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
243: {
245:   DMSNES         sdm;
246:   DMSNES_DA      *dmdasnes;

250:   DMGetDMSNESWrite(dm,&sdm);
251:   DMDASNESGetContext(dm,sdm,&dmdasnes);

253:   dmdasnes->residuallocalimode = imode;
254:   dmdasnes->residuallocal      = func;
255:   dmdasnes->residuallocalctx   = ctx;

257:   DMSNESSetFunction(dm,SNESComputeFunction_DMDA,dmdasnes);
258:   if (!sdm->ops->computejacobian) {  /* Call us for the Jacobian too, can be overridden by the user. */
259:     DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
260:   }
261:   return(0);
262: }

266: /*@C
267:    DMDASNESSetJacobianLocal - set a local Jacobian evaluation function

269:    Logically Collective

271:    Input Arguments:
272: +  dm - DM to associate callback with
273: .  func - local Jacobian evaluation
274: -  ctx - optional context for local Jacobian evaluation

276:    Calling sequence:
277:    For PetscErrorCode (*func)(DMDALocalInfo *info,void *x,Mat J,Mat M,void *ctx),
278: +  info - DMDALocalInfo defining the subdomain to evaluate the Jacobian at
279: .  x - dimensional pointer to state at which to evaluate Jacobian (e.g. PetscScalar *x or **x or ***x)
280: .  J - Mat object for the Jacobian
281: .  M - Mat object for the Jacobian preconditioner matrix
282: -  ctx - optional context passed above

284:    Level: beginner

286: .seealso: DMDASNESSetFunctionLocal(), DMSNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
287: @*/
288: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
289: {
291:   DMSNES         sdm;
292:   DMSNES_DA      *dmdasnes;

296:   DMGetDMSNESWrite(dm,&sdm);
297:   DMDASNESGetContext(dm,sdm,&dmdasnes);

299:   dmdasnes->jacobianlocal    = func;
300:   dmdasnes->jacobianlocalctx = ctx;

302:   DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
303:   return(0);
304: }


309: /*@C
310:    DMDASNESSetObjectiveLocal - set a local residual evaluation function

312:    Logically Collective

314:    Input Arguments:
315: +  dm - DM to associate callback with
316: .  func - local objective evaluation
317: -  ctx - optional context for local residual evaluation

319:    Calling sequence for func:
320: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
321: .  x - dimensional pointer to state at which to evaluate residual
322: .  ob - eventual objective value
323: -  ctx - optional context passed above

325:    Level: beginner

327: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
328: @*/
329: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
330: {
332:   DMSNES         sdm;
333:   DMSNES_DA      *dmdasnes;

337:   DMGetDMSNESWrite(dm,&sdm);
338:   DMDASNESGetContext(dm,sdm,&dmdasnes);

340:   dmdasnes->objectivelocal    = func;
341:   dmdasnes->objectivelocalctx = ctx;

343:   DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
344:   return(0);
345: }

349: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
350: {
352:   DM             dm;
353:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
354:   DMDALocalInfo  info;
355:   Vec            Xloc;
356:   void           *x,*f;

362:   if (!dmdasnes->rhsplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
363:   SNESGetDM(snes,&dm);
364:   DMGetLocalVector(dm,&Xloc);
365:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
366:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
367:   DMDAGetLocalInfo(dm,&info);
368:   DMDAVecGetArray(dm,Xloc,&x);
369:   switch (dmdasnes->residuallocalimode) {
370:   case INSERT_VALUES: {
371:     DMDAVecGetArray(dm,F,&f);
372:     CHKMEMQ;
373:     (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
374:     CHKMEMQ;
375:     DMDAVecRestoreArray(dm,F,&f);
376:   } break;
377:   case ADD_VALUES: {
378:     Vec Floc;
379:     DMGetLocalVector(dm,&Floc);
380:     VecZeroEntries(Floc);
381:     DMDAVecGetArray(dm,Floc,&f);
382:     CHKMEMQ;
383:     (*dmdasnes->rhsplocal)(&info,x,f,dmdasnes->picardlocalctx);
384:     CHKMEMQ;
385:     DMDAVecRestoreArray(dm,Floc,&f);
386:     VecZeroEntries(F);
387:     DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
388:     DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
389:     DMRestoreLocalVector(dm,&Floc);
390:   } break;
391:   default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
392:   }
393:   DMDAVecRestoreArray(dm,Xloc,&x);
394:   DMRestoreLocalVector(dm,&Xloc);
395:   return(0);
396: }

400: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
401: {
403:   DM             dm;
404:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
405:   DMDALocalInfo  info;
406:   Vec            Xloc;
407:   void           *x;

410:   if (!dmdasnes->jacobianplocal) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Corrupt context");
411:   SNESGetDM(snes,&dm);

413:   DMGetLocalVector(dm,&Xloc);
414:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
415:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
416:   DMDAGetLocalInfo(dm,&info);
417:   DMDAVecGetArray(dm,Xloc,&x);
418:   CHKMEMQ;
419:   (*dmdasnes->jacobianplocal)(&info,x,A,B,dmdasnes->picardlocalctx);
420:   CHKMEMQ;
421:   DMDAVecRestoreArray(dm,Xloc,&x);
422:   DMRestoreLocalVector(dm,&Xloc);
423:   if (A != B) {
424:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
425:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
426:   }
427:   return(0);
428: }

432: /*@C
433:    DMDASNESSetPicardLocal - set a local right hand side and matrix evaluation function for Picard iteration

435:    Logically Collective

437:    Input Arguments:
438: +  dm - DM to associate callback with
439: .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
440: .  func - local residual evaluation
441: -  ctx - optional context for local residual evaluation

443:    Calling sequence for func:
444: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
445: .  x - dimensional pointer to state at which to evaluate residual
446: .  f - dimensional pointer to residual, write the residual here
447: -  ctx - optional context passed above

449:    Notes:  The user must use
450:     extern PetscErrorCode  SNESPicardComputeFunction(SNES,Vec,Vec,void*);
451:     extern PetscErrorCode  SNESPicardComputeJacobian(SNES,Vec,Mat,Mat,MatStructure*,void*);
452:     SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
453:     in their code before calling this routine.


456:    Level: beginner

458: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
459: @*/
460: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
461:                                       PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
462: {
464:   DMSNES         sdm;
465:   DMSNES_DA      *dmdasnes;

469:   DMGetDMSNESWrite(dm,&sdm);
470:   DMDASNESGetContext(dm,sdm,&dmdasnes);

472:   dmdasnes->residuallocalimode = imode;
473:   dmdasnes->rhsplocal          = func;
474:   dmdasnes->jacobianplocal     = jac;
475:   dmdasnes->picardlocalctx     = ctx;

477:   DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
478:   return(0);
479: }