Actual source code: dmdats.c

petsc-master 2019-06-18
Report Typos and Errors
  1:  #include <petscdmda.h>
  2:  #include <petsc/private/dmimpl.h>
  3:  #include <petsc/private/tsimpl.h>
  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;

 20: static PetscErrorCode DMTSDestroy_DMDA(DMTS sdm)
 21: {

 25:   PetscFree(sdm->data);
 26:   return(0);
 27: }

 29: static PetscErrorCode DMTSDuplicate_DMDA(DMTS oldsdm,DMTS sdm)
 30: {

 34:   PetscNewLog(sdm,(DMTS_DA**)&sdm->data);
 35:   if (oldsdm->data) {PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMTS_DA));}
 36:   return(0);
 37: }

 39: static PetscErrorCode DMDATSGetContext(DM dm,DMTS sdm,DMTS_DA **dmdats)
 40: {

 44:   *dmdats = NULL;
 45:   if (!sdm->data) {
 46:     PetscNewLog(dm,(DMTS_DA**)&sdm->data);
 47:     sdm->ops->destroy   = DMTSDestroy_DMDA;
 48:     sdm->ops->duplicate = DMTSDuplicate_DMDA;
 49:   }
 50:   *dmdats = (DMTS_DA*)sdm->data;
 51:   return(0);
 52: }

 54: static PetscErrorCode TSComputeIFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *ctx)
 55: {
 57:   DM             dm;
 58:   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
 59:   DMDALocalInfo  info;
 60:   Vec            Xloc,Xdotloc;
 61:   void           *x,*f,*xdot;

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

109: static PetscErrorCode TSComputeIJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *ctx)
110: {
112:   DM             dm;
113:   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
114:   DMDALocalInfo  info;
115:   Vec            Xloc;
116:   void           *x,*xdot;

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

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

144: static PetscErrorCode TSComputeRHSFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec F,void *ctx)
145: {
147:   DM             dm;
148:   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
149:   DMDALocalInfo  info;
150:   Vec            Xloc;
151:   void           *x,*f;

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

193: static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat A,Mat B,void *ctx)
194: {
196:   DM             dm;
197:   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
198:   DMDALocalInfo  info;
199:   Vec            Xloc;
200:   void           *x;

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

206:   if (dmdats->rhsjacobianlocal) {
207:     DMGetLocalVector(dm,&Xloc);
208:     DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
209:     DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
210:     DMDAGetLocalInfo(dm,&info);
211:     DMDAVecGetArray(dm,Xloc,&x);
212:     CHKMEMQ;
213:     (*dmdats->rhsjacobianlocal)(&info,ptime,x,A,B,dmdats->rhsjacobianlocalctx);
214:     CHKMEMQ;
215:     DMDAVecRestoreArray(dm,Xloc,&x);
216:     DMRestoreLocalVector(dm,&Xloc);
217:   } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeRHSJacobian_DMDA() called without calling DMDATSSetRHSJacobian()");
218:   /* This will be redundant if the user called both, but it's too common to forget. */
219:   if (A != B) {
220:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
221:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
222:   }
223:   return(0);
224: }


227: /*@C
228:    DMDATSSetRHSFunctionLocal - set a local residual evaluation function

230:    Logically Collective

232:    Input Arguments:
233: +  dm - DM to associate callback with
234: .  imode - insert mode for the residual
235: .  func - local residual evaluation
236: -  ctx - optional context for local residual evaluation

238:    Calling sequence for func:

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

242: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
243: .  t - time at which to evaluate residual
244: .  x - array of local state information
245: .  f - output array of local residual information
246: -  ctx - optional user context

248:    Level: beginner

250: .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal()
251: @*/
252: PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx)
253: {
255:   DMTS           sdm;
256:   DMTS_DA        *dmdats;

260:   DMGetDMTSWrite(dm,&sdm);
261:   DMDATSGetContext(dm,sdm,&dmdats);
262:   dmdats->rhsfunctionlocalimode = imode;
263:   dmdats->rhsfunctionlocal      = func;
264:   dmdats->rhsfunctionlocalctx   = ctx;
265:   DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);
266:   return(0);
267: }

269: /*@C
270:    DMDATSSetRHSJacobianLocal - set a local residual evaluation function

272:    Logically Collective

274:    Input Arguments:
275: +  dm    - DM to associate callback with
276: .  func  - local RHS Jacobian evaluation routine
277: -  ctx   - optional context for local jacobian evaluation

279:    Calling sequence for func:

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

283: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
284: .  t    - time at which to evaluate residual
285: .  x    - array of local state information
286: .  J    - Jacobian matrix
287: .  B    - preconditioner matrix; often same as J
288: -  ctx  - optional context passed above

290:    Level: beginner

292: .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal()
293: @*/
294: PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx)
295: {
297:   DMTS           sdm;
298:   DMTS_DA        *dmdats;

302:   DMGetDMTSWrite(dm,&sdm);
303:   DMDATSGetContext(dm,sdm,&dmdats);
304:   dmdats->rhsjacobianlocal    = func;
305:   dmdats->rhsjacobianlocalctx = ctx;
306:   DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);
307:   return(0);
308: }


