Actual source code: dmdats.c
petsc-3.4.5 2014-06-29
1: #include <petscdmda.h> /*I "petscdmda.h" I*/
2: #include <petsc-private/dmimpl.h>
3: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
5: /* This structure holds the user-provided DMDA callbacks */
6: typedef struct {
7: PetscErrorCode (*ifunctionlocal)(DMDALocalInfo*,PetscReal,void*,void*,void*,void*);
8: PetscErrorCode (*rhsfunctionlocal)(DMDALocalInfo*,PetscReal,void*,void*,void*);
9: PetscErrorCode (*ijacobianlocal)(DMDALocalInfo*,PetscReal,void*,void*,PetscReal,Mat,Mat,MatStructure*,void*);
10: PetscErrorCode (*rhsjacobianlocal)(DMDALocalInfo*,PetscReal,void*,Mat,Mat,MatStructure*,void*);
11: void *ifunctionlocalctx;
12: void *ijacobianlocalctx;
13: void *rhsfunctionlocalctx;
14: void *rhsjacobianlocalctx;
15: InsertMode ifunctionlocalimode;
16: InsertMode rhsfunctionlocalimode;
17: } DMTS_DA;
21: static PetscErrorCode DMTSDestroy_DMDA(DMTS sdm)
22: {
26: PetscFree(sdm->data);
27: return(0);
28: }
32: static PetscErrorCode DMTSDuplicate_DMDA(DMTS oldsdm,DMTS sdm)
33: {
37: PetscNewLog(sdm,DMTS_DA,&sdm->data);
38: if (oldsdm->data) {PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMTS_DA));}
39: return(0);
40: }
44: static PetscErrorCode DMDATSGetContext(DM dm,DMTS sdm,DMTS_DA **dmdats)
45: {
49: *dmdats = NULL;
50: if (!sdm->data) {
51: PetscNewLog(dm,DMTS_DA,&sdm->data);
52: sdm->ops->destroy = DMTSDestroy_DMDA;
53: sdm->ops->duplicate = DMTSDuplicate_DMDA;
54: }
55: *dmdats = (DMTS_DA*)sdm->data;
56: return(0);
57: }
61: /*
62: This function should eventually replace:
63: DMDAComputeFunction() and DMDAComputeFunction1()
64: */
65: static PetscErrorCode TSComputeIFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *ctx)
66: {
68: DM dm;
69: DMTS_DA *dmdats = (DMTS_DA*)ctx;
70: DMDALocalInfo info;
71: Vec Xloc;
72: void *x,*f,*xdot;
78: if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
79: TSGetDM(ts,&dm);
80: DMGetLocalVector(dm,&Xloc);
81: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
82: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
83: DMDAGetLocalInfo(dm,&info);
84: DMDAVecGetArray(dm,Xloc,&x);
85: DMDAVecGetArray(dm,Xdot,&xdot);
86: switch (dmdats->ifunctionlocalimode) {
87: case INSERT_VALUES: {
88: DMDAVecGetArray(dm,F,&f);
89: CHKMEMQ;
90: (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);
91: CHKMEMQ;
92: DMDAVecRestoreArray(dm,F,&f);
93: } break;
94: case ADD_VALUES: {
95: Vec Floc;
96: DMGetLocalVector(dm,&Floc);
97: VecZeroEntries(Floc);
98: DMDAVecGetArray(dm,Floc,&f);
99: CHKMEMQ;
100: (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);
101: CHKMEMQ;
102: DMDAVecRestoreArray(dm,Floc,&f);
103: VecZeroEntries(F);
104: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
105: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
106: DMRestoreLocalVector(dm,&Floc);
107: } break;
108: default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode);
109: }
110: DMDAVecRestoreArray(dm,Xloc,&x);
111: DMRestoreLocalVector(dm,&Xloc);
112: DMDAVecRestoreArray(dm,Xdot,&xdot);
113: return(0);
114: }
118: static PetscErrorCode TSComputeIJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
119: {
121: DM dm;
122: DMTS_DA *dmdats = (DMTS_DA*)ctx;
123: DMDALocalInfo info;
124: Vec Xloc;
125: void *x,*xdot;
128: if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
129: TSGetDM(ts,&dm);
131: if (dmdats->ijacobianlocal) {
132: DMGetLocalVector(dm,&Xloc);
133: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
134: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
135: DMDAGetLocalInfo(dm,&info);
136: DMDAVecGetArray(dm,Xloc,&x);
137: DMDAVecGetArray(dm,Xdot,&xdot);
138: CHKMEMQ;
139: (*dmdats->ijacobianlocal)(&info,ptime,x,xdot,shift,*A,*B,mstr,dmdats->ijacobianlocalctx);
140: CHKMEMQ;
141: DMDAVecRestoreArray(dm,Xloc,&x);
142: DMDAVecRestoreArray(dm,Xdot,&xdot);
143: DMRestoreLocalVector(dm,&Xloc);
144: } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeIJacobian_DMDA() called without calling DMDATSSetIJacobian()");
145: /* This will be redundant if the user called both, but it's too common to forget. */
146: if (*A != *B) {
147: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
148: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
149: }
150: return(0);
151: }
155: /*
156: This function should eventually replace:
157: DMDAComputeFunction() and DMDAComputeFunction1()
158: */
159: static PetscErrorCode TSComputeRHSFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec F,void *ctx)
160: {
162: DM dm;
163: DMTS_DA *dmdats = (DMTS_DA*)ctx;
164: DMDALocalInfo info;
165: Vec Xloc;
166: void *x,*f;
172: if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
173: TSGetDM(ts,&dm);
174: DMGetLocalVector(dm,&Xloc);
175: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
176: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
177: DMDAGetLocalInfo(dm,&info);
178: DMDAVecGetArray(dm,Xloc,&x);
179: switch (dmdats->rhsfunctionlocalimode) {
180: case INSERT_VALUES: {
181: DMDAVecGetArray(dm,F,&f);
182: CHKMEMQ;
183: (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);
184: CHKMEMQ;
185: DMDAVecRestoreArray(dm,F,&f);
186: } break;
187: case ADD_VALUES: {
188: Vec Floc;
189: DMGetLocalVector(dm,&Floc);
190: VecZeroEntries(Floc);
191: DMDAVecGetArray(dm,Floc,&f);
192: CHKMEMQ;
193: (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);
194: CHKMEMQ;
195: DMDAVecRestoreArray(dm,Floc,&f);
196: VecZeroEntries(F);
197: DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);
198: DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);
199: DMRestoreLocalVector(dm,&Floc);
200: } break;
201: default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->rhsfunctionlocalimode);
202: }
203: DMDAVecRestoreArray(dm,Xloc,&x);
204: DMRestoreLocalVector(dm,&Xloc);
205: return(0);
206: }
210: static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
211: {
213: DM dm;
214: DMTS_DA *dmdats = (DMTS_DA*)ctx;
215: DMDALocalInfo info;
216: Vec Xloc;
217: void *x;
220: if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
221: TSGetDM(ts,&dm);
223: if (dmdats->rhsjacobianlocal) {
224: DMGetLocalVector(dm,&Xloc);
225: DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
226: DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);
227: DMDAGetLocalInfo(dm,&info);
228: DMDAVecGetArray(dm,Xloc,&x);
229: CHKMEMQ;
230: (*dmdats->rhsjacobianlocal)(&info,ptime,x,*A,*B,mstr,dmdats->rhsjacobianlocalctx);
231: CHKMEMQ;
232: DMDAVecRestoreArray(dm,Xloc,&x);
233: DMRestoreLocalVector(dm,&Xloc);
234: } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeRHSJacobian_DMDA() called without calling DMDATSSetRHSJacobian()");
235: /* This will be redundant if the user called both, but it's too common to forget. */
236: if (*A != *B) {
237: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
238: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
239: }
240: return(0);
241: }
246: /*@C
247: DMDATSSetRHSFunctionLocal - set a local residual evaluation function
249: Logically Collective
251: Input Arguments:
252: + dm - DM to associate callback with
253: . imode - insert mode for the residual
254: . func - local residual evaluation
255: - ctx - optional context for local residual evaluation
257: Calling sequence for func:
259: $ func(DMDALocalInfo info,PetscReal t,void *x,void *f,void *ctx)
261: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
262: . t - time at which to evaluate residual
263: . x - array of local state information
264: . f - output array of local residual information
265: - ctx - optional user context
267: Level: beginner
269: .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal()
270: @*/
271: PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx)
272: {
274: DMTS sdm;
275: DMTS_DA *dmdats;
279: DMGetDMTSWrite(dm,&sdm);
280: DMDATSGetContext(dm,sdm,&dmdats);
281: dmdats->rhsfunctionlocalimode = imode;
282: dmdats->rhsfunctionlocal = func;
283: dmdats->rhsfunctionlocalctx = ctx;
284: DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);
285: return(0);
286: }
290: /*@C
291: DMDATSSetRHSJacobianLocal - set a local residual evaluation function
293: Logically Collective
295: Input Arguments:
296: + dm - DM to associate callback with
297: . func - local RHS Jacobian evaluation routine
298: - ctx - optional context for local jacobian evaluation
300: Calling sequence for func:
302: $ func(DMDALocalInfo* info,PetscReal t,void* x,Mat J,Mat B,MatStructure *flg,void *ctx);
304: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
305: . t - time at which to evaluate residual
306: . x - array of local state information
307: . J - Jacobian matrix
308: . B - preconditioner matrix; often same as J
309: . flg - flag indicating information about the preconditioner matrix structure (same as flag in KSPSetOperators())
310: - ctx - optional context passed above
312: Level: beginner
314: .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal()
315: @*/
316: PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx)
317: {
319: DMTS sdm;
320: DMTS_DA *dmdats;
324: DMGetDMTSWrite(dm,&sdm);
325: DMDATSGetContext(dm,sdm,&dmdats);
326: dmdats->rhsjacobianlocal = func;
327: dmdats->rhsjacobianlocalctx = ctx;
328: DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);
329: return(0);
330: }
335: /*@C
336: DMDATSSetIFunctionLocal - set a local residual evaluation function
338: Logically Collective
340: Input Arguments:
341: + dm - DM to associate callback with
342: . func - local residual evaluation
343: - ctx - optional context for local residual evaluation
345: Calling sequence for func:
346: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
347: . t - time at which to evaluate residual
348: . x - array of local state information
349: . xdot - array of local time derivative information
350: . f - output array of local function evaluation information
351: - ctx - optional context passed above
353: Level: beginner
355: .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal()
356: @*/
357: PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx)
358: {
360: DMTS sdm;
361: DMTS_DA *dmdats;
365: DMGetDMTSWrite(dm,&sdm);
366: DMDATSGetContext(dm,sdm,&dmdats);
367: dmdats->ifunctionlocalimode = imode;
368: dmdats->ifunctionlocal = func;
369: dmdats->ifunctionlocalctx = ctx;
370: DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);
371: return(0);
372: }
376: /*@C
377: DMDATSSetIJacobianLocal - set a local residual evaluation function
379: Logically Collective
381: Input Arguments:
382: + dm - DM to associate callback with
383: . func - local residual evaluation
384: - ctx - optional context for local residual evaluation
386: Calling sequence for func:
388: $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,Mat J,Mat B,MatStructure *flg,void *ctx);
390: + info - DMDALocalInfo defining the subdomain to evaluate the residual on
391: . t - time at which to evaluate the jacobian
392: . x - array of local state information
393: . xdot - time derivative at this state
394: . J - Jacobian matrix
395: . B - preconditioner matrix; often same as J
396: . flg - flag indicating information about the preconditioner matrix structure (same as flag in KSPSetOperators())
397: - ctx - optional context passed above
399: Level: beginner
401: .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal()
402: @*/
403: PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx)
404: {
406: DMTS sdm;
407: DMTS_DA *dmdats;
411: DMGetDMTSWrite(dm,&sdm);
412: DMDATSGetContext(dm,sdm,&dmdats);
413: dmdats->ijacobianlocal = func;
414: dmdats->ijacobianlocalctx = ctx;
415: DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);
416: return(0);
417: }
421: PetscErrorCode TSMonitorDMDARayDestroy(void **mctx)
422: {
423: TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)*mctx;
424: PetscErrorCode ierr;
427: VecDestroy(&rayctx->ray);
428: VecScatterDestroy(&rayctx->scatter);
429: PetscViewerDestroy(&rayctx->viewer);
430: PetscFree(rayctx);
431: return(0);
432: }
436: PetscErrorCode TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx)
437: {
438: TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)mctx;
439: Vec solution;
440: PetscErrorCode ierr;
443: TSGetSolution(ts,&solution);
444: VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
445: VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);
446: if (rayctx->viewer) {
447: VecView(rayctx->ray,rayctx->viewer);
448: }
449: return(0);
450: }