Actual source code: dmdasnes.c

petsc-3.6.1 2015-07-22
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:     CHKMEMQ;
 90:     (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
 91:     CHKMEMQ;
 92:     DMDAVecRestoreArray(dm,F,&f);
 93:   } break;
 94:   case ADD_VALUES: {
 95:     Vec Floc;
 96:     DMGetLocalVector(dm,&Floc);
 97:     VecZeroEntries(Floc);
 98:     DMDAVecGetArray(dm,Floc,&f);
 99:     CHKMEMQ;
100:     (*dmdasnes->residuallocal)(&info,x,f,dmdasnes->residuallocalctx);
101:     CHKMEMQ;
102:     DMDAVecRestoreArray(dm,Floc,&f);
103:     VecZeroEntries(F);
104:     DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
105:     DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
106:     DMRestoreLocalVector(dm,&Floc);
107:   } break;
108:   default: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdasnes->residuallocalimode);
109:   }
110:   DMDAVecRestoreArray(dm,Xloc,&x);
111:   DMRestoreLocalVector(dm,&Xloc);
112:   if (snes->domainerror) {
113:     VecSetInf(F);
114:   }
115:   return(0);
116: }

120: static PetscErrorCode SNESComputeObjective_DMDA(SNES snes,Vec X,PetscReal *ob,void *ctx)
121: {
123:   DM             dm;
124:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
125:   DMDALocalInfo  info;
126:   Vec            Xloc;
127:   void           *x;

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


151: PetscErrorCode SNESComputeJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
152: {
154:   DM             dm;
155:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
156:   DMDALocalInfo  info;
157:   Vec            Xloc;
158:   void           *x;

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

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

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

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

216: /*@C
217:    DMDASNESSetFunctionLocal - set a local residual evaluation function

219:    Logically Collective

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

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

234:    Level: beginner

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

246:   DMGetDMSNESWrite(dm,&sdm);
247:   DMDASNESGetContext(dm,sdm,&dmdasnes);

249:   dmdasnes->residuallocalimode = imode;
250:   dmdasnes->residuallocal      = func;
251:   dmdasnes->residuallocalctx   = ctx;

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

262: /*@C
263:    DMDASNESSetJacobianLocal - set a local Jacobian evaluation function

265:    Logically Collective

267:    Input Arguments:
268: +  dm - DM to associate callback with
269: .  func - local Jacobian evaluation
270: -  ctx - optional context for local Jacobian evaluation

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

280:    Level: beginner

282: .seealso: DMDASNESSetFunctionLocal(), DMSNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
283: @*/
284: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
285: {
287:   DMSNES         sdm;
288:   DMSNES_DA      *dmdasnes;

292:   DMGetDMSNESWrite(dm,&sdm);
293:   DMDASNESGetContext(dm,sdm,&dmdasnes);

295:   dmdasnes->jacobianlocal    = func;
296:   dmdasnes->jacobianlocalctx = ctx;

298:   DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
299:   return(0);
300: }


305: /*@C
306:    DMDASNESSetObjectiveLocal - set a local residual evaluation function

308:    Logically Collective

310:    Input Arguments:
311: +  dm - DM to associate callback with
312: .  func - local objective evaluation
313: -  ctx - optional context for local residual evaluation

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

321:    Level: beginner

323: .seealso: DMSNESSetFunction(), DMDASNESSetJacobianLocal(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
324: @*/
325: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
326: {
328:   DMSNES         sdm;
329:   DMSNES_DA      *dmdasnes;

333:   DMGetDMSNESWrite(dm,&sdm);
334:   DMDASNESGetContext(dm,sdm,&dmdasnes);

336:   dmdasnes->objectivelocal    = func;
337:   dmdasnes->objectivelocalctx = ctx;

339:   DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
340:   return(0);
341: }

345: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
346: {
348:   DM             dm;
349:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
350:   DMDALocalInfo  info;
351:   Vec            Xloc;
352:   void           *x,*f;

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

396: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
397: {
399:   DM             dm;
400:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
401:   DMDALocalInfo  info;
402:   Vec            Xloc;
403:   void           *x;

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

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

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

431:    Logically Collective

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

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

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


452:    Level: beginner

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

465:   DMGetDMSNESWrite(dm,&sdm);
466:   DMDASNESGetContext(dm,sdm,&dmdasnes);

468:   dmdasnes->residuallocalimode = imode;
469:   dmdasnes->rhsplocal          = func;
470:   dmdasnes->jacobianplocal     = jac;
471:   dmdasnes->picardlocalctx     = ctx;

473:   DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
474:   return(0);
475: }