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