Actual source code: dmdasnes.c

petsc-3.5.2 2014-09-08
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:   return(0);
113: }

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

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


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

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

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

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

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

213: /*@C
214:    DMDASNESSetFunctionLocal - set a local residual evaluation function

216:    Logically Collective

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

224:    Calling sequence for func:
225: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
226: .  x - dimensional pointer to state at which to evaluate residual
227: .  f - dimensional pointer to residual, write the residual here
228: -  ctx - optional context passed above

230:    Level: beginner

232: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
233: @*/
234: PetscErrorCode DMDASNESSetFunctionLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),void *ctx)
235: {
237:   DMSNES         sdm;
238:   DMSNES_DA      *dmdasnes;

242:   DMGetDMSNESWrite(dm,&sdm);
243:   DMDASNESGetContext(dm,sdm,&dmdasnes);

245:   dmdasnes->residuallocalimode = imode;
246:   dmdasnes->residuallocal      = func;
247:   dmdasnes->residuallocalctx   = ctx;

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

258: /*@C
259:    DMDASNESSetJacobianLocal - set a local Jacobian evaluation function

261:    Logically Collective

263:    Input Arguments:
264: +  dm - DM to associate callback with
265: .  func - local residual evaluation
266: -  ctx - optional context for local residual evaluation

268:    Calling sequence for func:
269: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
270: .  x - dimensional pointer to state at which to evaluate residual
271: .  f - dimensional pointer to residual, write the residual here
272: -  ctx - optional context passed above

274:    Level: beginner

276: .seealso: DMSNESSetJacobian(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
277: @*/
278: PetscErrorCode DMDASNESSetJacobianLocal(DM dm,PetscErrorCode (*func)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
279: {
281:   DMSNES         sdm;
282:   DMSNES_DA      *dmdasnes;

286:   DMGetDMSNESWrite(dm,&sdm);
287:   DMDASNESGetContext(dm,sdm,&dmdasnes);

289:   dmdasnes->jacobianlocal    = func;
290:   dmdasnes->jacobianlocalctx = ctx;

292:   DMSNESSetJacobian(dm,SNESComputeJacobian_DMDA,dmdasnes);
293:   return(0);
294: }


299: /*@C
300:    DMDASNESSetObjectiveLocal - set a local residual evaluation function

302:    Logically Collective

304:    Input Arguments:
305: +  dm - DM to associate callback with
306: .  func - local objective evaluation
307: -  ctx - optional context for local residual evaluation

309:    Calling sequence for func:
310: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
311: .  x - dimensional pointer to state at which to evaluate residual
312: .  ob - eventual objective value
313: -  ctx - optional context passed above

315:    Level: beginner

317: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
318: @*/
319: PetscErrorCode DMDASNESSetObjectiveLocal(DM dm,DMDASNESObjective func,void *ctx)
320: {
322:   DMSNES         sdm;
323:   DMSNES_DA      *dmdasnes;

327:   DMGetDMSNESWrite(dm,&sdm);
328:   DMDASNESGetContext(dm,sdm,&dmdasnes);

330:   dmdasnes->objectivelocal    = func;
331:   dmdasnes->objectivelocalctx = ctx;

333:   DMSNESSetObjective(dm,SNESComputeObjective_DMDA,dmdasnes);
334:   return(0);
335: }

339: static PetscErrorCode SNESComputePicard_DMDA(SNES snes,Vec X,Vec F,void *ctx)
340: {
342:   DM             dm;
343:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
344:   DMDALocalInfo  info;
345:   Vec            Xloc;
346:   void           *x,*f;

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

390: static PetscErrorCode SNESComputePicardJacobian_DMDA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
391: {
393:   DM             dm;
394:   DMSNES_DA      *dmdasnes = (DMSNES_DA*)ctx;
395:   DMDALocalInfo  info;
396:   Vec            Xloc;
397:   void           *x;

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

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

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

425:    Logically Collective

427:    Input Arguments:
428: +  dm - DM to associate callback with
429: .  imode - INSERT_VALUES if local function computes owned part, ADD_VALUES if it contributes to ghosted part
430: .  func - local residual evaluation
431: -  ctx - optional context for local residual evaluation

433:    Calling sequence for func:
434: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
435: .  x - dimensional pointer to state at which to evaluate residual
436: .  f - dimensional pointer to residual, write the residual here
437: -  ctx - optional context passed above

439:    Notes:  The user must use
440:     extern PetscErrorCode  SNESPicardComputeFunction(SNES,Vec,Vec,void*);
441:     extern PetscErrorCode  SNESPicardComputeJacobian(SNES,Vec,Mat,Mat,MatStructure*,void*);
442:     SNESSetFunction(snes,NULL,SNESPicardComputeFunction,&user);
443:     in their code before calling this routine.


446:    Level: beginner

448: .seealso: DMSNESSetFunction(), DMDASNESSetJacobian(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d()
449: @*/
450: PetscErrorCode DMDASNESSetPicardLocal(DM dm,InsertMode imode,PetscErrorCode (*func)(DMDALocalInfo*,void*,void*,void*),
451:                                       PetscErrorCode (*jac)(DMDALocalInfo*,void*,Mat,Mat,void*),void *ctx)
452: {
454:   DMSNES         sdm;
455:   DMSNES_DA      *dmdasnes;

459:   DMGetDMSNESWrite(dm,&sdm);
460:   DMDASNESGetContext(dm,sdm,&dmdasnes);

462:   dmdasnes->residuallocalimode = imode;
463:   dmdasnes->rhsplocal          = func;
464:   dmdasnes->jacobianplocal     = jac;
465:   dmdasnes->picardlocalctx     = ctx;

467:   DMSNESSetPicard(dm,SNESComputePicard_DMDA,SNESComputePicardJacobian_DMDA,dmdasnes);
468:   return(0);
469: }