Actual source code: dmdats.c
petsc-3.5.4 2015-05-23
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: }