Actual source code: dmdats.c
petsc-3.7.3 2016-08-01
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,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: - ctx - optional context passed above
304: Level: beginner
306: .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal()
307: @*/
308: PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx)
309: {
311: DMTS sdm;
312: DMTS_DA *dmdats;
316: DMGetDMTSWrite(dm,&sdm);
317: DMDATSGetContext(dm,sdm,&dmdats);
318: dmdats->rhsjacobianlocal = func;
319: dmdats->rhsjacobianlocalctx = ctx;
320: DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);
321: return(0);
322: }
327: /*@C
328: DMDATSSetIFunctionLocal - set a local residual evaluation function
330: Logically Collective
332: Input Arguments:
333: + dm - DM to associate callback with
334: . func - local residual evaluation
335: - ctx - optional context for local residual evaluation
337: Calling sequence for func:
338: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
339: . t - time at which to evaluate residual
340: . x - array of local state information
341: . xdot - array of local time derivative information
342: . f - output array of local function evaluation information
343: - ctx - optional context passed above
345: Level: beginner
347: .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal()
348: @*/
349: PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx)
350: {
352: DMTS sdm;
353: DMTS_DA *dmdats;
357: DMGetDMTSWrite(dm,&sdm);
358: DMDATSGetContext(dm,sdm,&dmdats);
359: dmdats->ifunctionlocalimode = imode;
360: dmdats->ifunctionlocal = func;
361: dmdats->ifunctionlocalctx = ctx;
362: DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);
363: return(0);
364: }
368: /*@C
369: DMDATSSetIJacobianLocal - set a local residual evaluation function
371: Logically Collective
373: Input Arguments:
374: + dm - DM to associate callback with
375: . func - local residual evaluation
376: - ctx - optional context for local residual evaluation
378: Calling sequence for func:
380: $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,PetscScalar shift,Mat J,Mat B,void *ctx);
382: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
383: . t - time at which to evaluate the jacobian
384: . x - array of local state information
385: . xdot - time derivative at this state
386: . shift - see TSSetIJacobian() for the meaning of this parameter
387: . J - Jacobian matrix
388: . B - preconditioner matrix; often same as J
389: - ctx - optional context passed above
391: Level: beginner
393: .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal()
394: @*/
395: PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx)
396: {
398: DMTS sdm;
399: DMTS_DA *dmdats;
403: DMGetDMTSWrite(dm,&sdm);
404: DMDATSGetContext(dm,sdm,&dmdats);
405: dmdats->ijacobianlocal = func;
406: dmdats->ijacobianlocalctx = ctx;
407: DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);
408: return(0);
409: }
413: PetscErrorCode TSMonitorDMDARayDestroy(void **mctx)
414: {
415: TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) *mctx;
416: PetscErrorCode ierr;
419: if (rayctx->lgctx) {TSMonitorLGCtxDestroy(&rayctx->lgctx);}
420: VecDestroy(&rayctx->ray);
421: VecScatterDestroy(&rayctx->scatter);
422: PetscViewerDestroy(&rayctx->viewer);
423: PetscFree(rayctx);
424: return(0);
425: }
429: PetscErrorCode TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
430: {
431: TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)mctx;
432: Vec solution;
433: PetscErrorCode ierr;
436: TSGetSolution(ts,&solution);
437: VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
438: VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
439: if (rayctx->viewer) {
440: VecView(rayctx->ray,rayctx->viewer);
441: }
442: return(0);
443: }
447: PetscErrorCode TSMonitorLGDMDARay(TS ts, PetscInt step, PetscReal ptime, Vec u, void *ctx)
448: {
449: TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) ctx;
450: TSMonitorLGCtx lgctx = (TSMonitorLGCtx) rayctx->lgctx;
451: Vec v = rayctx->ray;
452: const PetscScalar *a;
453: PetscInt dim;
454: PetscErrorCode ierr;
457: VecScatterBegin(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
458: VecScatterEnd(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);
459: if (!step) {
460: PetscDrawAxis axis;
462: PetscDrawLGGetAxis(lgctx->lg, &axis);
463: PetscDrawAxisSetLabels(axis, "Solution Ray as function of time", "Time", "Solution");
464: VecGetLocalSize(rayctx->ray, &dim);
465: PetscDrawLGSetDimension(lgctx->lg, dim);
466: PetscDrawLGReset(lgctx->lg);
467: }
468: VecGetArrayRead(v, &a);
469: #if defined(PETSC_USE_COMPLEX)
470: {
471: PetscReal *areal;
472: PetscInt i,n;
473: VecGetLocalSize(v, &n);
474: PetscMalloc1(n, &areal);
475: for (i = 0; i < n; ++i) areal[i] = PetscRealPart(a[i]);
476: PetscDrawLGAddCommonPoint(lgctx->lg, ptime, areal);
477: PetscFree(areal);
478: }
479: #else
480: PetscDrawLGAddCommonPoint(lgctx->lg, ptime, a);
481: #endif
482: VecRestoreArrayRead(v, &a);
483: if (((lgctx->howoften > 0) && (!(step % lgctx->howoften))) || ((lgctx->howoften == -1) && ts->reason)) {
484: PetscDrawLGDraw(lgctx->lg);
485: PetscDrawLGSave(lgctx->lg);
486: }
487: return(0);
488: }