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