Actual source code: dmdats.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  1: #include <petscdmda.h>          /*I "petscdmda.h" I*/
  2: #include <petsc-private/dmimpl.h>
  3: #include <petsc-private/tsimpl.h>   /*I "petscts.h" I*/
  4: #include <petscdraw.h>

  6: /* This structure holds the user-provided DMDA callbacks */
  7: typedef struct {
  8:   PetscErrorCode (*ifunctionlocal)(DMDALocalInfo*,PetscReal,void*,void*,void*,void*);
  9:   PetscErrorCode (*rhsfunctionlocal)(DMDALocalInfo*,PetscReal,void*,void*,void*);
 10:   PetscErrorCode (*ijacobianlocal)(DMDALocalInfo*,PetscReal,void*,void*,PetscReal,Mat,Mat,void*);
 11:   PetscErrorCode (*rhsjacobianlocal)(DMDALocalInfo*,PetscReal,void*,Mat,Mat,void*);
 12:   void       *ifunctionlocalctx;
 13:   void       *ijacobianlocalctx;
 14:   void       *rhsfunctionlocalctx;
 15:   void       *rhsjacobianlocalctx;
 16:   InsertMode ifunctionlocalimode;
 17:   InsertMode rhsfunctionlocalimode;
 18: } DMTS_DA;

 22: static PetscErrorCode DMTSDestroy_DMDA(DMTS sdm)
 23: {

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

 33: static PetscErrorCode DMTSDuplicate_DMDA(DMTS oldsdm,DMTS sdm)
 34: {

 38:   PetscNewLog(sdm,(DMTS_DA**)&sdm->data);
 39:   if (oldsdm->data) {PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMTS_DA));}
 40:   return(0);
 41: }

 45: static PetscErrorCode DMDATSGetContext(DM dm,DMTS sdm,DMTS_DA **dmdats)
 46: {

 50:   *dmdats = NULL;
 51:   if (!sdm->data) {
 52:     PetscNewLog(dm,(DMTS_DA**)&sdm->data);
 53:     sdm->ops->destroy   = DMTSDestroy_DMDA;
 54:     sdm->ops->duplicate = DMTSDuplicate_DMDA;
 55:   }
 56:   *dmdats = (DMTS_DA*)sdm->data;
 57:   return(0);
 58: }

 62: static PetscErrorCode TSComputeIFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *ctx)
 63: {
 65:   DM             dm;
 66:   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
 67:   DMDALocalInfo  info;
 68:   Vec            Xloc;
 69:   void           *x,*f,*xdot;

 75:   if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
 76:   TSGetDM(ts,&dm);
 77:   DMGetLocalVector(dm,&Xloc);
 78:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
 79:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
 80:   DMDAGetLocalInfo(dm,&info);
 81:   DMDAVecGetArray(dm,Xloc,&x);
 82:   DMDAVecGetArray(dm,Xdot,&xdot);
 83:   switch (dmdats->ifunctionlocalimode) {
 84:   case INSERT_VALUES: {
 85:     DMDAVecGetArray(dm,F,&f);
 86:     CHKMEMQ;
 87:     (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);
 88:     CHKMEMQ;
 89:     DMDAVecRestoreArray(dm,F,&f);
 90:   } break;
 91:   case ADD_VALUES: {
 92:     Vec Floc;
 93:     DMGetLocalVector(dm,&Floc);
 94:     VecZeroEntries(Floc);
 95:     DMDAVecGetArray(dm,Floc,&f);
 96:     CHKMEMQ;
 97:     (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);
 98:     CHKMEMQ;
 99:     DMDAVecRestoreArray(dm,Floc,&f);
100:     VecZeroEntries(F);
101:     DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
102:     DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
103:     DMRestoreLocalVector(dm,&Floc);
104:   } break;
105:   default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode);
106:   }
107:   DMDAVecRestoreArray(dm,Xloc,&x);
108:   DMRestoreLocalVector(dm,&Xloc);
109:   DMDAVecRestoreArray(dm,Xdot,&xdot);
110:   return(0);
111: }