311: /*@C
312:    DMDATSSetIFunctionLocal - set a local residual evaluation function

314:    Logically Collective

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

321:    Calling sequence for func:
322: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
323: .  t    - time at which to evaluate residual
324: .  x    - array of local state information
325: .  xdot - array of local time derivative information
326: .  f    - output array of local function evaluation information
327: -  ctx - optional context passed above

329:    Level: beginner

331: .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal()
332: @*/
333: PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx)
334: {
336:   DMTS           sdm;
337:   DMTS_DA        *dmdats;

341:   DMGetDMTSWrite(dm,&sdm);
342:   DMDATSGetContext(dm,sdm,&dmdats);
343:   dmdats->ifunctionlocalimode = imode;
344:   dmdats->ifunctionlocal      = func;
345:   dmdats->ifunctionlocalctx   = ctx;
346:   DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);
347:   return(0);
348: }

350: /*@C
351:    DMDATSSetIJacobianLocal - set a local residual evaluation function

353:    Logically Collective

355:    Input Arguments:
356: +  dm   - DM to associate callback with
357: .  func - local residual evaluation
358: -  ctx   - optional context for local residual evaluation

360:    Calling sequence for func:

362: $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,PetscScalar shift,Mat J,Mat B,void *ctx);

364: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
365: .  t    - time at which to evaluate the jacobian
366: .  x    - array of local state information
367: .  xdot - time derivative at this state
368: .  shift - see TSSetIJacobian() for the meaning of this parameter
369: .  J    - Jacobian matrix
370: .  B    - preconditioner matrix; often same as J
371: -  ctx  - optional context passed above

373:    Level: beginner

375: .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal()
376: @*/
377: PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx)
378: {
380:   DMTS           sdm;
381:   DMTS_DA        *dmdats;

385:   DMGetDMTSWrite(dm,&sdm);
386:   DMDATSGetContext(dm,sdm,&dmdats);
387:   dmdats->ijacobianlocal    = func;
388:   dmdats->ijacobianlocalctx = ctx;
389:   DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);
390:   return(0);
391: }

393: PetscErrorCode TSMonitorDMDARayDestroy(void **mctx)
394: {
395:   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) *mctx;
396:   PetscErrorCode       ierr;

399:   if (rayctx->lgctx) {TSMonitorLGCtxDestroy(&rayctx->lgctx);}
400:   VecDestroy(&rayctx->ray);
401:   VecScatterDestroy(&rayctx->scatter);
402:   PetscViewerDestroy(&rayctx->viewer);
403:   PetscFree(rayctx);
404:   return(0);
405: }

407: PetscErrorCode TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
408: {
409:   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)mctx;
410:   Vec                 solution;
411:   PetscErrorCode      ierr;

414:   TSGetSolution(ts,&solution);
415:   VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
416:   VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
417:   if (rayctx->viewer) {
418:     VecView(rayctx->ray,rayctx->viewer);
419:   }
420:   return(0);
421: }

423: PetscErrorCode  TSMonitorLGDMDARay(TS ts, PetscInt step, PetscReal ptime, Vec u, void *ctx)
424: {
425:   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) ctx;
426:   TSMonitorLGCtx       lgctx  = (TSMonitorLGCtx) rayctx->lgctx;
427:   Vec                  v      = rayctx->ray;
428:   const PetscScalar   *a;
429:   PetscInt             dim;
430:   PetscErrorCode       ierr;

433:   VecScatterBegin(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
434:   VecScatterEnd(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
435:   if (!step) {
436:     PetscDrawAxis axis;

438:     PetscDrawLGGetAxis(lgctx->lg, &axis);
439:     PetscDrawAxisSetLabels(axis, "Solution Ray as function of time", "Time", "Solution");
440:     VecGetLocalSize(rayctx->ray, &dim);
441:     PetscDrawLGSetDimension(lgctx->lg, dim);
442:     PetscDrawLGReset(lgctx->lg);
443:   }
444:   VecGetArrayRead(v, &a);
445: #if defined(PETSC_USE_COMPLEX)
446:   {
447:     PetscReal *areal;
448:     PetscInt   i,n;
449:     VecGetLocalSize(v, &n);
450:     PetscMalloc1(n, &areal);
451:     for (i = 0; i < n; ++i) areal[i] = PetscRealPart(a[i]);
452:     PetscDrawLGAddCommonPoint(lgctx->lg, ptime, areal);
453:     PetscFree(areal);
454:   }
455: #else
456:   PetscDrawLGAddCommonPoint(lgctx->lg, ptime, a);
457: #endif
458:   VecRestoreArrayRead(v, &a);
459:   if (((lgctx->howoften > 0) && (!(step % lgctx->howoften))) || ((lgctx->howoften == -1) && ts->reason)) {
460:     PetscDrawLGDraw(lgctx->lg);
461:     PetscDrawLGSave(lgctx->lg);
462:   }
463:   return(0);
464: }