1 #define PRINTWHATTODO 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; /* The top of the stack */ 30 PetscInt max_cps; /* The maximum stack size */ 31 PetscInt numY; 32 PetscInt stride; 33 PetscInt max_steps; /* Max number of steps */ 34 MPI_Comm comm; 35 StackElement *stack; /* The storage */ 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 PRINTWHATTODO 47 static void printwhattodo(PetscInt whattodo,RevolveCTX *rctx) 48 { 49 switch(whattodo) { 50 case 1: 51 PetscPrintf(PETSC_COMM_WORLD,"Advance from %D to %D.\n",rctx->oldcapo,rctx->capo); 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 221 for (i=0;i<s->stride;i++) { 222 ierr = PetscCalloc1(1,&e); 223 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 224 ierr = VecDuplicate(Y[0],&e->X);CHKERRQ(ierr); 225 ierr = VecDuplicateVecs(Y[0],s->numY,&e->Y);CHKERRQ(ierr); 226 ierr = StackPush(s,e);CHKERRQ(ierr); 227 } 228 for (i=0;i<s->stride;i++) { 229 e = s->stack[i]; 230 ierr = PetscViewerBinaryRead(viewer,&e->stepnum,1,NULL,PETSC_INT);CHKERRQ(ierr); 231 ierr = VecLoad(e->X,viewer);CHKERRQ(ierr); 232 for (j=0;j<s->numY;j++) { 233 ierr = VecLoad(e->Y[j],viewer);CHKERRQ(ierr); 234 } 235 ierr = PetscViewerBinaryRead(viewer,&e->time,1,NULL,PETSC_REAL);CHKERRQ(ierr); 236 ierr = PetscViewerBinaryRead(viewer,&e->timeprev,1,NULL,PETSC_REAL);CHKERRQ(ierr); 237 } 238 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 239 PetscFunctionReturn(0); 240 } 241 242 #undef __FUNCT__ 243 #define __FUNCT__ "TSTrajectorySetStride_Memory" 244 PetscErrorCode TSTrajectorySetStride_Memory(TSTrajectory tj,TS ts,PetscInt stride) 245 { 246 Stack *s = (Stack*)tj->data; 247 248 PetscFunctionBegin; 249 s->stride = stride; 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "TSTrajectorySetMaxCheckpoints_Memory" 255 PetscErrorCode TSTrajectorySetMaxCheckpoints_Memory(TSTrajectory tj,TS ts,PetscInt max_cps) 256 { 257 Stack *s = (Stack*)tj->data; 258 259 PetscFunctionBegin; 260 s->max_cps = max_cps; 261 PetscFunctionReturn(0); 262 } 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "TSTrajectorySetFromOptions_Memory" 266 PetscErrorCode TSTrajectorySetFromOptions_Memory(PetscOptions *PetscOptionsObject,TSTrajectory tj) 267 { 268 Stack *s = (Stack*)tj->data; 269 PetscErrorCode ierr; 270 271 PetscFunctionBegin; 272 ierr = PetscOptionsHead(PetscOptionsObject,"Memory based TS trajectory options");CHKERRQ(ierr); 273 { 274 ierr = PetscOptionsInt("-tstrajectory_max_cps","Maximum number of checkpoints","TSTrajectorySetMaxCheckpoints_Memory",s->max_cps,&s->max_cps,NULL);CHKERRQ(ierr); 275 ierr = PetscOptionsInt("-tstrajectory_stride","Stride to save checkpoints to file","TSTrajectorySetStride_Memory",s->stride,&s->stride,NULL);CHKERRQ(ierr); 276 } 277 ierr = PetscOptionsTail();CHKERRQ(ierr); 278 PetscFunctionReturn(0); 279 } 280 281 #undef __FUNCT__ 282 #define __FUNCT__ "TSTrajectorySetUp_Memory" 283 PetscErrorCode TSTrajectorySetUp_Memory(TSTrajectory tj,TS ts) 284 { 285 Stack *s = (Stack*)tj->data; 286 RevolveCTX *rctx; 287 PetscInt numY; 288 PetscErrorCode ierr; 289 290 PetscFunctionBegin; 291 292 ierr = TSGetStages(ts,&numY,NULL);CHKERRQ(ierr); 293 s->max_steps = PetscMin(ts->max_steps,(PetscInt)(ceil(ts->max_time/ts->time_step))); 294 295 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->max_steps)) { 296 s->userevolve = PETSC_TRUE; 297 ierr = PetscCalloc1(1,&rctx);CHKERRQ(ierr); 298 s->rctx = rctx; 299 rctx->snaps_in = s->max_cps; /* for theta methods snaps_in=2*max_cps */ 300 rctx->reverseonestep = PETSC_FALSE; 301 rctx->check = -1; 302 rctx->oldcapo = 0; 303 rctx->capo = 0; 304 rctx->info = 2; 305 if (s->stride>1) rctx->fine = s->stride; 306 else rctx->fine = s->max_steps; 307 ierr = StackCreate(PetscObjectComm((PetscObject)ts),s,s->max_cps,numY);CHKERRQ(ierr); 308 } else { 309 s->userevolve = PETSC_FALSE; 310 if (s->stride>1) { 311 ierr = StackCreate(PetscObjectComm((PetscObject)ts),s,s->stride+1,numY);CHKERRQ(ierr); 312 } else { 313 ierr = StackCreate(PetscObjectComm((PetscObject)ts),s,ts->max_steps+1,numY);CHKERRQ(ierr); 314 } 315 } 316 PetscFunctionReturn(0); 317 } 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "TSTrajectorySet_Memory" 321 PetscErrorCode TSTrajectorySet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X) 322 { 323 PetscInt i; 324 Vec *Y; 325 PetscReal timeprev; 326 StackElement e; 327 Stack *s = (Stack*)tj->data; 328 RevolveCTX *rctx; 329 PetscInt whattodo,rank,localstepnum,id; 330 PetscErrorCode ierr; 331 332 PetscFunctionBegin; 333 if (stepnum<s->top) SETERRQ(s->comm,PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 334 335 if (s->stride>1) { 336 localstepnum = stepnum%s->stride; 337 if (s->userevolve && stepnum!=0 && localstepnum==0 && stepnum!=s->max_steps) { /* first turn point */ 338 id = stepnum/s->stride; 339 ierr = StackDumpAll(s,id);CHKERRQ(ierr); 340 s->top = -1; /* reset top */ 341 rctx = s->rctx; 342 rctx->check = -1; 343 rctx->capo = 0; 344 rctx->fine = s->stride; 345 } 346 } else { 347 localstepnum = stepnum; 348 } 349 350 if (s->userevolve) { 351 rctx = s->rctx; 352 if (rctx->reverseonestep) { 353 rctx->reverseonestep = PETSC_FALSE; 354 PetscFunctionReturn(0); 355 } 356 if (rctx->stepsleft==0) { /* let the controller determine what to do next */ 357 rctx->capo = localstepnum; 358 rctx->oldcapo = rctx->capo; 359 whattodo = wrap_revolve(&rctx->check,&rctx->capo,&rctx->fine,&rctx->snaps_in,&rctx->info,&rank); 360 #ifdef PRINTWHATTODO 361 printwhattodo(whattodo,rctx); 362 #endif 363 if (whattodo==-1) SETERRQ(s->comm,PETSC_ERR_MEMC,"Error in the controller"); 364 if (whattodo==1) { 365 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 366 PetscFunctionReturn(0); /* do not need to checkpoint */ 367 } 368 if (whattodo==3 || whattodo==4) { 369 rctx->reverseonestep = PETSC_TRUE; 370 PetscFunctionReturn(0); 371 } 372 if (whattodo==5) { 373 rctx->oldcapo = rctx->capo; 374 whattodo = wrap_revolve(&rctx->check,&rctx->capo,&rctx->fine,&rctx->snaps_in,&rctx->info,&rank); /* must return 1*/ 375 #ifdef PRINTWHATTODO 376 printwhattodo(whattodo,rctx); 377 #endif 378 rctx->stepsleft = rctx->capo-rctx->oldcapo; 379 PetscFunctionReturn(0); 380 } 381 if (whattodo==2) { 382 rctx->oldcapo = rctx->capo; 383 whattodo = wrap_revolve(&rctx->check,&rctx->capo,&rctx->fine,&rctx->snaps_in,&rctx->info,&rank); /* must return 1*/ 384 #ifdef PRINTWHATTODO 385 printwhattodo(whattodo,rctx); 386 #endif 387 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 388 } 389 } else { /* advance s->stepsleft time steps without checkpointing */ 390 rctx->stepsleft--; 391 PetscFunctionReturn(0); 392 } 393 } 394 395 if (!s->userevolve && stepnum==0) PetscFunctionReturn(0); 396 397 /* checkpoint to memmory */ 398 if (localstepnum==s->top) { /* overwrite the top checkpoint */ 399 ierr = StackTop(s,&e); 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); 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 if (stepnum == 0) { 421 e->timeprev = e->time - ts->time_step; /* for consistency */ 422 } else { 423 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 424 e->timeprev = timeprev; 425 } 426 ierr = StackPush(s,e);CHKERRQ(ierr); 427 } 428 if (!s->userevolve && stepnum!=0 && localstepnum==0 && stepnum!=s->max_steps) { 429 id = stepnum/s->stride; 430 ierr = StackDumpAll(s,id);CHKERRQ(ierr); 431 s->top = -1; /* reset top */ 432 } 433 PetscFunctionReturn(0); 434 } 435 436 #undef __FUNCT__ 437 #define __FUNCT__ "TSTrajectoryGet_Memory" 438 PetscErrorCode TSTrajectoryGet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t) 439 { 440 Vec *Y; 441 PetscInt i; 442 StackElement e; 443 Stack *s = (Stack*)tj->data; 444 RevolveCTX *rctx; 445 PetscReal stepsize; 446 PetscInt whattodo,rank,localstepnum,id; 447 PetscErrorCode ierr; 448 449 PetscFunctionBegin; 450 if (s->stride>1) { 451 localstepnum = stepnum%s->stride; 452 if (localstepnum==0 && stepnum!=0 && stepnum!=s->max_steps) { 453 id = stepnum/s->stride; 454 ierr = StackLoadAll(ts,s,id);CHKERRQ(ierr); 455 ierr = StackPop(s,&e);CHKERRQ(ierr); /* pop out stepnum 0 */ 456 s->top = s->stride-1; 457 if (s->userevolve) { 458 rctx = s->rctx; 459 rctx->reverseonestep = PETSC_TRUE; 460 rctx->check = s->top; 461 } 462 } 463 } 464 if (s->userevolve) rctx = s->rctx; 465 if (s->userevolve && rctx->reverseonestep) { 466 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 467 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); /* go backward */ 468 rctx->reverseonestep = PETSC_FALSE; 469 PetscFunctionReturn(0); 470 } 471 472 /* restore a checkpoint */ 473 ierr = StackTop(s,&e);CHKERRQ(ierr); 474 ierr = VecCopy(e->X,ts->vec_sol);CHKERRQ(ierr); 475 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 476 for (i=0;i<s->numY;i++) { 477 ierr = VecCopy(e->Y[i],Y[i]);CHKERRQ(ierr); 478 } 479 *t = e->time; 480 481 if (e->stepnum < stepnum) { /* need recomputation */ 482 rctx->capo = stepnum; 483 whattodo = wrap_revolve(&rctx->check,&rctx->capo,&rctx->fine,&rctx->snaps_in,&rctx->info,&rank); 484 #ifdef PRINTWHATTODO 485 printwhattodo(whattodo,rctx); 486 #endif 487 ierr = TSSetTimeStep(ts,(*t)-e->timeprev);CHKERRQ(ierr); 488 /* reset ts context */ 489 PetscInt steps = ts->steps; 490 ts->steps = e->stepnum; 491 ts->ptime = e->time; 492 ts->ptime_prev = e->timeprev; 493 for (i=e->stepnum;i<stepnum;i++) { /* assume fixed step size */ 494 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 495 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 496 ierr = TSStep(ts);CHKERRQ(ierr); 497 if (ts->event) { 498 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 499 } 500 if (!ts->steprollback) { 501 ierr = TSPostStep(ts);CHKERRQ(ierr); 502 } 503 } 504 /* reverseonestep must be true after the for loop */ 505 ts->steps = steps; 506 ts->total_steps = stepnum; 507 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 508 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); /* go backward */ 509 if (stepnum-e->stepnum==1) { 510 ierr = StackPop(s,&e);CHKERRQ(ierr); 511 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 512 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 513 ierr = PetscFree(e);CHKERRQ(ierr); 514 } 515 rctx->reverseonestep = PETSC_FALSE; 516 } else if (e->stepnum == stepnum) { 517 ierr = TSSetTimeStep(ts,-(*t)+e->timeprev);CHKERRQ(ierr); /* go backward */ 518 ierr = StackPop(s,&e);CHKERRQ(ierr); 519 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 520 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 521 ierr = PetscFree(e);CHKERRQ(ierr); 522 } else { 523 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); 524 } 525 526 PetscFunctionReturn(0); 527 } 528 529 #undef __FUNCT__ 530 #define __FUNCT__ "TSTrajectoryDestroy_Memory" 531 PETSC_EXTERN PetscErrorCode TSTrajectoryDestroy_Memory(TSTrajectory tj) 532 { 533 Stack *s = (Stack*)tj->data; 534 PetscErrorCode ierr; 535 536 PetscFunctionBegin; 537 ierr = StackDestroy(s);CHKERRQ(ierr); 538 PetscFunctionReturn(0); 539 } 540 541 /*MC 542 TSTRAJECTORYMEMORY - Stores each solution of the ODE/ADE in memory 543 544 Level: intermediate 545 546 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType() 547 548 M*/ 549 #undef __FUNCT__ 550 #define __FUNCT__ "TSTrajectoryCreate_Memory" 551 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory tj,TS ts) 552 { 553 Stack *s; 554 PetscErrorCode ierr; 555 556 PetscFunctionBegin; 557 tj->ops->set = TSTrajectorySet_Memory; 558 tj->ops->get = TSTrajectoryGet_Memory; 559 tj->ops->setup = TSTrajectorySetUp_Memory; 560 tj->ops->destroy = TSTrajectoryDestroy_Memory; 561 tj->ops->setfromoptions = TSTrajectorySetFromOptions_Memory; 562 563 ierr = PetscCalloc1(1,&s);CHKERRQ(ierr); 564 s->max_cps = -1; /* -1 indicates that it is not set */ 565 s->stride = 0; /* if not zero, two-level checkpointing will be used */ 566 tj->data = s; 567 PetscFunctionReturn(0); 568 } 569