115: static PetscErrorCode TSComputeIJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *ctx)
116: {
118:   DM             dm;
119:   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
120:   DMDALocalInfo  info;
121:   Vec            Xloc;
122:   void           *x,*xdot;

125:   if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
126:   TSGetDM(ts,&dm);

128:   if (dmdats->ijacobianlocal) {
129:     DMGetLocalVector(dm,&Xloc);
130:     DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
131:     DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
132:     DMDAGetLocalInfo(dm,&info);
133:     DMDAVecGetArray(dm,Xloc,&x);
134:     DMDAVecGetArray(dm,Xdot,&xdot);
135:     CHKMEMQ;
136:     (*dmdats->ijacobianlocal)(&info,ptime,x,xdot,shift,A,B,dmdats->ijacobianlocalctx);
137:     CHKMEMQ;
138:     DMDAVecRestoreArray(dm,Xloc,&x);
139:     DMDAVecRestoreArray(dm,Xdot,&xdot);
140:     DMRestoreLocalVector(dm,&Xloc);
141:   } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeIJacobian_DMDA() called without calling DMDATSSetIJacobian()");
142:   /* This will be redundant if the user called both, but it's too common to forget. */
143:   if (A != B) {
144:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
145:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
146:   }
147:   return(0);
148: }

152: static PetscErrorCode TSComputeRHSFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec F,void *ctx)
153: {
155:   DM             dm;
156:   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
157:   DMDALocalInfo  info;
158:   Vec            Xloc;
159:   void           *x,*f;

165:   if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
166:   TSGetDM(ts,&dm);
167:   DMGetLocalVector(dm,&Xloc);
168:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
169:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
170:   DMDAGetLocalInfo(dm,&info);
171:   DMDAVecGetArray(dm,Xloc,&x);
172:   switch (dmdats->rhsfunctionlocalimode) {
173:   case INSERT_VALUES: {
174:     DMDAVecGetArray(dm,F,&f);
175:     CHKMEMQ;
176:     (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);
177:     CHKMEMQ;
178:     DMDAVecRestoreArray(dm,F,&f);
179:   } break;
180:   case ADD_VALUES: {
181:     Vec Floc;
182:     DMGetLocalVector(dm,&Floc);
183:     VecZeroEntries(Floc);
184:     DMDAVecGetArray(dm,Floc,&f);
185:     CHKMEMQ;
186:     (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);
187:     CHKMEMQ;
188:     DMDAVecRestoreArray(dm,Floc,&f);
189:     VecZeroEntries(F);
190:     DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
191:     DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
192:     DMRestoreLocalVector(dm,&Floc);
193:   } break;
194:   default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->rhsfunctionlocalimode);
195:   }
196:   DMDAVecRestoreArray(dm,Xloc,&x);
197:   DMRestoreLocalVector(dm,&Xloc);
198:   return(0);
199: }

203: static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat A,Mat B,void *ctx)
204: {
206:   DM             dm;
207:   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
208:   DMDALocalInfo  info;
209:   Vec            Xloc;
210:   void           *x;

213:   if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
214:   TSGetDM(ts,&dm);

216:   if (dmdats->rhsjacobianlocal) {
217:     DMGetLocalVector(dm,&Xloc);
218:     DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
219:     DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
220:     DMDAGetLocalInfo(dm,&info);
221:     DMDAVecGetArray(dm,Xloc,&x);
222:     CHKMEMQ;
223:     (*dmdats->rhsjacobianlocal)(&info,ptime,x,A,B,dmdats->rhsjacobianlocalctx);
224:     CHKMEMQ;
225:     DMDAVecRestoreArray(dm,Xloc,&x);
226:     DMRestoreLocalVector(dm,&Xloc);
227:   } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeRHSJacobian_DMDA() called without calling DMDATSSetRHSJacobian()");
228:   /* This will be redundant if the user called both, but it's too common to forget. */
229:   if (A != B) {
230:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
231:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
232:   }
233:   return(0);
234: }


