Actual source code: dmdats.c

petsc-3.8.3 2017-12-09
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;
 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,&Xloc);
 70:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
 71:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
 72:   DMDAGetLocalInfo(dm,&info);
 73:   DMDAVecGetArray(dm,Xloc,&x);
 74:   DMDAVecGetArray(dm,Xdot,&xdot);
 75:   switch (dmdats->ifunctionlocalimode) {
 76:   case INSERT_VALUES: {
 77:     DMDAVecGetArray(dm,F,&f);
 78:     CHKMEMQ;
 79:     (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);
 80:     CHKMEMQ;
 81:     DMDAVecRestoreArray(dm,F,&f);
 82:   } break;
 83:   case ADD_VALUES: {
 84:     Vec Floc;
 85:     DMGetLocalVector(dm,&Floc);
 86:     VecZeroEntries(Floc);
 87:     DMDAVecGetArray(dm,Floc,&f);
 88:     CHKMEMQ;
 89:     (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);
 90:     CHKMEMQ;
 91:     DMDAVecRestoreArray(dm,Floc,&f);
 92:     VecZeroEntries(F);
 93:     DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
 94:     DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
 95:     DMRestoreLocalVector(dm,&Floc);
 96:   } break;
 97:   default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode);
 98:   }
 99:   DMDAVecRestoreArray(dm,Xloc,&x);
100:   DMRestoreLocalVector(dm,&Xloc);
101:   DMDAVecRestoreArray(dm,Xdot,&xdot);
102:   return(0);
103: }

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

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

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

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

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

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

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

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


223: /*@C
224:    DMDATSSetRHSFunctionLocal - set a local residual evaluation function

226:    Logically Collective

228:    Input Arguments:
229: +  dm - DM to associate callback with
230: .  imode - insert mode for the residual
231: .  func - local residual evaluation
232: -  ctx - optional context for local residual evaluation

234:    Calling sequence for func:

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

238: +  info - DMDALocalInfo defining the subdomain to evaluate the residual on
239: .  t - time at which to evaluate residual
240: .  x - array of local state information
241: .  f - output array of local residual information
242: -  ctx - optional user context

244:    Level: beginner

246: .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal()
247: @*/
248: PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx)
249: {
251:   DMTS           sdm;
252:   DMTS_DA        *dmdats;

256:   DMGetDMTSWrite(dm,&sdm);
257:   DMDATSGetContext(dm,sdm,&dmdats);
258:   dmdats->rhsfunctionlocalimode = imode;
259:   dmdats->rhsfunctionlocal      = func;
260:   dmdats->rhsfunctionlocalctx   = ctx;
261:   DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);
262:   return(0);
263: }

265: /*@C
266:    DMDATSSetRHSJacobianLocal - set a local residual evaluation function

268:    Logically Collective

270:    Input Arguments:
271: +  dm    - DM to associate callback with
272: .  func  - local RHS Jacobian evaluation routine
273: -  ctx   - optional context for local jacobian evaluation

275:    Calling sequence for func:

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

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

286:    Level: beginner

288: .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal()
289: @*/
290: PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx)
291: {
293:   DMTS           sdm;
294:   DMTS_DA        *dmdats;

298:   DMGetDMTSWrite(dm,&sdm);
299:   DMDATSGetContext(dm,sdm,&dmdats);
300:   dmdats->rhsjacobianlocal    = func;
301:   dmdats->rhsjacobianlocalctx = ctx;
302:   DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);
303:   return(0);
304: }


307: /*@C
308:    DMDATSSetIFunctionLocal - set a local residual evaluation function

310:    Logically Collective

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

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

325:    Level: beginner

327: .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal()
328: @*/
329: PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx)
330: {
332:   DMTS           sdm;
333:   DMTS_DA        *dmdats;

337:   DMGetDMTSWrite(dm,&sdm);
338:   DMDATSGetContext(dm,sdm,&dmdats);
339:   dmdats->ifunctionlocalimode = imode;
340:   dmdats->ifunctionlocal      = func;
341:   dmdats->ifunctionlocalctx   = ctx;
342:   DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);
343:   return(0);
344: }

346: /*@C
347:    DMDATSSetIJacobianLocal - set a local residual evaluation function

349:    Logically Collective

351:    Input Arguments:
352: +  dm   - DM to associate callback with
353: .  func - local residual evaluation
354: -  ctx   - optional context for local residual evaluation

356:    Calling sequence for func:

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

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

369:    Level: beginner

371: .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal()
372: @*/
373: PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx)
374: {
376:   DMTS           sdm;
377:   DMTS_DA        *dmdats;

381:   DMGetDMTSWrite(dm,&sdm);
382:   DMDATSGetContext(dm,sdm,&dmdats);
383:   dmdats->ijacobianlocal    = func;
384:   dmdats->ijacobianlocalctx = ctx;
385:   DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);
386:   return(0);
387: }

389: PetscErrorCode TSMonitorDMDARayDestroy(void **mctx)
390: {
391:   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) *mctx;
392:   PetscErrorCode       ierr;

395:   if (rayctx->lgctx) {TSMonitorLGCtxDestroy(&rayctx->lgctx);}
396:   VecDestroy(&rayctx->ray);
397:   VecScatterDestroy(&rayctx->scatter);
398:   PetscViewerDestroy(&rayctx->viewer);
399:   PetscFree(rayctx);
400:   return(0);
401: }

403: PetscErrorCode TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
404: {
405:   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)mctx;
406:   Vec                 solution;
407:   PetscErrorCode      ierr;

410:   TSGetSolution(ts,&solution);
411:   VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
412:   VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
413:   if (rayctx->viewer) {
414:     VecView(rayctx->ray,rayctx->viewer);
415:   }
416:   return(0);
417: }

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

429:   VecScatterBegin(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
430:   VecScatterEnd(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
431:   if (!step) {
432:     PetscDrawAxis axis;

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