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