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