Actual source code: dmdats.c

petsc-3.4.4 2014-03-13
  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: }