xref: /petsc/src/ts/utils/dmdats.c (revision 6ac5842e34eedc6428162d8d42bedaaf46eae34c)
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 
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;
18 
19 #undef __FUNCT__
20 #define __FUNCT__ "DMTSDestroy_DMDA"
21 static PetscErrorCode DMTSDestroy_DMDA(DMTS sdm)
22 {
23   PetscErrorCode ierr;
24 
25   PetscFunctionBegin;
26   ierr = PetscFree(sdm->data);CHKERRQ(ierr);
27   PetscFunctionReturn(0);
28 }
29 
30 #undef __FUNCT__
31 #define __FUNCT__ "DMTSDuplicate_DMDA"
32 static PetscErrorCode DMTSDuplicate_DMDA(DMTS oldsdm,DMTS sdm)
33 {
34   PetscErrorCode ierr;
35 
36   PetscFunctionBegin;
37   ierr = PetscNewLog(sdm,DMTS_DA,&sdm->data);CHKERRQ(ierr);
38   if (oldsdm->data) {ierr = PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMTS_DA));CHKERRQ(ierr);}
39   PetscFunctionReturn(0);
40 }
41 
42 #undef __FUNCT__
43 #define __FUNCT__ "DMDATSGetContext"
44 static PetscErrorCode DMDATSGetContext(DM dm,DMTS sdm,DMTS_DA **dmdats)
45 {
46   PetscErrorCode ierr;
47 
48   PetscFunctionBegin;
49   *dmdats = NULL;
50   if (!sdm->data) {
51     ierr = PetscNewLog(dm,DMTS_DA,&sdm->data);CHKERRQ(ierr);
52     sdm->ops->destroy   = DMTSDestroy_DMDA;
53     sdm->ops->duplicate = DMTSDuplicate_DMDA;
54   }
55   *dmdats = (DMTS_DA*)sdm->data;
56   PetscFunctionReturn(0);
57 }
58 
59 #undef __FUNCT__
60 #define __FUNCT__ "TSComputeIFunction_DMDA"
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 {
67   PetscErrorCode ierr;
68   DM             dm;
69   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
70   DMDALocalInfo  info;
71   Vec            Xloc;
72   void           *x,*f,*xdot;
73 
74   PetscFunctionBegin;
75   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
76   PetscValidHeaderSpecific(X,VEC_CLASSID,2);
77   PetscValidHeaderSpecific(F,VEC_CLASSID,3);
78   if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
79   ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
80   ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
81   ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
82   ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
83   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
84   ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
85   ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr);
86   switch (dmdats->ifunctionlocalimode) {
87   case INSERT_VALUES: {
88     ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);
89     CHKMEMQ;
90     ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr);
91     CHKMEMQ;
92     ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);
93   } break;
94   case ADD_VALUES: {
95     Vec Floc;
96     ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
97     ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
98     ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);
99     CHKMEMQ;
100     ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr);
101     CHKMEMQ;
102     ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);
103     ierr = VecZeroEntries(F);CHKERRQ(ierr);
104     ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
105     ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
106     ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
107   } break;
108   default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode);
109   }
110   ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
111   ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
112   ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr);
113   PetscFunctionReturn(0);
114 }
115 
116 #undef __FUNCT__
117 #define __FUNCT__ "TSComputeIJacobian_DMDA"
118 static PetscErrorCode TSComputeIJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,PetscReal shift,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
119 {
120   PetscErrorCode ierr;
121   DM             dm;
122   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
123   DMDALocalInfo  info;
124   Vec            Xloc;
125   void           *x,*xdot;
126 
127   PetscFunctionBegin;
128   if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
129   ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
130 
131   if (dmdats->ijacobianlocal) {
132     ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
133     ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
134     ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
135     ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
136     ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
137     ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr);
138     CHKMEMQ;
139     ierr = (*dmdats->ijacobianlocal)(&info,ptime,x,xdot,shift,*A,*B,mstr,dmdats->ijacobianlocalctx);CHKERRQ(ierr);
140     CHKMEMQ;
141     ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
142     ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr);
143     ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
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     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
148     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
149   }
150   PetscFunctionReturn(0);
151 }
152 
153 #undef __FUNCT__
154 #define __FUNCT__ "TSComputeRHSFunction_DMDA"
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 {
161   PetscErrorCode ierr;
162   DM             dm;
163   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
164   DMDALocalInfo  info;
165   Vec            Xloc;
166   void           *x,*f;
167 
168   PetscFunctionBegin;
169   PetscValidHeaderSpecific(ts,TS_CLASSID,1);
170   PetscValidHeaderSpecific(X,VEC_CLASSID,2);
171   PetscValidHeaderSpecific(F,VEC_CLASSID,3);
172   if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
173   ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
174   ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
175   ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
176   ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
177   ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
178   ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
179   switch (dmdats->rhsfunctionlocalimode) {
180   case INSERT_VALUES: {
181     ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);
182     CHKMEMQ;
183     ierr = (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);CHKERRQ(ierr);
184     CHKMEMQ;
185     ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);
186   } break;
187   case ADD_VALUES: {
188     Vec Floc;
189     ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
190     ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
191     ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);
192     CHKMEMQ;
193     ierr = (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);CHKERRQ(ierr);
194     CHKMEMQ;
195     ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);
196     ierr = VecZeroEntries(F);CHKERRQ(ierr);
197     ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
198     ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
199     ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
200   } break;
201   default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->rhsfunctionlocalimode);
202   }
203   ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
204   ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
205   PetscFunctionReturn(0);
206 }
207 
208 #undef __FUNCT__
209 #define __FUNCT__ "TSComputeRHSJacobian_DMDA"
210 static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
211 {
212   PetscErrorCode ierr;
213   DM             dm;
214   DMTS_DA        *dmdats = (DMTS_DA*)ctx;
215   DMDALocalInfo  info;
216   Vec            Xloc;
217   void           *x;
218 
219   PetscFunctionBegin;
220   if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
221   ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
222 
223   if (dmdats->rhsjacobianlocal) {
224     ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
225     ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
226     ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
227     ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
228     ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
229     CHKMEMQ;
230     ierr = (*dmdats->rhsjacobianlocal)(&info,ptime,x,*A,*B,mstr,dmdats->rhsjacobianlocalctx);CHKERRQ(ierr);
231     CHKMEMQ;
232     ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
233     ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
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     ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
238     ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
239   }
240   PetscFunctionReturn(0);
241 }
242 
243 
244 #undef __FUNCT__
245 #define __FUNCT__ "DMDATSSetRHSFunctionLocal"
246 /*@C
247    DMDATSSetRHSFunctionLocal - set a local residual evaluation function
248 
249    Logically Collective
250 
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
256 
257    Calling sequence for func:
258 
259 $ func(DMDALocalInfo info,PetscReal t,void *x,void *f,void *ctx)
260 
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
266 
267    Level: beginner
268 
269 .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal()
270 @*/
271 PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx)
272 {
273   PetscErrorCode ierr;
274   DMTS           sdm;
275   DMTS_DA        *dmdats;
276 
277   PetscFunctionBegin;
278   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
279   ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr);
280   ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr);
281   dmdats->rhsfunctionlocalimode = imode;
282   dmdats->rhsfunctionlocal      = func;
283   dmdats->rhsfunctionlocalctx   = ctx;
284   ierr = DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);CHKERRQ(ierr);
285   PetscFunctionReturn(0);
286 }
287 
288 #undef __FUNCT__
289 #define __FUNCT__ "DMDATSSetRHSJacobianLocal"
290 /*@C
291    DMDATSSetRHSJacobianLocal - set a local residual evaluation function
292 
293    Logically Collective
294 
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
299 
300    Calling sequence for func:
301 
302 $ func(DMDALocalInfo* info,PetscReal t,void* x,Mat J,Mat B,MatStructure *flg,void *ctx);
303 
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
311 
312    Level: beginner
313 
314 .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal()
315 @*/
316 PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx)
317 {
318   PetscErrorCode ierr;
319   DMTS           sdm;
320   DMTS_DA        *dmdats;
321 
322   PetscFunctionBegin;
323   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
324   ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr);
325   ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr);
326   dmdats->rhsjacobianlocal    = func;
327   dmdats->rhsjacobianlocalctx = ctx;
328   ierr = DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);CHKERRQ(ierr);
329   PetscFunctionReturn(0);
330 }
331 
332 
333 #undef __FUNCT__
334 #define __FUNCT__ "DMDATSSetIFunctionLocal"
335 /*@C
336    DMDATSSetIFunctionLocal - set a local residual evaluation function
337 
338    Logically Collective
339 
340    Input Arguments:
341 +  dm   - DM to associate callback with
342 .  func - local residual evaluation
343 -  ctx  - optional context for local residual evaluation
344 
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
352 
353    Level: beginner
354 
355 .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal()
356 @*/
357 PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx)
358 {
359   PetscErrorCode ierr;
360   DMTS           sdm;
361   DMTS_DA        *dmdats;
362 
363   PetscFunctionBegin;
364   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
365   ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr);
366   ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr);
367   dmdats->ifunctionlocalimode = imode;
368   dmdats->ifunctionlocal      = func;
369   dmdats->ifunctionlocalctx   = ctx;
370   ierr = DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);CHKERRQ(ierr);
371   PetscFunctionReturn(0);
372 }
373 
374 #undef __FUNCT__
375 #define __FUNCT__ "DMDATSSetIJacobianLocal"
376 /*@C
377    DMDATSSetIJacobianLocal - set a local residual evaluation function
378 
379    Logically Collective
380 
381    Input Arguments:
382 +  dm   - DM to associate callback with
383 .  func - local residual evaluation
384 -  ctx   - optional context for local residual evaluation
385 
386    Calling sequence for func:
387 
388 $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,Mat J,Mat B,MatStructure *flg,void *ctx);
389 
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
398 
399    Level: beginner
400 
401 .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal()
402 @*/
403 PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx)
404 {
405   PetscErrorCode ierr;
406   DMTS           sdm;
407   DMTS_DA        *dmdats;
408 
409   PetscFunctionBegin;
410   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
411   ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr);
412   ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr);
413   dmdats->ijacobianlocal    = func;
414   dmdats->ijacobianlocalctx = ctx;
415   ierr = DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);CHKERRQ(ierr);
416   PetscFunctionReturn(0);
417 }
418 
419 #undef __FUNCT__
420 #define __FUNCT__ "TSMonitorDMDARayDestroy"
421 PetscErrorCode TSMonitorDMDARayDestroy(void **mctx)
422 {
423   TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)*mctx;
424   PetscErrorCode      ierr;
425 
426   PetscFunctionBegin;
427   ierr = VecDestroy(&rayctx->ray);CHKERRQ(ierr);
428   ierr = VecScatterDestroy(&rayctx->scatter);CHKERRQ(ierr);
429   ierr = PetscViewerDestroy(&rayctx->viewer);CHKERRQ(ierr);
430   ierr = PetscFree(rayctx);CHKERRQ(ierr);
431   PetscFunctionReturn(0);
432 }
433 
434 #undef __FUNCT__
435 #define __FUNCT__ "TSMonitorDMDARay"
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;
441 
442   PetscFunctionBegin;
443   ierr = TSGetSolution(ts,&solution);CHKERRQ(ierr);
444   ierr = VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
445   ierr = VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
446   if (rayctx->viewer) {
447     ierr = VecView(rayctx->ray,rayctx->viewer);CHKERRQ(ierr);
448   }
449   PetscFunctionReturn(0);
450 }
451