Actual source code: dmdats.c
petsc-3.14.6 2021-03-30
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: }