Actual source code: dmdats.c
petsc-3.8.4 2018-03-24
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: }