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>
5
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;
19
DMTSDestroy_DMDA(DMTS sdm)20 static PetscErrorCode DMTSDestroy_DMDA(DMTS sdm)
21 {
22 PetscErrorCode ierr;
23
24 PetscFunctionBegin;
25 ierr = PetscFree(sdm->data);CHKERRQ(ierr);
26 PetscFunctionReturn(0);
27 }
28
DMTSDuplicate_DMDA(DMTS oldsdm,DMTS sdm)29 static PetscErrorCode DMTSDuplicate_DMDA(DMTS oldsdm,DMTS sdm)
30 {
31 PetscErrorCode ierr;
32
33 PetscFunctionBegin;
34 ierr = PetscNewLog(sdm,(DMTS_DA**)&sdm->data);CHKERRQ(ierr);
35 if (oldsdm->data) {ierr = PetscMemcpy(sdm->data,oldsdm->data,sizeof(DMTS_DA));CHKERRQ(ierr);}
36 PetscFunctionReturn(0);
37 }
38
DMDATSGetContext(DM dm,DMTS sdm,DMTS_DA ** dmdats)39 static PetscErrorCode DMDATSGetContext(DM dm,DMTS sdm,DMTS_DA **dmdats)
40 {
41 PetscErrorCode ierr;
42
43 PetscFunctionBegin;
44 *dmdats = NULL;
45 if (!sdm->data) {
46 ierr = PetscNewLog(dm,(DMTS_DA**)&sdm->data);CHKERRQ(ierr);
47 sdm->ops->destroy = DMTSDestroy_DMDA;
48 sdm->ops->duplicate = DMTSDuplicate_DMDA;
49 }
50 *dmdats = (DMTS_DA*)sdm->data;
51 PetscFunctionReturn(0);
52 }
53
TSComputeIFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void * ctx)54 static PetscErrorCode TSComputeIFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *ctx)
55 {
56 PetscErrorCode ierr;
57 DM dm;
58 DMTS_DA *dmdats = (DMTS_DA*)ctx;
59 DMDALocalInfo info;
60 Vec Xloc,Xdotloc;
61 void *x,*f,*xdot;
62
63 PetscFunctionBegin;
64 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
65 PetscValidHeaderSpecific(X,VEC_CLASSID,2);
66 PetscValidHeaderSpecific(F,VEC_CLASSID,3);
67 if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
68 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
69 ierr = DMGetLocalVector(dm,&Xdotloc);CHKERRQ(ierr);
70 ierr = DMGlobalToLocalBegin(dm,Xdot,INSERT_VALUES,Xdotloc);CHKERRQ(ierr);
71 ierr = DMGlobalToLocalEnd(dm,Xdot,INSERT_VALUES,Xdotloc);CHKERRQ(ierr);
72 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
73 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
74 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
75 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
76 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
77 ierr = DMDAVecGetArray(dm,Xdotloc,&xdot);CHKERRQ(ierr);
78 switch (dmdats->ifunctionlocalimode) {
79 case INSERT_VALUES: {
80 ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);
81 CHKMEMQ;
82 ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr);
83 CHKMEMQ;
84 ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);
85 } break;
86 case ADD_VALUES: {
87 Vec Floc;
88 ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
89 ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
90 ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);
91 CHKMEMQ;
92 ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr);
93 CHKMEMQ;
94 ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);
95 ierr = VecZeroEntries(F);CHKERRQ(ierr);
96 ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
97 ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
98 ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
99 } break;
100 default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode);
101 }
102 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
103 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
104 ierr = DMDAVecRestoreArray(dm,Xdotloc,&xdot);CHKERRQ(ierr);
105 ierr = DMRestoreLocalVector(dm,&Xdotloc);CHKERRQ(ierr);
106 PetscFunctionReturn(0);
107 }
108
TSComputeIJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void * ctx)109 static PetscErrorCode TSComputeIJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *ctx)
110 {
111 PetscErrorCode ierr;
112 DM dm;
113 DMTS_DA *dmdats = (DMTS_DA*)ctx;
114 DMDALocalInfo info;
115 Vec Xloc;
116 void *x,*xdot;
117
118 PetscFunctionBegin;
119 if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
120 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
121
122 if (dmdats->ijacobianlocal) {
123 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
124 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
125 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
126 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
127 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
128 ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr);
129 CHKMEMQ;
130 ierr = (*dmdats->ijacobianlocal)(&info,ptime,x,xdot,shift,A,B,dmdats->ijacobianlocalctx);CHKERRQ(ierr);
131 CHKMEMQ;
132 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
133 ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr);
134 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
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 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
139 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
140 }
141 PetscFunctionReturn(0);
142 }
143
TSComputeRHSFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec F,void * ctx)144 static PetscErrorCode TSComputeRHSFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec F,void *ctx)
145 {
146 PetscErrorCode ierr;
147 DM dm;
148 DMTS_DA *dmdats = (DMTS_DA*)ctx;
149 DMDALocalInfo info;
150 Vec Xloc;
151 void *x,*f;
152
153 PetscFunctionBegin;
154 PetscValidHeaderSpecific(ts,TS_CLASSID,1);
155 PetscValidHeaderSpecific(X,VEC_CLASSID,2);
156 PetscValidHeaderSpecific(F,VEC_CLASSID,3);
157 if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
158 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
159 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
160 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
161 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
162 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
163 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
164 switch (dmdats->rhsfunctionlocalimode) {
165 case INSERT_VALUES: {
166 ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr);
167 CHKMEMQ;
168 ierr = (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);CHKERRQ(ierr);
169 CHKMEMQ;
170 ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr);
171 } break;
172 case ADD_VALUES: {
173 Vec Floc;
174 ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr);
175 ierr = VecZeroEntries(Floc);CHKERRQ(ierr);
176 ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr);
177 CHKMEMQ;
178 ierr = (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);CHKERRQ(ierr);
179 CHKMEMQ;
180 ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr);
181 ierr = VecZeroEntries(F);CHKERRQ(ierr);
182 ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
183 ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr);
184 ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr);
185 } break;
186 default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->rhsfunctionlocalimode);
187 }
188 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
189 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
190 PetscFunctionReturn(0);
191 }
192
TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat A,Mat B,void * ctx)193 static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat A,Mat B,void *ctx)
194 {
195 PetscErrorCode ierr;
196 DM dm;
197 DMTS_DA *dmdats = (DMTS_DA*)ctx;
198 DMDALocalInfo info;
199 Vec Xloc;
200 void *x;
201
202 PetscFunctionBegin;
203 if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
204 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
205
206 if (dmdats->rhsjacobianlocal) {
207 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
208 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
209 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
210 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
211 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
212 CHKMEMQ;
213 ierr = (*dmdats->rhsjacobianlocal)(&info,ptime,x,A,B,dmdats->rhsjacobianlocalctx);CHKERRQ(ierr);
214 CHKMEMQ;
215 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
216 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
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 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
221 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
222 }
223 PetscFunctionReturn(0);
224 }
225
226
227 /*@C
228 DMDATSSetRHSFunctionLocal - set a local residual evaluation function
229
230 Logically Collective
231
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
237
238 Calling sequence for func:
239
240 $ func(DMDALocalInfo info,PetscReal t,void *x,void *f,void *ctx)
241
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
247
248 Level: beginner
249
250 .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal()
251 @*/
DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void * ctx)252 PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx)
253 {
254 PetscErrorCode ierr;
255 DMTS sdm;
256 DMTS_DA *dmdats;
257
258 PetscFunctionBegin;
259 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
260 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr);
261 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr);
262 dmdats->rhsfunctionlocalimode = imode;
263 dmdats->rhsfunctionlocal = func;
264 dmdats->rhsfunctionlocalctx = ctx;
265 ierr = DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);CHKERRQ(ierr);
266 PetscFunctionReturn(0);
267 }
268
269 /*@C
270 DMDATSSetRHSJacobianLocal - set a local residual evaluation function
271
272 Logically Collective
273
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
278
279 Calling sequence for func:
280
281 $ func(DMDALocalInfo* info,PetscReal t,void* x,Mat J,Mat B,void *ctx);
282
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
289
290 Level: beginner
291
292 .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal()
293 @*/
DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void * ctx)294 PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx)
295 {
296 PetscErrorCode ierr;
297 DMTS sdm;
298 DMTS_DA *dmdats;
299
300 PetscFunctionBegin;
301 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
302 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr);
303 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr);
304 dmdats->rhsjacobianlocal = func;
305 dmdats->rhsjacobianlocalctx = ctx;
306 ierr = DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);CHKERRQ(ierr);
307 PetscFunctionReturn(0);
308 }
309
310
311 /*@C
312 DMDATSSetIFunctionLocal - set a local residual evaluation function
313
314 Logically Collective
315
316 Input Arguments:
317 + dm - DM to associate callback with
318 . func - local residual evaluation
319 - ctx - optional context for local residual evaluation
320
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
328
329 Level: beginner
330
331 .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal()
332 @*/
DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void * ctx)333 PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx)
334 {
335 PetscErrorCode ierr;
336 DMTS sdm;
337 DMTS_DA *dmdats;
338
339 PetscFunctionBegin;
340 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
341 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr);
342 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr);
343 dmdats->ifunctionlocalimode = imode;
344 dmdats->ifunctionlocal = func;
345 dmdats->ifunctionlocalctx = ctx;
346 ierr = DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);CHKERRQ(ierr);
347 PetscFunctionReturn(0);
348 }
349
350 /*@C
351 DMDATSSetIJacobianLocal - set a local residual evaluation function
352
353 Logically Collective
354
355 Input Arguments:
356 + dm - DM to associate callback with
357 . func - local residual evaluation
358 - ctx - optional context for local residual evaluation
359
360 Calling sequence for func:
361
362 $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,PetscScalar shift,Mat J,Mat B,void *ctx);
363
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
372
373 Level: beginner
374
375 .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal()
376 @*/
DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void * ctx)377 PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx)
378 {
379 PetscErrorCode ierr;
380 DMTS sdm;
381 DMTS_DA *dmdats;
382
383 PetscFunctionBegin;
384 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
385 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr);
386 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr);
387 dmdats->ijacobianlocal = func;
388 dmdats->ijacobianlocalctx = ctx;
389 ierr = DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);CHKERRQ(ierr);
390 PetscFunctionReturn(0);
391 }
392
TSMonitorDMDARayDestroy(void ** mctx)393 PetscErrorCode TSMonitorDMDARayDestroy(void **mctx)
394 {
395 TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) *mctx;
396 PetscErrorCode ierr;
397
398 PetscFunctionBegin;
399 if (rayctx->lgctx) {ierr = TSMonitorLGCtxDestroy(&rayctx->lgctx);CHKERRQ(ierr);}
400 ierr = VecDestroy(&rayctx->ray);CHKERRQ(ierr);
401 ierr = VecScatterDestroy(&rayctx->scatter);CHKERRQ(ierr);
402 ierr = PetscViewerDestroy(&rayctx->viewer);CHKERRQ(ierr);
403 ierr = PetscFree(rayctx);CHKERRQ(ierr);
404 PetscFunctionReturn(0);
405 }
406
TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void * mctx)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;
412
413 PetscFunctionBegin;
414 ierr = TSGetSolution(ts,&solution);CHKERRQ(ierr);
415 ierr = VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
416 ierr = VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
417 if (rayctx->viewer) {
418 ierr = VecView(rayctx->ray,rayctx->viewer);CHKERRQ(ierr);
419 }
420 PetscFunctionReturn(0);
421 }
422
TSMonitorLGDMDARay(TS ts,PetscInt step,PetscReal ptime,Vec u,void * ctx)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;
431
432 PetscFunctionBegin;
433 ierr = VecScatterBegin(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
434 ierr = VecScatterEnd(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
435 if (!step) {
436 PetscDrawAxis axis;
437
438 ierr = PetscDrawLGGetAxis(lgctx->lg, &axis);CHKERRQ(ierr);
439 ierr = PetscDrawAxisSetLabels(axis, "Solution Ray as function of time", "Time", "Solution");CHKERRQ(ierr);
440 ierr = VecGetLocalSize(rayctx->ray, &dim);CHKERRQ(ierr);
441 ierr = PetscDrawLGSetDimension(lgctx->lg, dim);CHKERRQ(ierr);
442 ierr = PetscDrawLGReset(lgctx->lg);CHKERRQ(ierr);
443 }
444 ierr = VecGetArrayRead(v, &a);CHKERRQ(ierr);
445 #if defined(PETSC_USE_COMPLEX)
446 {
447 PetscReal *areal;
448 PetscInt i,n;
449 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
450 ierr = PetscMalloc1(n, &areal);CHKERRQ(ierr);
451 for (i = 0; i < n; ++i) areal[i] = PetscRealPart(a[i]);
452 ierr = PetscDrawLGAddCommonPoint(lgctx->lg, ptime, areal);CHKERRQ(ierr);
453 ierr = PetscFree(areal);CHKERRQ(ierr);
454 }
455 #else
456 ierr = PetscDrawLGAddCommonPoint(lgctx->lg, ptime, a);CHKERRQ(ierr);
457 #endif
458 ierr = VecRestoreArrayRead(v, &a);CHKERRQ(ierr);
459 if (((lgctx->howoften > 0) && (!(step % lgctx->howoften))) || ((lgctx->howoften == -1) && ts->reason)) {
460 ierr = PetscDrawLGDraw(lgctx->lg);CHKERRQ(ierr);
461 ierr = PetscDrawLGSave(lgctx->lg);CHKERRQ(ierr);
462 }
463 PetscFunctionReturn(0);
464 }
465