239: /*@C
240:    DMDATSSetRHSFunctionLocal - set a local residual evaluation function

242:    Logically Collective

244:    Input Arguments:
245: +  dm - DM to associate callback with
246: .  imode - insert mode for the residual
247: .  func - local residual evaluation
248: -  ctx - optional context for local residual evaluation

250:    Calling sequence for func:

252: $ func(DMDALocalInfo info,PetscReal t,void *x,void *f,void *ctx)

254: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
255: .  t - time at which to evaluate residual
256: .  x - array of local state information
257: .  f - output array of local residual information
258: -  ctx - optional user context

260:    Level: beginner

262: .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal()
263: @*/
264: PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx)
265: {
267:   DMTS           sdm;
268:   DMTS_DA        *dmdats;

272:   DMGetDMTSWrite(dm,&sdm);
273:   DMDATSGetContext(dm,sdm,&dmdats);
274:   dmdats->rhsfunctionlocalimode = imode;
275:   dmdats->rhsfunctionlocal      = func;
276:   dmdats->rhsfunctionlocalctx   = ctx;
277:   DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);
278:   return(0);
279: }

283: /*@C
284:    DMDATSSetRHSJacobianLocal - set a local residual evaluation function

286:    Logically Collective

288:    Input Arguments:
289: +  dm    - DM to associate callback with
290: .  func  - local RHS Jacobian evaluation routine
291: -  ctx   - optional context for local jacobian evaluation

293:    Calling sequence for func:

295: $ func(DMDALocalInfo* info,PetscReal t,void* x,Mat J,Mat B,MatStructure *flg,void *ctx);

297: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
298: .  t    - time at which to evaluate residual
299: .  x    - array of local state information
300: .  J    - Jacobian matrix
301: .  B    - preconditioner matrix; often same as J
302: .  flg  - flag indicating information about the preconditioner matrix structure (same as flag in KSPSetOperators())
303: -  ctx  - optional context passed above

305:    Level: beginner

307: .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal()
308: @*/
309: PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx)
310: {
312:   DMTS           sdm;
313:   DMTS_DA        *dmdats;

317:   DMGetDMTSWrite(dm,&sdm);
318:   DMDATSGetContext(dm,sdm,&dmdats);
319:   dmdats->rhsjacobianlocal    = func;
320:   dmdats->rhsjacobianlocalctx = ctx;
321:   DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);
322:   return(0);
323: }


328: /*@C
329:    DMDATSSetIFunctionLocal - set a local residual evaluation function

331:    Logically Collective

333:    Input Arguments:
334: +  dm   - DM to associate callback with
335: .  func - local residual evaluation
336: -  ctx  - optional context for local residual evaluation

338:    Calling sequence for func:
339: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
340: .  t    - time at which to evaluate residual
341: .  x    - array of local state information
342: .  xdot - array of local time derivative information
343: .  f    - output array of local function evaluation information
344: -  ctx - optional context passed above

346:    Level: beginner

348: .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal()
349: @*/
350: PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx)
351: {
353:   DMTS           sdm;
354:   DMTS_DA        *dmdats;

358:   DMGetDMTSWrite(dm,&sdm);
359:   DMDATSGetContext(dm,sdm,&dmdats);
360:   dmdats->ifunctionlocalimode = imode;
361:   dmdats->ifunctionlocal      = func;
362:   dmdats->ifunctionlocalctx   = ctx;
363:   DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);
364:   return(0);
365: }

