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 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 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 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 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,3); 66 PetscValidHeaderSpecific(F,VEC_CLASSID,5); 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 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 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,3); 156 PetscValidHeaderSpecific(F,VEC_CLASSID,4); 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 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 /*@C 227 DMDATSSetRHSFunctionLocal - set a local residual evaluation function 228 229 Logically Collective 230 231 Input Arguments: 232 + dm - DM to associate callback with 233 . imode - insert mode for the residual 234 . func - local residual evaluation 235 - ctx - optional context for local residual evaluation 236 237 Calling sequence for func: 238 239 $ func(DMDALocalInfo info,PetscReal t,void *x,void *f,void *ctx) 240 241 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 242 . t - time at which to evaluate residual 243 . x - array of local state information 244 . f - output array of local residual information 245 - ctx - optional user context 246 247 Level: beginner 248 249 .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal() 250 @*/ 251 PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx) 252 { 253 PetscErrorCode ierr; 254 DMTS sdm; 255 DMTS_DA *dmdats; 256 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 259 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 260 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 261 dmdats->rhsfunctionlocalimode = imode; 262 dmdats->rhsfunctionlocal = func; 263 dmdats->rhsfunctionlocalctx = ctx; 264 ierr = DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 268 /*@C 269 DMDATSSetRHSJacobianLocal - set a local residual evaluation function 270 271 Logically Collective 272 273 Input Arguments: 274 + dm - DM to associate callback with 275 . func - local RHS Jacobian evaluation routine 276 - ctx - optional context for local jacobian evaluation 277 278 Calling sequence for func: 279 280 $ func(DMDALocalInfo* info,PetscReal t,void* x,Mat J,Mat B,void *ctx); 281 282 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 283 . t - time at which to evaluate residual 284 . x - array of local state information 285 . J - Jacobian matrix 286 . B - preconditioner matrix; often same as J 287 - ctx - optional context passed above 288 289 Level: beginner 290 291 .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal() 292 @*/ 293 PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx) 294 { 295 PetscErrorCode ierr; 296 DMTS sdm; 297 DMTS_DA *dmdats; 298 299 PetscFunctionBegin; 300 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 301 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 302 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 303 dmdats->rhsjacobianlocal = func; 304 dmdats->rhsjacobianlocalctx = ctx; 305 ierr = DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);CHKERRQ(ierr); 306 PetscFunctionReturn(0); 307 } 308 309 /*@C 310 DMDATSSetIFunctionLocal - set a local residual evaluation function 311 312 Logically Collective 313 314 Input Arguments: 315 + dm - DM to associate callback with 316 . func - local residual evaluation 317 - ctx - optional context for local residual evaluation 318 319 Calling sequence for func: 320 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 321 . t - time at which to evaluate residual 322 . x - array of local state information 323 . xdot - array of local time derivative information 324 . f - output array of local function evaluation information 325 - ctx - optional context passed above 326 327 Level: beginner 328 329 .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal() 330 @*/ 331 PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx) 332 { 333 PetscErrorCode ierr; 334 DMTS sdm; 335 DMTS_DA *dmdats; 336 337 PetscFunctionBegin; 338 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 339 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 340 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 341 dmdats->ifunctionlocalimode = imode; 342 dmdats->ifunctionlocal = func; 343 dmdats->ifunctionlocalctx = ctx; 344 ierr = DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);CHKERRQ(ierr); 345 PetscFunctionReturn(0); 346 } 347 348 /*@C 349 DMDATSSetIJacobianLocal - set a local residual evaluation function 350 351 Logically Collective 352 353 Input Arguments: 354 + dm - DM to associate callback with 355 . func - local residual evaluation 356 - ctx - optional context for local residual evaluation 357 358 Calling sequence for func: 359 360 $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,PetscScalar shift,Mat J,Mat B,void *ctx); 361 362 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 363 . t - time at which to evaluate the jacobian 364 . x - array of local state information 365 . xdot - time derivative at this state 366 . shift - see TSSetIJacobian() for the meaning of this parameter 367 . J - Jacobian matrix 368 . B - preconditioner matrix; often same as J 369 - ctx - optional context passed above 370 371 Level: beginner 372 373 .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal() 374 @*/ 375 PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx) 376 { 377 PetscErrorCode ierr; 378 DMTS sdm; 379 DMTS_DA *dmdats; 380 381 PetscFunctionBegin; 382 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 383 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 384 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 385 dmdats->ijacobianlocal = func; 386 dmdats->ijacobianlocalctx = ctx; 387 ierr = DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);CHKERRQ(ierr); 388 PetscFunctionReturn(0); 389 } 390 391 PetscErrorCode TSMonitorDMDARayDestroy(void **mctx) 392 { 393 TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) *mctx; 394 PetscErrorCode ierr; 395 396 PetscFunctionBegin; 397 if (rayctx->lgctx) {ierr = TSMonitorLGCtxDestroy(&rayctx->lgctx);CHKERRQ(ierr);} 398 ierr = VecDestroy(&rayctx->ray);CHKERRQ(ierr); 399 ierr = VecScatterDestroy(&rayctx->scatter);CHKERRQ(ierr); 400 ierr = PetscViewerDestroy(&rayctx->viewer);CHKERRQ(ierr); 401 ierr = PetscFree(rayctx);CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 405 PetscErrorCode TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 406 { 407 TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)mctx; 408 Vec solution; 409 PetscErrorCode ierr; 410 411 PetscFunctionBegin; 412 ierr = TSGetSolution(ts,&solution);CHKERRQ(ierr); 413 ierr = VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 414 ierr = VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 415 if (rayctx->viewer) { 416 ierr = VecView(rayctx->ray,rayctx->viewer);CHKERRQ(ierr); 417 } 418 PetscFunctionReturn(0); 419 } 420 421 PetscErrorCode TSMonitorLGDMDARay(TS ts, PetscInt step, PetscReal ptime, Vec u, void *ctx) 422 { 423 TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) ctx; 424 TSMonitorLGCtx lgctx = (TSMonitorLGCtx) rayctx->lgctx; 425 Vec v = rayctx->ray; 426 const PetscScalar *a; 427 PetscInt dim; 428 PetscErrorCode ierr; 429 430 PetscFunctionBegin; 431 ierr = VecScatterBegin(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 432 ierr = VecScatterEnd(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 433 if (!step) { 434 PetscDrawAxis axis; 435 436 ierr = PetscDrawLGGetAxis(lgctx->lg, &axis);CHKERRQ(ierr); 437 ierr = PetscDrawAxisSetLabels(axis, "Solution Ray as function of time", "Time", "Solution");CHKERRQ(ierr); 438 ierr = VecGetLocalSize(rayctx->ray, &dim);CHKERRQ(ierr); 439 ierr = PetscDrawLGSetDimension(lgctx->lg, dim);CHKERRQ(ierr); 440 ierr = PetscDrawLGReset(lgctx->lg);CHKERRQ(ierr); 441 } 442 ierr = VecGetArrayRead(v, &a);CHKERRQ(ierr); 443 #if defined(PETSC_USE_COMPLEX) 444 { 445 PetscReal *areal; 446 PetscInt i,n; 447 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 448 ierr = PetscMalloc1(n, &areal);CHKERRQ(ierr); 449 for (i = 0; i < n; ++i) areal[i] = PetscRealPart(a[i]); 450 ierr = PetscDrawLGAddCommonPoint(lgctx->lg, ptime, areal);CHKERRQ(ierr); 451 ierr = PetscFree(areal);CHKERRQ(ierr); 452 } 453 #else 454 ierr = PetscDrawLGAddCommonPoint(lgctx->lg, ptime, a);CHKERRQ(ierr); 455 #endif 456 ierr = VecRestoreArrayRead(v, &a);CHKERRQ(ierr); 457 if (((lgctx->howoften > 0) && (!(step % lgctx->howoften))) || ((lgctx->howoften == -1) && ts->reason)) { 458 ierr = PetscDrawLGDraw(lgctx->lg);CHKERRQ(ierr); 459 ierr = PetscDrawLGSave(lgctx->lg);CHKERRQ(ierr); 460 } 461 PetscFunctionReturn(0); 462 } 463