Actual source code: dmdats.c
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: }
226: /*@C
227: DMDATSSetRHSFunctionLocal - set a local residual evaluation function
229: Logically Collective
231: Input Parameters:
232: + dm - DM to associate callback with
233: . imode - insert mode for the residual
234: . func - local residual evaluation
235: - ctx - optional context for local residual evaluation
237: Calling sequence for func:
239: $ func(DMDALocalInfo info,PetscReal t,void *x,void *f,void *ctx)
241: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
242: . t - time at which to evaluate residual
243: . x - array of local state information
244: . f - output array of local residual information
245: - ctx - optional user context
247: Level: beginner
249: .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal()
250: @*/
251: PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx)
252: {
254: DMTS sdm;
255: DMTS_DA *dmdats;
259: DMGetDMTSWrite(dm,&sdm);
260: DMDATSGetContext(dm,sdm,&dmdats);
261: dmdats->rhsfunctionlocalimode = imode;
262: dmdats->rhsfunctionlocal = func;
263: dmdats->rhsfunctionlocalctx = ctx;
264: DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);
265: return(0);
266: }
268: /*@C
269: DMDATSSetRHSJacobianLocal - set a local residual evaluation function
271: Logically Collective
273: Input Parameters:
274: + dm - DM to associate callback with
275: . func - local RHS Jacobian evaluation routine
276: - ctx - optional context for local jacobian evaluation
278: Calling sequence for func:
280: $ func(DMDALocalInfo* info,PetscReal t,void* x,Mat J,Mat B,void *ctx);
282: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
283: . t - time at which to evaluate residual
284: . x - array of local state information
285: . J - Jacobian matrix
286: . B - preconditioner matrix; often same as J
287: - ctx - optional context passed above
289: Level: beginner
291: .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal()
292: @*/
293: PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx)
294: {
296: DMTS sdm;
297: DMTS_DA *dmdats;
301: DMGetDMTSWrite(dm,&sdm);
302: DMDATSGetContext(dm,sdm,&dmdats);
303: dmdats->rhsjacobianlocal = func;
304: dmdats->rhsjacobianlocalctx = ctx;
305: DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);
306: return(0);
307: }
309: /*@C
310: DMDATSSetIFunctionLocal - set a local residual evaluation function
312: Logically Collective
314: Input Parameters:
315: + dm - DM to associate callback with
316: . func - local residual evaluation
317: - ctx - optional context for local residual evaluation
319: Calling sequence for func:
320: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
321: . t - time at which to evaluate residual
322: . x - array of local state information
323: . xdot - array of local time derivative information
324: . f - output array of local function evaluation information
325: - ctx - optional context passed above
327: Level: beginner
329: .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal()
330: @*/
331: PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx)
332: {
334: DMTS sdm;
335: DMTS_DA *dmdats;
339: DMGetDMTSWrite(dm,&sdm);
340: DMDATSGetContext(dm,sdm,&dmdats);
341: dmdats->ifunctionlocalimode = imode;
342: dmdats->ifunctionlocal = func;
343: dmdats->ifunctionlocalctx = ctx;
344: DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);
345: return(0);
346: }
348: /*@C
349: DMDATSSetIJacobianLocal - set a local residual evaluation function
351: Logically Collective
353: Input Parameters:
354: + dm - DM to associate callback with
355: . func - local residual evaluation
356: - ctx - optional context for local residual evaluation
358: Calling sequence for func:
360: $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,PetscScalar shift,Mat J,Mat B,void *ctx);
362: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
363: . t - time at which to evaluate the jacobian
364: . x - array of local state information
365: . xdot - time derivative at this state
366: . shift - see TSSetIJacobian() for the meaning of this parameter
367: . J - Jacobian matrix
368: . B - preconditioner matrix; often same as J
369: - ctx - optional context passed above
371: Level: beginner
373: .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal()
374: @*/
375: PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx)
376: {
378: DMTS sdm;
379: DMTS_DA *dmdats;
383: DMGetDMTSWrite(dm,&sdm);
384: DMDATSGetContext(dm,sdm,&dmdats);
385: dmdats->ijacobianlocal = func;
386: dmdats->ijacobianlocalctx = ctx;
387: DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);
388: return(0);
389: }
391: PetscErrorCode TSMonitorDMDARayDestroy(void **mctx)
392: {
393: TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) *mctx;
394: PetscErrorCode ierr;
397: if (rayctx->lgctx) {TSMonitorLGCtxDestroy(&rayctx->lgctx);}
398: VecDestroy(&rayctx->ray);
399: VecScatterDestroy(&rayctx->scatter);
400: PetscViewerDestroy(&rayctx->viewer);
401: PetscFree(rayctx);
402: return(0);
403: }
405: PetscErrorCode TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
406: {
407: TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)mctx;
408: Vec solution;
409: PetscErrorCode ierr;
412: TSGetSolution(ts,&solution);
413: VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
414: VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
415: if (rayctx->viewer) {
416: VecView(rayctx->ray,rayctx->viewer);
417: }
418: return(0);
419: }
421: PetscErrorCode TSMonitorLGDMDARay(TS ts, PetscInt step, PetscReal ptime, Vec u, void *ctx)
422: {
423: TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) ctx;
424: TSMonitorLGCtx lgctx = (TSMonitorLGCtx) rayctx->lgctx;
425: Vec v = rayctx->ray;
426: const PetscScalar *a;
427: PetscInt dim;
428: PetscErrorCode ierr;
431: VecScatterBegin(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
432: VecScatterEnd(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
433: if (!step) {
434: PetscDrawAxis axis;
436: PetscDrawLGGetAxis(lgctx->lg, &axis);
437: PetscDrawAxisSetLabels(axis, "Solution Ray as function of time", "Time", "Solution");
438: VecGetLocalSize(rayctx->ray, &dim);
439: PetscDrawLGSetDimension(lgctx->lg, dim);
440: PetscDrawLGReset(lgctx->lg);
441: }
442: VecGetArrayRead(v, &a);
443: #if defined(PETSC_USE_COMPLEX)
444: {
445: PetscReal *areal;
446: PetscInt i,n;
447: VecGetLocalSize(v, &n);
448: PetscMalloc1(n, &areal);
449: for (i = 0; i < n; ++i) areal[i] = PetscRealPart(a[i]);
450: PetscDrawLGAddCommonPoint(lgctx->lg, ptime, areal);
451: PetscFree(areal);
452: }
453: #else
454: PetscDrawLGAddCommonPoint(lgctx->lg, ptime, a);
455: #endif
456: VecRestoreArrayRead(v, &a);
457: if (((lgctx->howoften > 0) && (!(step % lgctx->howoften))) || ((lgctx->howoften == -1) && ts->reason)) {
458: PetscDrawLGDraw(lgctx->lg);
459: PetscDrawLGSave(lgctx->lg);
460: }
461: return(0);
462: }