369: /*@C
370:    DMDATSSetIJacobianLocal - set a local residual evaluation function

372:    Logically Collective

374:    Input Arguments:
375: +  dm   - DM to associate callback with
376: .  func - local residual evaluation
377: -  ctx   - optional context for local residual evaluation

379:    Calling sequence for func:

381: $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,Mat J,Mat B,MatStructure *flg,void *ctx);

383: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
384: .  t    - time at which to evaluate the jacobian
385: .  x    - array of local state information
386: .  xdot - time derivative at this state
387: .  J    - Jacobian matrix
388: .  B    - preconditioner matrix; often same as J
389: .  flg  - flag indicating information about the preconditioner matrix structure (same as flag in KSPSetOperators())
390: -  ctx  - optional context passed above

392:    Level: beginner

394: .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal()
395: @*/
396: PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx)
397: {
399:   DMTS           sdm;
400:   DMTS_DA        *dmdats;

404:   DMGetDMTSWrite(dm,&sdm);
405:   DMDATSGetContext(dm,sdm,&dmdats);
406:   dmdats->ijacobianlocal    = func;
407:   dmdats->ijacobianlocalctx = ctx;
408:   DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);
409:   return(0);
410: }

414: PetscErrorCode TSMonitorDMDARayDestroy(void **mctx)
415: {
416:   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) *mctx;
417:   PetscErrorCode       ierr;

420:   if (rayctx->lgctx) {TSMonitorLGCtxDestroy(&rayctx->lgctx);}
421:   VecDestroy(&rayctx->ray);
422:   VecScatterDestroy(&rayctx->scatter);
423:   PetscViewerDestroy(&rayctx->viewer);
424:   PetscFree(rayctx);
425:   return(0);
426: }

430: PetscErrorCode TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
431: {
432:   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)mctx;
433:   Vec                 solution;
434:   PetscErrorCode      ierr;

437:   TSGetSolution(ts,&solution);
438:   VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
439:   VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
440:   if (rayctx->viewer) {
441:     VecView(rayctx->ray,rayctx->viewer);
442:   }
443:   return(0);
444: }

448: PetscErrorCode  TSMonitorLGDMDARay(TS ts, PetscInt step, PetscReal ptime, Vec u, void *ctx)
449: {
450:   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) ctx;
451:   TSMonitorLGCtx       lgctx  = (TSMonitorLGCtx) rayctx->lgctx;
452:   Vec                  v      = rayctx->ray;
453:   const PetscScalar   *a;
454:   PetscInt             dim;
455:   PetscErrorCode       ierr;

458:   VecScatterBegin(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
459:   VecScatterEnd(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
460:   if (!step) {
461:     PetscDrawAxis axis;

463:     PetscDrawLGGetAxis(lgctx->lg, &axis);
464:     PetscDrawAxisSetLabels(axis, "Solution Ray as function of time", "Time", "Solution");
465:     VecGetLocalSize(rayctx->ray, &dim);
466:     PetscDrawLGSetDimension(lgctx->lg, dim);
467:     PetscDrawLGReset(lgctx->lg);
468:   }
469:   VecGetArrayRead(v, &a);
470: #if defined(PETSC_USE_COMPLEX)
471:   {
472:     PetscReal *areal;
473:     PetscInt   i,n;
474:     VecGetLocalSize(v, &n);
475:     PetscMalloc1(n, &areal);
476:     for (i = 0; i < n; ++i) areal[i] = PetscRealPart(a[i]);
477:     PetscDrawLGAddCommonPoint(lgctx->lg, ptime, areal);
478:     PetscFree(areal);
479:   }
480: #else
481:   PetscDrawLGAddCommonPoint(lgctx->lg, ptime, a);
482: #endif
483:   VecRestoreArrayRead(v, &a);
484:   if (((lgctx->howoften > 0) && (!(step % lgctx->howoften))) || ((lgctx->howoften == -1) && ts->reason)) {
485:     PetscDrawLGDraw(lgctx->lg);
486:   }
487:   return(0);
488: }