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; 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,&Xloc);CHKERRQ(ierr); 70 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 71 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 72 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 73 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 74 ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr); 75 switch (dmdats->ifunctionlocalimode) { 76 case INSERT_VALUES: { 77 ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr); 78 CHKMEMQ; 79 ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr); 80 CHKMEMQ; 81 ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr); 82 } break; 83 case ADD_VALUES: { 84 Vec Floc; 85 ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr); 86 ierr = VecZeroEntries(Floc);CHKERRQ(ierr); 87 ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr); 88 CHKMEMQ; 89 ierr = (*dmdats->ifunctionlocal)(&info,ptime,x,xdot,f,dmdats->ifunctionlocalctx);CHKERRQ(ierr); 90 CHKMEMQ; 91 ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr); 92 ierr = VecZeroEntries(F);CHKERRQ(ierr); 93 ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 94 ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 95 ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr); 96 } break; 97 default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->ifunctionlocalimode); 98 } 99 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 100 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 101 ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr); 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode TSComputeIJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,void *ctx) 106 { 107 PetscErrorCode ierr; 108 DM dm; 109 DMTS_DA *dmdats = (DMTS_DA*)ctx; 110 DMDALocalInfo info; 111 Vec Xloc; 112 void *x,*xdot; 113 114 PetscFunctionBegin; 115 if (!dmdats->ifunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context"); 116 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 117 118 if (dmdats->ijacobianlocal) { 119 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 120 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 121 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 122 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 123 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 124 ierr = DMDAVecGetArray(dm,Xdot,&xdot);CHKERRQ(ierr); 125 CHKMEMQ; 126 ierr = (*dmdats->ijacobianlocal)(&info,ptime,x,xdot,shift,A,B,dmdats->ijacobianlocalctx);CHKERRQ(ierr); 127 CHKMEMQ; 128 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 129 ierr = DMDAVecRestoreArray(dm,Xdot,&xdot);CHKERRQ(ierr); 130 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 131 } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeIJacobian_DMDA() called without calling DMDATSSetIJacobian()"); 132 /* This will be redundant if the user called both, but it's too common to forget. */ 133 if (A != B) { 134 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 135 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 136 } 137 PetscFunctionReturn(0); 138 } 139 140 static PetscErrorCode TSComputeRHSFunction_DMDA(TS ts,PetscReal ptime,Vec X,Vec F,void *ctx) 141 { 142 PetscErrorCode ierr; 143 DM dm; 144 DMTS_DA *dmdats = (DMTS_DA*)ctx; 145 DMDALocalInfo info; 146 Vec Xloc; 147 void *x,*f; 148 149 PetscFunctionBegin; 150 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 151 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 152 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 153 if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context"); 154 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 155 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 156 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 157 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 158 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 159 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 160 switch (dmdats->rhsfunctionlocalimode) { 161 case INSERT_VALUES: { 162 ierr = DMDAVecGetArray(dm,F,&f);CHKERRQ(ierr); 163 CHKMEMQ; 164 ierr = (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);CHKERRQ(ierr); 165 CHKMEMQ; 166 ierr = DMDAVecRestoreArray(dm,F,&f);CHKERRQ(ierr); 167 } break; 168 case ADD_VALUES: { 169 Vec Floc; 170 ierr = DMGetLocalVector(dm,&Floc);CHKERRQ(ierr); 171 ierr = VecZeroEntries(Floc);CHKERRQ(ierr); 172 ierr = DMDAVecGetArray(dm,Floc,&f);CHKERRQ(ierr); 173 CHKMEMQ; 174 ierr = (*dmdats->rhsfunctionlocal)(&info,ptime,x,f,dmdats->rhsfunctionlocalctx);CHKERRQ(ierr); 175 CHKMEMQ; 176 ierr = DMDAVecRestoreArray(dm,Floc,&f);CHKERRQ(ierr); 177 ierr = VecZeroEntries(F);CHKERRQ(ierr); 178 ierr = DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 179 ierr = DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F);CHKERRQ(ierr); 180 ierr = DMRestoreLocalVector(dm,&Floc);CHKERRQ(ierr); 181 } break; 182 default: SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Cannot use imode=%d",(int)dmdats->rhsfunctionlocalimode); 183 } 184 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 185 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat A,Mat B,void *ctx) 190 { 191 PetscErrorCode ierr; 192 DM dm; 193 DMTS_DA *dmdats = (DMTS_DA*)ctx; 194 DMDALocalInfo info; 195 Vec Xloc; 196 void *x; 197 198 PetscFunctionBegin; 199 if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context"); 200 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 201 202 if (dmdats->rhsjacobianlocal) { 203 ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr); 204 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 205 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 206 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 207 ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr); 208 CHKMEMQ; 209 ierr = (*dmdats->rhsjacobianlocal)(&info,ptime,x,A,B,dmdats->rhsjacobianlocalctx);CHKERRQ(ierr); 210 CHKMEMQ; 211 ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr); 212 ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr); 213 } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeRHSJacobian_DMDA() called without calling DMDATSSetRHSJacobian()"); 214 /* This will be redundant if the user called both, but it's too common to forget. */ 215 if (A != B) { 216 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 217 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 218 } 219 PetscFunctionReturn(0); 220 } 221 222 223 /*@C 224 DMDATSSetRHSFunctionLocal - set a local residual evaluation function 225 226 Logically Collective 227 228 Input Arguments: 229 + dm - DM to associate callback with 230 . imode - insert mode for the residual 231 . func - local residual evaluation 232 - ctx - optional context for local residual evaluation 233 234 Calling sequence for func: 235 236 $ func(DMDALocalInfo info,PetscReal t,void *x,void *f,void *ctx) 237 238 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 239 . t - time at which to evaluate residual 240 . x - array of local state information 241 . f - output array of local residual information 242 - ctx - optional user context 243 244 Level: beginner 245 246 .seealso: DMTSSetRHSFunction(), DMDATSSetRHSJacobianLocal(), DMDASNESSetFunctionLocal() 247 @*/ 248 PetscErrorCode DMDATSSetRHSFunctionLocal(DM dm,InsertMode imode,DMDATSRHSFunctionLocal func,void *ctx) 249 { 250 PetscErrorCode ierr; 251 DMTS sdm; 252 DMTS_DA *dmdats; 253 254 PetscFunctionBegin; 255 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 256 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 257 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 258 dmdats->rhsfunctionlocalimode = imode; 259 dmdats->rhsfunctionlocal = func; 260 dmdats->rhsfunctionlocalctx = ctx; 261 ierr = DMTSSetRHSFunction(dm,TSComputeRHSFunction_DMDA,dmdats);CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 /*@C 266 DMDATSSetRHSJacobianLocal - set a local residual evaluation function 267 268 Logically Collective 269 270 Input Arguments: 271 + dm - DM to associate callback with 272 . func - local RHS Jacobian evaluation routine 273 - ctx - optional context for local jacobian evaluation 274 275 Calling sequence for func: 276 277 $ func(DMDALocalInfo* info,PetscReal t,void* x,Mat J,Mat B,void *ctx); 278 279 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 280 . t - time at which to evaluate residual 281 . x - array of local state information 282 . J - Jacobian matrix 283 . B - preconditioner matrix; often same as J 284 - ctx - optional context passed above 285 286 Level: beginner 287 288 .seealso: DMTSSetRHSJacobian(), DMDATSSetRHSFunctionLocal(), DMDASNESSetJacobianLocal() 289 @*/ 290 PetscErrorCode DMDATSSetRHSJacobianLocal(DM dm,DMDATSRHSJacobianLocal func,void *ctx) 291 { 292 PetscErrorCode ierr; 293 DMTS sdm; 294 DMTS_DA *dmdats; 295 296 PetscFunctionBegin; 297 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 298 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 299 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 300 dmdats->rhsjacobianlocal = func; 301 dmdats->rhsjacobianlocalctx = ctx; 302 ierr = DMTSSetRHSJacobian(dm,TSComputeRHSJacobian_DMDA,dmdats);CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 307 /*@C 308 DMDATSSetIFunctionLocal - set a local residual evaluation function 309 310 Logically Collective 311 312 Input Arguments: 313 + dm - DM to associate callback with 314 . func - local residual evaluation 315 - ctx - optional context for local residual evaluation 316 317 Calling sequence for func: 318 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 319 . t - time at which to evaluate residual 320 . x - array of local state information 321 . xdot - array of local time derivative information 322 . f - output array of local function evaluation information 323 - ctx - optional context passed above 324 325 Level: beginner 326 327 .seealso: DMTSSetIFunction(), DMDATSSetIJacobianLocal(), DMDASNESSetFunctionLocal() 328 @*/ 329 PetscErrorCode DMDATSSetIFunctionLocal(DM dm,InsertMode imode,DMDATSIFunctionLocal func,void *ctx) 330 { 331 PetscErrorCode ierr; 332 DMTS sdm; 333 DMTS_DA *dmdats; 334 335 PetscFunctionBegin; 336 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 337 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 338 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 339 dmdats->ifunctionlocalimode = imode; 340 dmdats->ifunctionlocal = func; 341 dmdats->ifunctionlocalctx = ctx; 342 ierr = DMTSSetIFunction(dm,TSComputeIFunction_DMDA,dmdats);CHKERRQ(ierr); 343 PetscFunctionReturn(0); 344 } 345 346 /*@C 347 DMDATSSetIJacobianLocal - set a local residual evaluation function 348 349 Logically Collective 350 351 Input Arguments: 352 + dm - DM to associate callback with 353 . func - local residual evaluation 354 - ctx - optional context for local residual evaluation 355 356 Calling sequence for func: 357 358 $ func(DMDALocalInfo* info,PetscReal t,void* x,void *xdot,PetscScalar shift,Mat J,Mat B,void *ctx); 359 360 + info - DMDALocalInfo defining the subdomain to evaluate the residual on 361 . t - time at which to evaluate the jacobian 362 . x - array of local state information 363 . xdot - time derivative at this state 364 . shift - see TSSetIJacobian() for the meaning of this parameter 365 . J - Jacobian matrix 366 . B - preconditioner matrix; often same as J 367 - ctx - optional context passed above 368 369 Level: beginner 370 371 .seealso: DMTSSetJacobian(), DMDATSSetIFunctionLocal(), DMDASNESSetJacobianLocal() 372 @*/ 373 PetscErrorCode DMDATSSetIJacobianLocal(DM dm,DMDATSIJacobianLocal func,void *ctx) 374 { 375 PetscErrorCode ierr; 376 DMTS sdm; 377 DMTS_DA *dmdats; 378 379 PetscFunctionBegin; 380 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 381 ierr = DMGetDMTSWrite(dm,&sdm);CHKERRQ(ierr); 382 ierr = DMDATSGetContext(dm,sdm,&dmdats);CHKERRQ(ierr); 383 dmdats->ijacobianlocal = func; 384 dmdats->ijacobianlocalctx = ctx; 385 ierr = DMTSSetIJacobian(dm,TSComputeIJacobian_DMDA,dmdats);CHKERRQ(ierr); 386 PetscFunctionReturn(0); 387 } 388 389 PetscErrorCode TSMonitorDMDARayDestroy(void **mctx) 390 { 391 TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) *mctx; 392 PetscErrorCode ierr; 393 394 PetscFunctionBegin; 395 if (rayctx->lgctx) {ierr = TSMonitorLGCtxDestroy(&rayctx->lgctx);CHKERRQ(ierr);} 396 ierr = VecDestroy(&rayctx->ray);CHKERRQ(ierr); 397 ierr = VecScatterDestroy(&rayctx->scatter);CHKERRQ(ierr); 398 ierr = PetscViewerDestroy(&rayctx->viewer);CHKERRQ(ierr); 399 ierr = PetscFree(rayctx);CHKERRQ(ierr); 400 PetscFunctionReturn(0); 401 } 402 403 PetscErrorCode TSMonitorDMDARay(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 404 { 405 TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx*)mctx; 406 Vec solution; 407 PetscErrorCode ierr; 408 409 PetscFunctionBegin; 410 ierr = TSGetSolution(ts,&solution);CHKERRQ(ierr); 411 ierr = VecScatterBegin(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 412 ierr = VecScatterEnd(rayctx->scatter,solution,rayctx->ray,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 413 if (rayctx->viewer) { 414 ierr = VecView(rayctx->ray,rayctx->viewer);CHKERRQ(ierr); 415 } 416 PetscFunctionReturn(0); 417 } 418 419 PetscErrorCode TSMonitorLGDMDARay(TS ts, PetscInt step, PetscReal ptime, Vec u, void *ctx) 420 { 421 TSMonitorDMDARayCtx *rayctx = (TSMonitorDMDARayCtx *) ctx; 422 TSMonitorLGCtx lgctx = (TSMonitorLGCtx) rayctx->lgctx; 423 Vec v = rayctx->ray; 424 const PetscScalar *a; 425 PetscInt dim; 426 PetscErrorCode ierr; 427 428 PetscFunctionBegin; 429 ierr = VecScatterBegin(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 430 ierr = VecScatterEnd(rayctx->scatter, u, v, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 431 if (!step) { 432 PetscDrawAxis axis; 433 434 ierr = PetscDrawLGGetAxis(lgctx->lg, &axis);CHKERRQ(ierr); 435 ierr = PetscDrawAxisSetLabels(axis, "Solution Ray as function of time", "Time", "Solution");CHKERRQ(ierr); 436 ierr = VecGetLocalSize(rayctx->ray, &dim);CHKERRQ(ierr); 437 ierr = PetscDrawLGSetDimension(lgctx->lg, dim);CHKERRQ(ierr); 438 ierr = PetscDrawLGReset(lgctx->lg);CHKERRQ(ierr); 439 } 440 ierr = VecGetArrayRead(v, &a);CHKERRQ(ierr); 441 #if defined(PETSC_USE_COMPLEX) 442 { 443 PetscReal *areal; 444 PetscInt i,n; 445 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr); 446 ierr = PetscMalloc1(n, &areal);CHKERRQ(ierr); 447 for (i = 0; i < n; ++i) areal[i] = PetscRealPart(a[i]); 448 ierr = PetscDrawLGAddCommonPoint(lgctx->lg, ptime, areal);CHKERRQ(ierr); 449 ierr = PetscFree(areal);CHKERRQ(ierr); 450 } 451 #else 452 ierr = PetscDrawLGAddCommonPoint(lgctx->lg, ptime, a);CHKERRQ(ierr); 453 #endif 454 ierr = VecRestoreArrayRead(v, &a);CHKERRQ(ierr); 455 if (((lgctx->howoften > 0) && (!(step % lgctx->howoften))) || ((lgctx->howoften == -1) && ts->reason)) { 456 ierr = PetscDrawLGDraw(lgctx->lg);CHKERRQ(ierr); 457 ierr = PetscDrawLGSave(lgctx->lg);CHKERRQ(ierr); 458 } 459 PetscFunctionReturn(0); 460 } 461