1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscsys.h> 3 4 #ifdef PETSC_HAVE_REVOLVE 5 #include <wrap_revolve.h> 6 #endif 7 8 typedef struct _StackElement { 9 PetscInt stepnum; 10 Vec X; 11 Vec *Y; 12 PetscReal time; 13 PetscReal timeprev; 14 } *StackElement; 15 16 typedef struct _RevolveCTX { 17 PetscBool reverseonestep; 18 PetscInt snaps_in; 19 PetscInt stepsleft; 20 PetscInt check; 21 PetscInt oldcapo; 22 PetscInt capo; 23 PetscInt fine; 24 PetscInt info; 25 } RevolveCTX; 26 27 typedef struct _Stack { 28 PetscBool userevolve; 29 PetscBool recompute; 30 MPI_Comm comm; 31 RevolveCTX *rctx; 32 PetscInt max_cps; /* maximum stack size */ 33 PetscInt stride; 34 PetscInt total_steps; /* total number of steps */ 35 PetscInt numY; 36 PetscInt stacksize; 37 PetscInt top; /* top of the stack */ 38 StackElement *stack; /* container */ 39 } Stack; 40 41 static PetscErrorCode StackCreate(MPI_Comm,Stack *,PetscInt,PetscInt); 42 static PetscErrorCode StackDestroy(Stack*); 43 static PetscErrorCode StackPush(Stack*,StackElement); 44 static PetscErrorCode StackPop(Stack*,StackElement*); 45 static PetscErrorCode StackTop(Stack*,StackElement*); 46 static PetscErrorCode StackDumpAll(TS,Stack*,PetscInt); 47 static PetscErrorCode StackLoadAll(TS,Stack*,PetscInt); 48 49 #ifdef PETSC_HAVE_REVOLVE 50 static void printwhattodo(PetscInt whattodo,RevolveCTX *rctx,PetscInt shift) 51 { 52 switch(whattodo) { 53 case 1: 54 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mAdvance from %D to %D\033[0m\n",rctx->oldcapo+shift,rctx->capo+shift); 55 break; 56 case 2: 57 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mStore in checkpoint number %D\n\033[0m",rctx->check); 58 break; 59 case 3: 60 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mFirst turn: Initialize adjoints and reverse first step\033[0m\n"); 61 break; 62 case 4: 63 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mForward and reverse one step\033[0m\n"); 64 break; 65 case 5: 66 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mRestore in checkpoint number %D\033[0m\n",rctx->check); 67 break; 68 case -1: 69 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mError!"); 70 break; 71 } 72 } 73 #endif 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "StackCreate" 77 static PetscErrorCode StackCreate(MPI_Comm comm,Stack *s,PetscInt size,PetscInt ny) 78 { 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 s->top = -1; 83 s->comm = comm; 84 s->numY = ny; 85 86 ierr = PetscMalloc1(size*sizeof(StackElement),&s->stack);CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "StackDestroy" 92 static PetscErrorCode StackDestroy(Stack *s) 93 { 94 PetscInt i; 95 PetscErrorCode ierr; 96 97 PetscFunctionBegin; 98 if (s->top>-1) { 99 for (i=0;i<=s->top;i++) { 100 ierr = VecDestroy(&s->stack[i]->X);CHKERRQ(ierr); 101 ierr = VecDestroyVecs(s->numY,&s->stack[i]->Y);CHKERRQ(ierr); 102 ierr = PetscFree(s->stack[i]);CHKERRQ(ierr); 103 } 104 } 105 ierr = PetscFree(s->stack);CHKERRQ(ierr); 106 if (s->userevolve) { 107 ierr = PetscFree(s->rctx);CHKERRQ(ierr); 108 } 109 ierr = PetscFree(s);CHKERRQ(ierr); 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "StackPush" 115 static PetscErrorCode StackPush(Stack *s,StackElement e) 116 { 117 PetscFunctionBegin; 118 if (s->top+1 >= s->stacksize) SETERRQ1(s->comm,PETSC_ERR_MEMC,"Maximum stack size (%D) exceeded",s->stacksize); 119 s->stack[++s->top] = e; 120 PetscFunctionReturn(0); 121 } 122 123 #undef __FUNCT__ 124 #define __FUNCT__ "StackPop" 125 static PetscErrorCode StackPop(Stack *s,StackElement *e) 126 { 127 PetscFunctionBegin; 128 if (s->top == -1) SETERRQ(s->comm,PETSC_ERR_MEMC,"Emptry stack"); 129 *e = s->stack[s->top--]; 130 PetscFunctionReturn(0); 131 } 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "StackTop" 135 static PetscErrorCode StackTop(Stack *s,StackElement *e) 136 { 137 PetscFunctionBegin; 138 *e = s->stack[s->top]; 139 PetscFunctionReturn(0); 140 } 141 142 #undef __FUNCT__ 143 #define __FUNCT__ "OutputBIN" 144 static PetscErrorCode OutputBIN(const char *filename,PetscViewer *viewer) 145 { 146 PetscErrorCode ierr; 147 148 PetscFunctionBegin; 149 ierr = PetscViewerCreate(PETSC_COMM_WORLD,viewer);CHKERRQ(ierr); 150 ierr = PetscViewerSetType(*viewer,PETSCVIEWERBINARY);CHKERRQ(ierr); 151 ierr = PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);CHKERRQ(ierr); 152 ierr = PetscViewerFileSetName(*viewer,filename);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 #undef __FUNCT__ 157 #define __FUNCT__ "StackDumpAll" 158 static PetscErrorCode StackDumpAll(TS ts,Stack *s,PetscInt id) 159 { 160 PetscInt i,j; 161 StackElement e; 162 PetscViewer viewer; 163 char filename[PETSC_MAX_PATH_LEN]; 164 Vec *Y; 165 PetscErrorCode ierr; 166 167 PetscFunctionBegin; 168 if (id == 1) { 169 #if defined(PETSC_HAVE_POPEN) 170 PetscMPIInt rank; 171 ierr = MPI_Comm_rank(s->comm,&rank);CHKERRQ(ierr); 172 if (!rank) { 173 char command[PETSC_MAX_PATH_LEN]; 174 FILE *fd; 175 int err; 176 177 ierr = PetscMemzero(command,sizeof(command));CHKERRQ(ierr); 178 ierr = PetscSNPrintf(command,PETSC_MAX_PATH_LEN,"rm -fr %s","SA-data");CHKERRQ(ierr); 179 ierr = PetscPOpen(PETSC_COMM_SELF,NULL,command,"r",&fd);CHKERRQ(ierr); 180 ierr = PetscPClose(PETSC_COMM_SELF,fd,&err);CHKERRQ(ierr); 181 ierr = PetscSNPrintf(command,PETSC_MAX_PATH_LEN,"mkdir %s","SA-data");CHKERRQ(ierr); 182 ierr = PetscPOpen(PETSC_COMM_SELF,NULL,command,"r",&fd);CHKERRQ(ierr); 183 ierr = PetscPClose(PETSC_COMM_SELF,fd,&err);CHKERRQ(ierr); 184 } 185 #endif 186 } 187 ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-STACK%06d.bin",id);CHKERRQ(ierr); 188 ierr = OutputBIN(filename,&viewer);CHKERRQ(ierr); 189 for (i=0;i<s->stacksize;i++) { 190 e = s->stack[i]; 191 ierr = PetscViewerBinaryWrite(viewer,&e->stepnum,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 192 ierr = VecView(e->X,viewer);CHKERRQ(ierr); 193 for (j=0;j<s->numY;j++) { 194 ierr = VecView(e->Y[j],viewer);CHKERRQ(ierr); 195 } 196 ierr = PetscViewerBinaryWrite(viewer,&e->time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 197 ierr = PetscViewerBinaryWrite(viewer,&e->timeprev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 198 } 199 /* dump the state inside TS from the current step */ 200 ierr = PetscViewerBinaryWrite(viewer,&ts->total_steps,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 201 ierr = VecView(ts->vec_sol,viewer);CHKERRQ(ierr); 202 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 203 for (j=0;j<s->numY;j++) { 204 ierr = VecView(Y[j],viewer);CHKERRQ(ierr); 205 } 206 ierr = PetscViewerBinaryWrite(viewer,&ts->ptime,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 207 ierr = PetscViewerBinaryWrite(viewer,&ts->ptime_prev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 208 for (i=0;i<s->stacksize;i++) { 209 ierr = StackPop(s,&e);CHKERRQ(ierr); 210 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 211 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 212 ierr = PetscFree(e);CHKERRQ(ierr); 213 } 214 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "StackLoadAll" 220 static PetscErrorCode StackLoadAll(TS ts,Stack *s,PetscInt id) 221 { 222 Vec *Y; 223 PetscInt i,j; 224 StackElement e; 225 PetscViewer viewer; 226 char filename[PETSC_MAX_PATH_LEN]; 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-STACK%06d.bin",id);CHKERRQ(ierr); 231 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 232 for (i=0;i<s->stacksize;i++) { 233 ierr = PetscCalloc1(1,&e); 234 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 235 ierr = VecDuplicate(Y[0],&e->X);CHKERRQ(ierr); 236 ierr = VecDuplicateVecs(Y[0],s->numY,&e->Y);CHKERRQ(ierr); 237 ierr = StackPush(s,e);CHKERRQ(ierr); 238 } 239 for (i=0;i<s->stacksize;i++) { 240 e = s->stack[i]; 241 ierr = PetscViewerBinaryRead(viewer,&e->stepnum,1,NULL,PETSC_INT);CHKERRQ(ierr); 242 ierr = VecLoad(e->X,viewer);CHKERRQ(ierr); 243 for (j=0;j<s->numY;j++) { 244 ierr = VecLoad(e->Y[j],viewer);CHKERRQ(ierr); 245 } 246 ierr = PetscViewerBinaryRead(viewer,&e->time,1,NULL,PETSC_REAL);CHKERRQ(ierr); 247 ierr = PetscViewerBinaryRead(viewer,&e->timeprev,1,NULL,PETSC_REAL);CHKERRQ(ierr); 248 } 249 /* load the additioinal state into TS */ 250 ierr = PetscViewerBinaryRead(viewer,&ts->total_steps,1,NULL,PETSC_INT);CHKERRQ(ierr); 251 ierr = VecLoad(ts->vec_sol,viewer);CHKERRQ(ierr); 252 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 253 for (j=0;j<s->numY;j++) { 254 ierr = VecLoad(Y[j],viewer);CHKERRQ(ierr); 255 } 256 ierr = PetscViewerBinaryRead(viewer,&ts->ptime,1,NULL,PETSC_REAL);CHKERRQ(ierr); 257 ierr = PetscViewerBinaryRead(viewer,&ts->ptime_prev,1,NULL,PETSC_REAL);CHKERRQ(ierr); 258 ierr = TSSetTimeStep(ts,ts->ptime-ts->ptime_prev);CHKERRQ(ierr); 259 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 263 #undef __FUNCT__ 264 #define __FUNCT__ "TSTrajectorySetStride_Memory" 265 PetscErrorCode TSTrajectorySetStride_Memory(TSTrajectory tj,TS ts,PetscInt stride) 266 { 267 Stack *s = (Stack*)tj->data; 268 269 PetscFunctionBegin; 270 s->stride = stride; 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNCT__ 275 #define __FUNCT__ "TSTrajectorySetMaxCheckpoints_Memory" 276 PetscErrorCode TSTrajectorySetMaxCheckpoints_Memory(TSTrajectory tj,TS ts,PetscInt max_cps) 277 { 278 Stack *s = (Stack*)tj->data; 279 280 PetscFunctionBegin; 281 s->max_cps = max_cps; 282 PetscFunctionReturn(0); 283 } 284 285 #undef __FUNCT__ 286 #define __FUNCT__ "TSTrajectorySetFromOptions_Memory" 287 PetscErrorCode TSTrajectorySetFromOptions_Memory(PetscOptions *PetscOptionsObject,TSTrajectory tj) 288 { 289 Stack *s = (Stack*)tj->data; 290 PetscErrorCode ierr; 291 292 PetscFunctionBegin; 293 ierr = PetscOptionsHead(PetscOptionsObject,"Memory based TS trajectory options");CHKERRQ(ierr); 294 { 295 ierr = PetscOptionsInt("-tstrajectory_max_cps","Maximum number of checkpoints","TSTrajectorySetMaxCheckpoints_Memory",s->max_cps,&s->max_cps,NULL);CHKERRQ(ierr); 296 ierr = PetscOptionsInt("-tstrajectory_stride","Stride to save checkpoints to file","TSTrajectorySetStride_Memory",s->stride,&s->stride,NULL);CHKERRQ(ierr); 297 } 298 ierr = PetscOptionsTail();CHKERRQ(ierr); 299 PetscFunctionReturn(0); 300 } 301 302 #undef __FUNCT__ 303 #define __FUNCT__ "TSTrajectorySetUp_Memory" 304 PetscErrorCode TSTrajectorySetUp_Memory(TSTrajectory tj,TS ts) 305 { 306 Stack *s = (Stack*)tj->data; 307 RevolveCTX *rctx; 308 PetscInt numY; 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 ierr = TSGetStages(ts,&numY,NULL);CHKERRQ(ierr); 313 s->total_steps = PetscMin(ts->max_steps,(PetscInt)(ceil(ts->max_time/ts->time_step))); 314 315 if ((s->stride>1 && s->max_cps>1 && s->max_cps<s->stride-1)||(s->stride<=1 && s->max_cps>1 && s->max_cps<s->total_steps-1)) { 316 #ifdef PETSC_HAVE_REVOLVE 317 s->userevolve = PETSC_TRUE; 318 #else 319 SETERRQ(s->comm,PETSC_ERR_SUP,"revolve is needed when there is not enought memory to checkpoint all time steps according to the user's settings, please reconfigure with the additional option --download-revolve."); 320 #endif 321 s->recompute = PETSC_FALSE; 322 ierr = PetscCalloc1(1,&rctx);CHKERRQ(ierr); 323 rctx->snaps_in = s->max_cps; /* for theta methods snaps_in=2*max_cps */ 324 rctx->reverseonestep = PETSC_FALSE; 325 rctx->check = -1; 326 rctx->oldcapo = 0; 327 rctx->capo = 0; 328 rctx->info = 2; 329 s->rctx = rctx; 330 if (s->stride>1) rctx->fine = s->stride; 331 else rctx->fine = s->total_steps; 332 s->stacksize = s->max_cps; 333 } else { 334 s->userevolve = PETSC_FALSE; 335 if (s->stride>1) s->stacksize = s->stride-1; 336 else s->stacksize = s->total_steps-1; 337 } 338 ierr = StackCreate(PetscObjectComm((PetscObject)ts),s,s->stacksize,numY);CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "TSTrajectorySet_Memory" 344 PetscErrorCode TSTrajectorySet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X) 345 { 346 PetscInt i; 347 Vec *Y; 348 PetscReal timeprev; 349 StackElement e; 350 Stack *s = (Stack*)tj->data; 351 PetscInt localstepnum,id; 352 RevolveCTX *rctx; 353 #ifdef PETSC_HAVE_REVOLVE 354 PetscInt whattodo,rank,shift; 355 #endif 356 PetscErrorCode ierr; 357 358 PetscFunctionBegin; 359 if (!s->recompute) { /* use global stepnum in the forward sweep */ 360 ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr); 361 } 362 if (s->stride>1) { /* multilevel mode */ 363 localstepnum = stepnum%s->stride; 364 if (stepnum!=0 && stepnum!=s->total_steps && localstepnum==0 && !s->recompute) { /* never need to recompute localstepnum=0 */ 365 #ifdef PETSC_HAVE_REVOLVE 366 PetscPrintf(PETSC_COMM_WORLD,"\x1B[33mDump stack to file\033[0m\n"); 367 #endif 368 id = stepnum/s->stride; 369 ierr = StackDumpAll(ts,s,id);CHKERRQ(ierr); 370 s->top = -1; /* reset top */ 371 #ifdef PETSC_HAVE_REVOLVE 372 if (s->userevolve) { 373 rctx = s->rctx; 374 rctx->check = -1; 375 rctx->capo = 0; 376 rctx->fine = s->stride; 377 } 378 wrap_revolve_reset(); 379 #endif 380 } 381 } else { /* in-memory mode */ 382 localstepnum = stepnum; 383 } 384 385 if (s->userevolve) { 386 #ifdef PETSC_HAVE_REVOLVE 387 rctx = s->rctx; 388 if (rctx->reverseonestep && stepnum==s->total_steps) { /* intermediate information is ready inside TS, this happens at last time step */ 389 rctx->reverseonestep = PETSC_FALSE; 390 PetscFunctionReturn(0); 391 } 392 if (rctx->stepsleft!=0) { /* advance the solution without checkpointing anything as Revolve requires */ 393 rctx->stepsleft--; 394 PetscFunctionReturn(0); 395 } 396 397 /* let Revolve determine what to do next */ 398 shift = stepnum-localstepnum; 399 rctx->capo = localstepnum; 400 rctx->oldcapo = rctx->capo; 401 whattodo = wrap_revolve(&rctx->check,&rctx->capo,&rctx->fine,&rctx->snaps_in,&rctx->info,&rank); 402 printwhattodo(whattodo,rctx,shift); 403 if (whattodo==-1) SETERRQ(s->comm,PETSC_ERR_LIB,"Error in the Revolve library"); 404 if (whattodo==1) { /* advance some time steps */ 405 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 406 PetscFunctionReturn(0); 407 } 408 if (whattodo==3 || whattodo==4) { /* ready for a reverse step */ 409 rctx->reverseonestep = PETSC_TRUE; 410 PetscFunctionReturn(0); 411 } 412 if (whattodo==5) { /* restore a checkpoint and ask Revolve what to do next */ 413 rctx->oldcapo = rctx->capo; 414 whattodo = wrap_revolve(&rctx->check,&rctx->capo,&rctx->fine,&rctx->snaps_in,&rctx->info,&rank); /* must return 1*/ 415 printwhattodo(whattodo,rctx,shift); 416 rctx->stepsleft = rctx->capo-rctx->oldcapo; 417 PetscFunctionReturn(0); 418 } 419 if (whattodo==2) { /* store a checkpoint and ask Revolve how many time steps to advance next */ 420 rctx->oldcapo = rctx->capo; 421 whattodo = wrap_revolve(&rctx->check,&rctx->capo,&rctx->fine,&rctx->snaps_in,&rctx->info,&rank); /* must return 1*/ 422 printwhattodo(whattodo,rctx,shift); 423 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 424 } 425 #endif 426 } else { /* Revolve is not used */ 427 /* skip the first and the last steps of each stride or the whole interval */ 428 if (localstepnum==0 || stepnum==s->total_steps) PetscFunctionReturn(0); 429 } 430 431 /* checkpoint to memmory */ 432 if (localstepnum==s->top) { /* overwrite the top checkpoint, this might happen when the time interval is split into several smaller ones, each corresponding to a call of TSSolve() */ 433 ierr = StackTop(s,&e);CHKERRQ(ierr); 434 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 435 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 436 for (i=0;i<s->numY;i++) { 437 ierr = VecCopy(Y[i],e->Y[i]);CHKERRQ(ierr); 438 } 439 e->stepnum = stepnum; 440 e->time = time; 441 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 442 e->timeprev = timeprev; 443 } else { 444 if (localstepnum<s->top) SETERRQ(s->comm,PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 445 ierr = PetscCalloc1(1,&e);CHKERRQ(ierr); 446 ierr = VecDuplicate(X,&e->X);CHKERRQ(ierr); 447 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 448 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 449 ierr = VecDuplicateVecs(Y[0],s->numY,&e->Y);CHKERRQ(ierr); 450 for (i=0;i<s->numY;i++) { 451 ierr = VecCopy(Y[i],e->Y[i]);CHKERRQ(ierr); 452 } 453 e->stepnum = stepnum; 454 e->time = time; 455 /* for consistency */ 456 if (stepnum == 0) { 457 e->timeprev = e->time - ts->time_step; 458 } else { 459 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 460 e->timeprev = timeprev; 461 } 462 ierr = StackPush(s,e);CHKERRQ(ierr); 463 } 464 PetscFunctionReturn(0); 465 } 466 467 #undef __FUNCT__ 468 #define __FUNCT__ "TSTrajectoryGet_Memory" 469 PetscErrorCode TSTrajectoryGet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t) 470 { 471 Vec *Y; 472 PetscInt i; 473 StackElement e; 474 Stack *s = (Stack*)tj->data; 475 PetscReal stepsize; 476 PetscInt localstepnum,id; 477 #ifdef PETSC_HAVE_REVOLVE 478 PetscInt whattodo,rank,shift; 479 #endif 480 PetscErrorCode ierr; 481 482 PetscFunctionBegin; 483 ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr); 484 if (s->stride>1) { /* multilevel mode */ 485 localstepnum = stepnum%s->stride; 486 if (localstepnum==0 && stepnum!=0 && stepnum!=s->total_steps) { 487 #ifdef PETSC_HAVE_REVOLVE 488 PetscPrintf(PETSC_COMM_WORLD,"\x1B[33mLoad stack from file\033[0m\n"); 489 #endif 490 id = stepnum/s->stride; 491 ierr = StackLoadAll(ts,s,id);CHKERRQ(ierr); 492 s->top = s->stacksize-1; 493 #ifdef PETSC_HAVE_REVOLVE 494 if (s->userevolve) { 495 wrap_revolve_reset(); 496 s->rctx->check = -1; 497 s->rctx->capo = 0; 498 s->rctx->fine = s->stride; 499 whattodo = 0; 500 while(whattodo!=3) { /* stupid revolve */ 501 whattodo = wrap_revolve(&s->rctx->check,&s->rctx->capo,&s->rctx->fine,&s->rctx->snaps_in,&s->rctx->info,&rank); 502 } 503 } 504 #endif 505 } 506 } else { 507 localstepnum = stepnum; 508 } 509 #ifdef PETSC_HAVE_REVOLVE 510 if (s->userevolve && s->rctx->reverseonestep) { /* ready for the reverse step */ 511 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 512 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); 513 s->rctx->reverseonestep = PETSC_FALSE; 514 PetscFunctionReturn(0); 515 } 516 #endif 517 if (stepnum==s->total_steps || localstepnum==0) { 518 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 519 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); 520 PetscFunctionReturn(0); 521 } 522 /* restore a checkpoint */ 523 ierr = StackTop(s,&e);CHKERRQ(ierr); 524 ierr = VecCopy(e->X,ts->vec_sol);CHKERRQ(ierr); 525 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 526 for (i=0;i<s->numY;i++) { 527 ierr = VecCopy(e->Y[i],Y[i]);CHKERRQ(ierr); 528 } 529 *t = e->time; 530 531 if (e->stepnum < stepnum) { /* need recomputation */ 532 s->recompute = PETSC_TRUE; 533 shift = stepnum-localstepnum; 534 #ifdef PETSC_HAVE_REVOLVE 535 s->rctx->capo = localstepnum; 536 whattodo = wrap_revolve(&s->rctx->check,&s->rctx->capo,&s->rctx->fine,&s->rctx->snaps_in,&s->rctx->info,&rank); 537 printwhattodo(whattodo,s->rctx,shift); 538 #endif 539 ierr = TSSetTimeStep(ts,(*t)-e->timeprev);CHKERRQ(ierr); 540 /* reset ts context */ 541 PetscInt steps = ts->steps; 542 ts->steps = e->stepnum; 543 ts->ptime = e->time; 544 ts->ptime_prev = e->timeprev; 545 for (i=e->stepnum;i<stepnum;i++) { /* assume fixed step size */ 546 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 547 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 548 ierr = TSStep(ts);CHKERRQ(ierr); 549 if (ts->event) { 550 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 551 } 552 if (!ts->steprollback) { 553 ierr = TSPostStep(ts);CHKERRQ(ierr); 554 } 555 } 556 /* reverseonestep must be true after the for loop */ 557 ts->steps = steps; 558 ts->total_steps = stepnum; 559 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 560 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); 561 if (stepnum-e->stepnum==1) { /* the restored checkpoint can be deleted now */ 562 ierr = StackPop(s,&e);CHKERRQ(ierr); 563 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 564 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 565 ierr = PetscFree(e);CHKERRQ(ierr); 566 } 567 #ifdef PETSC_HAVE_REVOLVE 568 s->rctx->reverseonestep = PETSC_FALSE; 569 #endif 570 } else if (e->stepnum == stepnum) { /* restore the data directly from checkpoints */ 571 ierr = TSSetTimeStep(ts,-(*t)+e->timeprev);CHKERRQ(ierr); 572 ierr = StackPop(s,&e);CHKERRQ(ierr); 573 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 574 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 575 ierr = PetscFree(e);CHKERRQ(ierr); 576 } else { 577 SETERRQ2(s->comm,PETSC_ERR_ARG_OUTOFRANGE,"The current step no. is %D, but the step number at top of the stack is %D",stepnum,e->stepnum); 578 } 579 580 PetscFunctionReturn(0); 581 } 582 583 #undef __FUNCT__ 584 #define __FUNCT__ "TSTrajectoryDestroy_Memory" 585 PETSC_EXTERN PetscErrorCode TSTrajectoryDestroy_Memory(TSTrajectory tj) 586 { 587 Stack *s = (Stack*)tj->data; 588 PetscErrorCode ierr; 589 590 PetscFunctionBegin; 591 ierr = StackDestroy(s);CHKERRQ(ierr); 592 PetscFunctionReturn(0); 593 } 594 595 /*MC 596 TSTRAJECTORYMEMORY - Stores each solution of the ODE/ADE in memory 597 598 Level: intermediate 599 600 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType() 601 602 M*/ 603 #undef __FUNCT__ 604 #define __FUNCT__ "TSTrajectoryCreate_Memory" 605 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory tj,TS ts) 606 { 607 Stack *s; 608 PetscErrorCode ierr; 609 610 PetscFunctionBegin; 611 tj->ops->set = TSTrajectorySet_Memory; 612 tj->ops->get = TSTrajectoryGet_Memory; 613 tj->ops->setup = TSTrajectorySetUp_Memory; 614 tj->ops->destroy = TSTrajectoryDestroy_Memory; 615 tj->ops->setfromoptions = TSTrajectorySetFromOptions_Memory; 616 617 ierr = PetscCalloc1(1,&s);CHKERRQ(ierr); 618 s->max_cps = -1; /* -1 indicates that it is not set */ 619 s->stride = 0; /* if not zero, two-level checkpointing will be used */ 620 tj->data = s; 621 PetscFunctionReturn(0); 622 } 623