Actual source code: dmdasnes.c

petsc-3.7.4 2016-10-02
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: /* Routine is called by example, hence must be labeled PETSC_EXTERN */
156: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
157: {
159:   DM             dm;
160:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
161:   DMDALocalInfo  info;
162:   Vec            Xloc;
163:   void           *x;

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

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

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

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

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

224:    Logically Collective

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

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

239:    Level: beginner

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

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

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

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

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

270:    Logically Collective

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

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

285:    Level: beginner

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

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

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

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


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

313:    Logically Collective

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

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

326:    Level: beginner

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

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

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

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

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

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

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

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

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

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

436:    Logically Collective

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

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

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


457:    Level: beginner

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

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

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

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