1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscsys.h> 3 #if defined(PETSC_HAVE_REVOLVE) 4 #include <revolve_c.h> 5 #endif 6 7 PetscLogEvent TSTrajectory_DiskWrite, TSTrajectory_DiskRead; 8 static PetscErrorCode TSTrajectorySet_Memory(TSTrajectory,TS,PetscInt,PetscReal,Vec); 9 10 typedef enum {NONE,TWO_LEVEL_NOREVOLVE,TWO_LEVEL_REVOLVE,TWO_LEVEL_TWO_REVOLVE,REVOLVE_OFFLINE,REVOLVE_ONLINE,REVOLVE_MULTISTAGE} SchedulerType; 11 12 typedef struct _StackElement { 13 PetscInt stepnum; 14 Vec X; 15 Vec *Y; 16 PetscReal time; 17 PetscReal timeprev; /* for no solution_only mode */ 18 PetscReal timenext; /* for solution_only mode */ 19 } *StackElement; 20 21 #if defined(PETSC_HAVE_REVOLVE) 22 typedef struct _RevolveCTX { 23 PetscBool reverseonestep; 24 PetscInt where; 25 PetscInt snaps_in; 26 PetscInt stepsleft; 27 PetscInt check; 28 PetscInt oldcapo; 29 PetscInt capo; 30 PetscInt fine; 31 PetscInt info; 32 } RevolveCTX; 33 #endif 34 35 typedef struct _Stack { 36 PetscInt stacksize; 37 PetscInt top; 38 StackElement *container; 39 PetscInt numY; 40 PetscBool solution_only; 41 PetscBool use_dram; 42 } Stack; 43 44 typedef struct _DiskStack { 45 PetscInt stacksize; 46 PetscInt top; 47 PetscInt *container; 48 } DiskStack; 49 50 typedef struct _TJScheduler { 51 SchedulerType stype; 52 #if defined(PETSC_HAVE_REVOLVE) 53 RevolveCTX *rctx,*rctx2; 54 PetscBool use_online; 55 PetscInt store_stride; 56 #endif 57 PetscBool recompute; 58 PetscBool skip_trajectory; 59 PetscBool save_stack; 60 PetscInt max_cps_ram; /* maximum checkpoints in RAM */ 61 PetscInt max_cps_disk; /* maximum checkpoints on disk */ 62 PetscInt stride; 63 PetscInt total_steps; /* total number of steps */ 64 Stack stack; 65 DiskStack diskstack; 66 } TJScheduler; 67 68 static PetscErrorCode TurnForwardWithStepsize(TS ts,PetscReal nextstepsize) 69 { 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 /* reverse the direction */ 74 ierr = TSSetTimeStep(ts,nextstepsize);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 static PetscErrorCode TurnForward(TS ts) 79 { 80 PetscReal stepsize; 81 PetscErrorCode ierr; 82 83 PetscFunctionBegin; 84 /* reverse the direction */ 85 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 86 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 static PetscErrorCode TurnBackward(TS ts) 91 { 92 PetscReal stepsize; 93 PetscErrorCode ierr; 94 95 PetscFunctionBegin; 96 if (!ts->trajectory->adjoint_solve_mode) PetscFunctionReturn(0); 97 /* reverse the direction */ 98 stepsize = ts->ptime_prev-ts->ptime; 99 ierr = TSSetTimeStep(ts,stepsize);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 103 static PetscErrorCode ElementCreate(TS ts,Stack *stack,StackElement *e) 104 { 105 Vec X; 106 Vec *Y; 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 if (stack->use_dram) { 111 ierr = PetscMallocSetDRAM();CHKERRQ(ierr); 112 } 113 ierr = PetscCalloc1(1,e);CHKERRQ(ierr); 114 ierr = TSGetSolution(ts,&X);CHKERRQ(ierr); 115 ierr = VecDuplicate(X,&(*e)->X);CHKERRQ(ierr); 116 if (stack->numY > 0 && !stack->solution_only) { 117 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 118 ierr = VecDuplicateVecs(Y[0],stack->numY,&(*e)->Y);CHKERRQ(ierr); 119 } 120 if (stack->use_dram) { 121 ierr = PetscMallocResetDRAM();CHKERRQ(ierr); 122 } 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode ElementSet(TS ts,Stack *stack,StackElement *e,PetscInt stepnum,PetscReal time,Vec X) 127 { 128 Vec *Y; 129 PetscInt i; 130 PetscReal timeprev; 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 ierr = VecCopy(X,(*e)->X);CHKERRQ(ierr); 135 if (stack->numY > 0 && !stack->solution_only) { 136 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 137 for (i=0;i<stack->numY;i++) { 138 ierr = VecCopy(Y[i],(*e)->Y[i]);CHKERRQ(ierr); 139 } 140 } 141 (*e)->stepnum = stepnum; 142 (*e)->time = time; 143 /* for consistency */ 144 if (stepnum == 0) { 145 (*e)->timeprev = (*e)->time - ts->time_step; 146 } else { 147 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 148 (*e)->timeprev = timeprev; 149 } 150 PetscFunctionReturn(0); 151 } 152 153 static PetscErrorCode ElementDestroy(Stack *stack,StackElement e) 154 { 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 if (stack->use_dram) { 159 ierr = PetscMallocSetDRAM();CHKERRQ(ierr); 160 } 161 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 162 if (stack->numY > 0 && !stack->solution_only) { 163 ierr = VecDestroyVecs(stack->numY,&e->Y);CHKERRQ(ierr); 164 } 165 ierr = PetscFree(e);CHKERRQ(ierr); 166 if (stack->use_dram) { 167 ierr = PetscMallocResetDRAM();CHKERRQ(ierr); 168 } 169 PetscFunctionReturn(0); 170 } 171 172 static PetscErrorCode StackResize(Stack *stack,PetscInt newsize) 173 { 174 StackElement *newcontainer; 175 PetscInt i; 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 ierr = PetscMalloc1(newsize*sizeof(StackElement),&newcontainer);CHKERRQ(ierr); 180 for (i=0;i<stack->stacksize;i++) { 181 newcontainer[i] = stack->container[i]; 182 } 183 ierr = PetscFree(stack->container);CHKERRQ(ierr); 184 stack->container = newcontainer; 185 stack->stacksize = newsize; 186 PetscFunctionReturn(0); 187 } 188 189 static PetscErrorCode StackPush(Stack *stack,StackElement e) 190 { 191 PetscFunctionBegin; 192 if (stack->top+1 >= stack->stacksize) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MEMC,"Maximum stack size (%D) exceeded",stack->stacksize); 193 stack->container[++stack->top] = e; 194 PetscFunctionReturn(0); 195 } 196 197 static PetscErrorCode StackPop(Stack *stack,StackElement *e) 198 { 199 PetscFunctionBegin; 200 *e = NULL; 201 if (stack->top == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEMC,"Empty stack"); 202 *e = stack->container[stack->top--]; 203 PetscFunctionReturn(0); 204 } 205 206 static PetscErrorCode StackTop(Stack *stack,StackElement *e) 207 { 208 PetscFunctionBegin; 209 *e = stack->container[stack->top]; 210 PetscFunctionReturn(0); 211 } 212 213 static PetscErrorCode StackCreate(Stack *stack,PetscInt size,PetscInt ny) 214 { 215 PetscErrorCode ierr; 216 217 PetscFunctionBegin; 218 stack->top = -1; 219 stack->numY = ny; 220 221 ierr = PetscMalloc1(size*sizeof(StackElement),&stack->container);CHKERRQ(ierr); 222 PetscFunctionReturn(0); 223 } 224 225 static PetscErrorCode StackDestroy(Stack *stack) 226 { 227 PetscInt i,n; 228 StackElement e = NULL; 229 PetscErrorCode ierr; 230 231 PetscFunctionBegin; 232 if (!stack->container) PetscFunctionReturn(0); 233 if (stack->top > -1) { 234 n = stack->top+1; /* number of elements in the stack */ 235 for (i=0;i<n;i++) { 236 ierr = StackPop(stack,&e);CHKERRQ(ierr); 237 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 238 } 239 } 240 ierr = PetscFree(stack->container);CHKERRQ(ierr); 241 PetscFunctionReturn(0); 242 } 243 244 static PetscErrorCode StackFind(Stack *stack,StackElement *e,PetscInt index) 245 { 246 PetscFunctionBegin; 247 *e = NULL; 248 if (index < 0 || index > stack->top) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid index %D",index); 249 *e = stack->container[index]; 250 PetscFunctionReturn(0); 251 } 252 253 static PetscErrorCode OutputBIN(MPI_Comm comm,const char *filename,PetscViewer *viewer) 254 { 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 ierr = PetscViewerCreate(comm,viewer);CHKERRQ(ierr); 259 ierr = PetscViewerSetType(*viewer,PETSCVIEWERBINARY);CHKERRQ(ierr); 260 ierr = PetscViewerPushFormat(*viewer,PETSC_VIEWER_NATIVE);CHKERRQ(ierr); 261 ierr = PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);CHKERRQ(ierr); 262 ierr = PetscViewerFileSetName(*viewer,filename);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 static PetscErrorCode WriteToDisk(PetscInt stepnum,PetscReal time,PetscReal timeprev,Vec X,Vec *Y,PetscInt numY,PetscBool solution_only,PetscViewer viewer) 267 { 268 PetscInt i; 269 PetscErrorCode ierr; 270 271 PetscFunctionBegin; 272 ierr = PetscViewerBinaryWrite(viewer,&stepnum,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 273 ierr = VecView(X,viewer);CHKERRQ(ierr); 274 for (i=0;!solution_only && i<numY;i++) { 275 ierr = VecView(Y[i],viewer);CHKERRQ(ierr); 276 } 277 ierr = PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 278 ierr = PetscViewerBinaryWrite(viewer,&timeprev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 279 PetscFunctionReturn(0); 280 } 281 282 static PetscErrorCode ReadFromDisk(PetscInt *stepnum,PetscReal *time,PetscReal *timeprev,Vec X,Vec *Y,PetscInt numY,PetscBool solution_only,PetscViewer viewer) 283 { 284 PetscInt i; 285 PetscErrorCode ierr; 286 287 PetscFunctionBegin; 288 ierr = PetscViewerBinaryRead(viewer,stepnum,1,NULL,PETSC_INT);CHKERRQ(ierr); 289 ierr = VecLoad(X,viewer);CHKERRQ(ierr); 290 for (i=0;!solution_only && i<numY;i++) { 291 ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr); 292 } 293 ierr = PetscViewerBinaryRead(viewer,time,1,NULL,PETSC_REAL);CHKERRQ(ierr); 294 ierr = PetscViewerBinaryRead(viewer,timeprev,1,NULL,PETSC_REAL);CHKERRQ(ierr); 295 PetscFunctionReturn(0); 296 } 297 298 static PetscErrorCode StackDumpAll(TSTrajectory tj,TS ts,Stack *stack,PetscInt id) 299 { 300 Vec *Y; 301 PetscInt i; 302 StackElement e = NULL; 303 PetscViewer viewer; 304 char filename[PETSC_MAX_PATH_LEN]; 305 PetscErrorCode ierr; 306 MPI_Comm comm; 307 308 PetscFunctionBegin; 309 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 310 if (tj->monitor) { 311 ierr = PetscViewerASCIIPushTab(tj->monitor);CHKERRQ(ierr); 312 ierr = PetscViewerASCIIPrintf(tj->monitor,"Dump stack id %D to file\n",id);CHKERRQ(ierr); 313 ierr = PetscViewerASCIIPopTab(tj->monitor);CHKERRQ(ierr); 314 } 315 if (id == 1) { 316 PetscMPIInt rank; 317 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 318 if (!rank) { 319 ierr = PetscRMTree("SA-data");CHKERRQ(ierr); 320 ierr = PetscMkdir("SA-data");CHKERRQ(ierr); 321 } 322 } 323 ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-STACK%06d.bin",id);CHKERRQ(ierr); 324 ierr = OutputBIN(comm,filename,&viewer);CHKERRQ(ierr); 325 for (i=0;i<stack->stacksize;i++) { 326 e = stack->container[i]; 327 ierr = PetscLogEventBegin(TSTrajectory_DiskWrite,tj,ts,0,0);CHKERRQ(ierr); 328 ierr = WriteToDisk(e->stepnum,e->time,e->timeprev,e->X,e->Y,stack->numY,stack->solution_only,viewer);CHKERRQ(ierr); 329 ierr = PetscLogEventEnd(TSTrajectory_DiskWrite,tj,ts,0,0);CHKERRQ(ierr); 330 ts->trajectory->diskwrites++; 331 } 332 /* save the last step for restart, the last step is in memory when using single level schemes, but not necessarily the case for multi level schemes */ 333 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 334 ierr = PetscLogEventBegin(TSTrajectory_DiskWrite,tj,ts,0,0);CHKERRQ(ierr); 335 ierr = WriteToDisk(ts->steps,ts->ptime,ts->ptime_prev,ts->vec_sol,Y,stack->numY,stack->solution_only,viewer);CHKERRQ(ierr); 336 ierr = PetscLogEventEnd(TSTrajectory_DiskWrite,tj,ts,0,0);CHKERRQ(ierr); 337 ts->trajectory->diskwrites++; 338 for (i=0;i<stack->stacksize;i++) { 339 ierr = StackPop(stack,&e);CHKERRQ(ierr); 340 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 341 } 342 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 343 PetscFunctionReturn(0); 344 } 345 346 static PetscErrorCode StackLoadAll(TSTrajectory tj,TS ts,Stack *stack,PetscInt id) 347 { 348 Vec *Y; 349 PetscInt i; 350 StackElement e; 351 PetscViewer viewer; 352 char filename[PETSC_MAX_PATH_LEN]; 353 PetscErrorCode ierr; 354 355 PetscFunctionBegin; 356 if (tj->monitor) { 357 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 358 ierr = PetscViewerASCIIPrintf(tj->monitor,"Load stack from file\n");CHKERRQ(ierr); 359 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 360 } 361 ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-STACK%06d.bin",id);CHKERRQ(ierr); 362 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 363 for (i=0;i<stack->stacksize;i++) { 364 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 365 ierr = StackPush(stack,e);CHKERRQ(ierr); 366 ierr = PetscLogEventBegin(TSTrajectory_DiskRead,tj,ts,0,0);CHKERRQ(ierr); 367 ierr = ReadFromDisk(&e->stepnum,&e->time,&e->timeprev,e->X,e->Y,stack->numY,stack->solution_only,viewer);CHKERRQ(ierr); 368 ierr = PetscLogEventEnd(TSTrajectory_DiskRead,tj,ts,0,0);CHKERRQ(ierr); 369 ts->trajectory->diskreads++; 370 } 371 /* load the last step into TS */ 372 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 373 ierr = PetscLogEventBegin(TSTrajectory_DiskRead,tj,ts,0,0);CHKERRQ(ierr); 374 ierr = ReadFromDisk(&ts->steps,&ts->ptime,&ts->ptime_prev,ts->vec_sol,Y,stack->numY,stack->solution_only,viewer);CHKERRQ(ierr); 375 ierr = PetscLogEventEnd(TSTrajectory_DiskRead,tj,ts,0,0);CHKERRQ(ierr); 376 ts->trajectory->diskreads++; 377 ierr = TurnBackward(ts);CHKERRQ(ierr); 378 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 379 PetscFunctionReturn(0); 380 } 381 382 #if defined(PETSC_HAVE_REVOLVE) 383 static PetscErrorCode StackLoadLast(TSTrajectory tj,TS ts,Stack *stack,PetscInt id) 384 { 385 Vec *Y; 386 PetscInt size; 387 PetscViewer viewer; 388 char filename[PETSC_MAX_PATH_LEN]; 389 #if defined(PETSC_HAVE_MPIIO) 390 PetscBool usempiio; 391 #endif 392 int fd; 393 off_t off,offset; 394 PetscErrorCode ierr; 395 396 PetscFunctionBegin; 397 if (tj->monitor) { 398 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 399 ierr = PetscViewerASCIIPrintf(tj->monitor,"Load last stack element from file\n");CHKERRQ(ierr); 400 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 401 } 402 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 403 ierr = VecGetSize(Y[0],&size);CHKERRQ(ierr); 404 /* VecView writes to file two extra int's for class id and number of rows */ 405 off = -((stack->solution_only?0:stack->numY)+1)*(size*PETSC_BINARY_SCALAR_SIZE+2*PETSC_BINARY_INT_SIZE)-PETSC_BINARY_INT_SIZE-2*PETSC_BINARY_SCALAR_SIZE; 406 407 ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-STACK%06d.bin",id);CHKERRQ(ierr); 408 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 409 #if defined(PETSC_HAVE_MPIIO) 410 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&usempiio);CHKERRQ(ierr); 411 if (usempiio) { 412 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,(MPI_File*)&fd);CHKERRQ(ierr); 413 ierr = PetscBinarySynchronizedSeek(PETSC_COMM_WORLD,fd,off,PETSC_BINARY_SEEK_END,&offset);CHKERRQ(ierr); 414 } else { 415 #endif 416 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 417 ierr = PetscBinarySeek(fd,off,PETSC_BINARY_SEEK_END,&offset);CHKERRQ(ierr); 418 #if defined(PETSC_HAVE_MPIIO) 419 } 420 #endif 421 /* load the last step into TS */ 422 ierr = PetscLogEventBegin(TSTrajectory_DiskRead,tj,ts,0,0);CHKERRQ(ierr); 423 ierr = ReadFromDisk(&ts->steps,&ts->ptime,&ts->ptime_prev,ts->vec_sol,Y,stack->numY,stack->solution_only,viewer);CHKERRQ(ierr); 424 ierr = PetscLogEventEnd(TSTrajectory_DiskRead,tj,ts,0,0);CHKERRQ(ierr); 425 ts->trajectory->diskreads++; 426 ierr = TurnBackward(ts);CHKERRQ(ierr); 427 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 428 PetscFunctionReturn(0); 429 } 430 #endif 431 432 static PetscErrorCode DumpSingle(TSTrajectory tj,TS ts,Stack *stack,PetscInt id) 433 { 434 Vec *Y; 435 PetscInt stepnum; 436 PetscViewer viewer; 437 char filename[PETSC_MAX_PATH_LEN]; 438 PetscErrorCode ierr; 439 MPI_Comm comm; 440 441 PetscFunctionBegin; 442 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 443 if (tj->monitor) { 444 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 445 ierr = PetscViewerASCIIPrintf(tj->monitor,"Dump a single point from file\n");CHKERRQ(ierr); 446 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 447 } 448 ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr); 449 if (id == 1) { 450 PetscMPIInt rank; 451 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 452 if (!rank) { 453 ierr = PetscRMTree("SA-data");CHKERRQ(ierr); 454 ierr = PetscMkdir("SA-data");CHKERRQ(ierr); 455 } 456 } 457 ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-CPS%06d.bin",id);CHKERRQ(ierr); 458 ierr = OutputBIN(comm,filename,&viewer);CHKERRQ(ierr); 459 460 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 461 ierr = PetscLogEventBegin(TSTrajectory_DiskWrite,tj,ts,0,0);CHKERRQ(ierr); 462 ierr = WriteToDisk(stepnum,ts->ptime,ts->ptime_prev,ts->vec_sol,Y,stack->numY,stack->solution_only,viewer);CHKERRQ(ierr); 463 ierr = PetscLogEventEnd(TSTrajectory_DiskWrite,tj,ts,0,0);CHKERRQ(ierr); 464 ts->trajectory->diskwrites++; 465 466 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 467 PetscFunctionReturn(0); 468 } 469 470 static PetscErrorCode LoadSingle(TSTrajectory tj,TS ts,Stack *stack,PetscInt id) 471 { 472 Vec *Y; 473 PetscViewer viewer; 474 char filename[PETSC_MAX_PATH_LEN]; 475 PetscErrorCode ierr; 476 477 PetscFunctionBegin; 478 if (tj->monitor) { 479 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 480 ierr = PetscViewerASCIIPrintf(tj->monitor,"Load a single point from file\n");CHKERRQ(ierr); 481 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 482 } 483 ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-CPS%06d.bin",id);CHKERRQ(ierr); 484 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 485 486 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 487 ierr = PetscLogEventBegin(TSTrajectory_DiskRead,tj,ts,0,0);CHKERRQ(ierr); 488 ierr = ReadFromDisk(&ts->steps,&ts->ptime,&ts->ptime_prev,ts->vec_sol,Y,stack->numY,stack->solution_only,viewer);CHKERRQ(ierr); 489 ierr = PetscLogEventEnd(TSTrajectory_DiskRead,tj,ts,0,0);CHKERRQ(ierr); 490 ts->trajectory->diskreads++; 491 492 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 493 PetscFunctionReturn(0); 494 } 495 496 static PetscErrorCode UpdateTS(TS ts,Stack *stack,StackElement e, PetscBool adjoint_mode) 497 { 498 Vec *Y; 499 PetscInt i; 500 PetscErrorCode ierr; 501 502 PetscFunctionBegin; 503 ierr = VecCopy(e->X,ts->vec_sol);CHKERRQ(ierr); 504 if (!stack->solution_only) { 505 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 506 for (i=0;i<stack->numY;i++) { 507 ierr = VecCopy(e->Y[i],Y[i]);CHKERRQ(ierr); 508 } 509 } 510 if (adjoint_mode) { 511 ierr = TSSetTimeStep(ts,e->timeprev-e->time);CHKERRQ(ierr); /* stepsize will be negative */ 512 } else { 513 ierr = TSSetTimeStep(ts,e->time-e->timeprev);CHKERRQ(ierr); /* stepsize will be positive */ 514 } 515 ts->ptime = e->time; 516 ts->ptime_prev = e->timeprev; 517 PetscFunctionReturn(0); 518 } 519 520 static PetscErrorCode ReCompute(TS ts,TJScheduler *tjsch,PetscInt stepnumbegin,PetscInt stepnumend) 521 { 522 Stack *stack = &tjsch->stack; 523 PetscInt i; 524 PetscErrorCode ierr; 525 526 PetscFunctionBegin; 527 ierr = TSSetStepNumber(ts,stepnumbegin);CHKERRQ(ierr);/* global step number */ 528 for (i=stepnumbegin;i<stepnumend;i++) { /* assume fixed step size */ 529 if (stack->solution_only && !tjsch->skip_trajectory) { /* revolve online need this */ 530 /* don't use the public interface as it will update the TSHistory: this need a better fix */ 531 ierr = TSTrajectorySet_Memory(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 532 } 533 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 534 ierr = TSStep(ts);CHKERRQ(ierr); 535 if (!stack->solution_only && !tjsch->skip_trajectory) { 536 /* don't use the public interface as it will update the TSHistory: this need a better fix */ 537 ierr = TSTrajectorySet_Memory(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 538 } 539 ierr = TSEventHandler(ts);CHKERRQ(ierr); 540 if (!ts->steprollback) { 541 ierr = TSPostStep(ts);CHKERRQ(ierr); 542 } 543 } 544 ierr = TurnBackward(ts);CHKERRQ(ierr); 545 ts->trajectory->recomps += stepnumend-stepnumbegin; /* recomputation counter */ 546 ierr = TSSetStepNumber(ts,stepnumend);CHKERRQ(ierr); 547 PetscFunctionReturn(0); 548 } 549 550 static PetscErrorCode TopLevelStore(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscInt localstepnum,PetscInt laststridesize,PetscBool *done) 551 { 552 Stack *stack = &tjsch->stack; 553 DiskStack *diskstack = &tjsch->diskstack; 554 PetscInt stridenum; 555 PetscErrorCode ierr; 556 557 PetscFunctionBegin; 558 *done = PETSC_FALSE; 559 stridenum = stepnum/tjsch->stride; 560 /* make sure saved checkpoint id starts from 1 561 skip last stride when using stridenum+1 562 skip first stride when using stridenum */ 563 if (stack->solution_only) { 564 if (tjsch->save_stack) { 565 if (localstepnum == tjsch->stride-1 && stepnum < tjsch->total_steps-laststridesize) { /* current step will be saved without going through stack */ 566 ierr = StackDumpAll(tj,ts,stack,stridenum+1);CHKERRQ(ierr); 567 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) diskstack->container[++diskstack->top] = stridenum+1; 568 *done = PETSC_TRUE; 569 } 570 } else { 571 if (localstepnum == 0 && stepnum < tjsch->total_steps-laststridesize) { 572 ierr = DumpSingle(tj,ts,stack,stridenum+1);CHKERRQ(ierr); 573 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) diskstack->container[++diskstack->top] = stridenum+1; 574 *done = PETSC_TRUE; 575 } 576 } 577 } else { 578 if (tjsch->save_stack) { 579 if (localstepnum == 0 && stepnum < tjsch->total_steps && stepnum != 0) { /* skip the first stride */ 580 ierr = StackDumpAll(tj,ts,stack,stridenum);CHKERRQ(ierr); 581 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) diskstack->container[++diskstack->top] = stridenum; 582 *done = PETSC_TRUE; 583 } 584 } else { 585 if (localstepnum == 1 && stepnum < tjsch->total_steps-laststridesize) { 586 ierr = DumpSingle(tj,ts,stack,stridenum+1);CHKERRQ(ierr); 587 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) diskstack->container[++diskstack->top] = stridenum+1; 588 *done = PETSC_TRUE; 589 } 590 } 591 } 592 PetscFunctionReturn(0); 593 } 594 595 static PetscErrorCode SetTrajN(TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 596 { 597 Stack *stack = &tjsch->stack; 598 StackElement e; 599 PetscErrorCode ierr; 600 601 PetscFunctionBegin; 602 /* skip the last step */ 603 if (ts->reason) { /* only affect the forward run */ 604 /* update total_steps in the end of forward run */ 605 if (stepnum != tjsch->total_steps) tjsch->total_steps = stepnum; 606 if (stack->solution_only) { 607 /* get rid of the solution at second last step */ 608 ierr = StackPop(stack,&e);CHKERRQ(ierr); 609 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 610 } 611 PetscFunctionReturn(0); 612 } 613 /* do not save trajectory at the recompute stage for solution_only mode */ 614 if (stack->solution_only && tjsch->recompute) PetscFunctionReturn(0); 615 /* skip the first step for no_solution_only mode */ 616 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 617 618 /* resize the stack */ 619 if (stack->top+1 == stack->stacksize) { 620 ierr = StackResize(stack,2*stack->stacksize);CHKERRQ(ierr); 621 } 622 /* update timenext for the previous step; necessary for step adaptivity */ 623 if (stack->top > -1) { 624 ierr = StackTop(stack,&e);CHKERRQ(ierr); 625 e->timenext = ts->ptime; 626 } 627 if (stepnum < stack->top) { 628 SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 629 } 630 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 631 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 632 ierr = StackPush(stack,e);CHKERRQ(ierr); 633 PetscFunctionReturn(0); 634 } 635 636 static PetscErrorCode SetTrajN_2(TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 637 { 638 Stack *stack = &tjsch->stack; 639 StackElement e; 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 if (stack->top+1 == stack->stacksize) { 644 ierr = StackResize(stack,2*stack->stacksize);CHKERRQ(ierr); 645 } 646 /* update timenext for the previous step; necessary for step adaptivity */ 647 if (stack->top > -1) { 648 ierr = StackTop(stack,&e);CHKERRQ(ierr); 649 e->timenext = ts->ptime; 650 } 651 if (stepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 652 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 653 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 654 ierr = StackPush(stack,e);CHKERRQ(ierr); 655 PetscFunctionReturn(0); 656 } 657 658 static PetscErrorCode GetTrajN(TS ts,TJScheduler *tjsch,PetscInt stepnum) 659 { 660 Stack *stack = &tjsch->stack; 661 StackElement e; 662 PetscInt ns; 663 PetscErrorCode ierr; 664 665 PetscFunctionBegin; 666 if (stepnum == tjsch->total_steps) { 667 ierr = TurnBackward(ts);CHKERRQ(ierr); 668 PetscFunctionReturn(0); 669 } 670 /* restore a checkpoint */ 671 ierr = StackTop(stack,&e);CHKERRQ(ierr); 672 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 673 ierr = TSGetStages(ts,&ns,NULL);CHKERRQ(ierr); 674 if (stack->solution_only && ns) { /* recompute one step */ 675 tjsch->recompute = PETSC_TRUE; 676 ierr = TurnForwardWithStepsize(ts,e->timenext-e->time);CHKERRQ(ierr); 677 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 678 } 679 ierr = StackPop(stack,&e);CHKERRQ(ierr); 680 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 681 PetscFunctionReturn(0); 682 } 683 684 static PetscErrorCode GetTrajN_2(TS ts,TJScheduler *tjsch,PetscInt stepnum) 685 { 686 Stack *stack = &tjsch->stack; 687 StackElement e = NULL; 688 PetscErrorCode ierr; 689 690 PetscFunctionBegin; 691 ierr = StackFind(stack,&e,stepnum);CHKERRQ(ierr); 692 if (stepnum != e->stepnum) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Inconsistent steps! %D != %D",stepnum,e->stepnum); 693 ierr = UpdateTS(ts,stack,e,PETSC_FALSE);CHKERRQ(ierr); 694 PetscFunctionReturn(0); 695 } 696 697 static PetscErrorCode SetTrajTLNR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 698 { 699 Stack *stack = &tjsch->stack; 700 PetscInt localstepnum,laststridesize; 701 StackElement e; 702 PetscBool done; 703 PetscErrorCode ierr; 704 705 PetscFunctionBegin; 706 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 707 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 708 if (tjsch->save_stack && tjsch->recompute) PetscFunctionReturn(0); 709 710 localstepnum = stepnum%tjsch->stride; 711 /* (stridesize-1) checkpoints are saved in each stride; an extra point is added by StackDumpAll() */ 712 laststridesize = tjsch->total_steps%tjsch->stride; 713 if (!laststridesize) laststridesize = tjsch->stride; 714 715 if (!tjsch->recompute) { 716 ierr = TopLevelStore(tj,ts,tjsch,stepnum,localstepnum,laststridesize,&done);CHKERRQ(ierr); 717 if (!tjsch->save_stack && stepnum < tjsch->total_steps-laststridesize) PetscFunctionReturn(0); 718 } 719 if (!stack->solution_only && localstepnum == 0) PetscFunctionReturn(0); /* skip last point in each stride at recompute stage or last stride */ 720 if (stack->solution_only && localstepnum == tjsch->stride-1) PetscFunctionReturn(0); /* skip last step in each stride at recompute stage or last stride */ 721 722 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 723 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 724 ierr = StackPush(stack,e);CHKERRQ(ierr); 725 PetscFunctionReturn(0); 726 } 727 728 static PetscErrorCode GetTrajTLNR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 729 { 730 Stack *stack = &tjsch->stack; 731 PetscInt id,localstepnum,laststridesize; 732 StackElement e; 733 PetscErrorCode ierr; 734 735 PetscFunctionBegin; 736 if (stepnum == tjsch->total_steps) { 737 ierr = TurnBackward(ts);CHKERRQ(ierr); 738 PetscFunctionReturn(0); 739 } 740 741 localstepnum = stepnum%tjsch->stride; 742 laststridesize = tjsch->total_steps%tjsch->stride; 743 if (!laststridesize) laststridesize = tjsch->stride; 744 if (stack->solution_only) { 745 /* fill stack with info */ 746 if (localstepnum == 0 && tjsch->total_steps-stepnum >= laststridesize) { 747 id = stepnum/tjsch->stride; 748 if (tjsch->save_stack) { 749 ierr = StackLoadAll(tj,ts,stack,id);CHKERRQ(ierr); 750 tjsch->recompute = PETSC_TRUE; 751 tjsch->skip_trajectory = PETSC_TRUE; 752 ierr = TurnForward(ts);CHKERRQ(ierr); 753 ierr = ReCompute(ts,tjsch,id*tjsch->stride-1,id*tjsch->stride);CHKERRQ(ierr); 754 tjsch->skip_trajectory = PETSC_FALSE; 755 } else { 756 ierr = LoadSingle(tj,ts,stack,id);CHKERRQ(ierr); 757 tjsch->recompute = PETSC_TRUE; 758 ierr = TurnForward(ts);CHKERRQ(ierr); 759 ierr = ReCompute(ts,tjsch,(id-1)*tjsch->stride,id*tjsch->stride);CHKERRQ(ierr); 760 } 761 PetscFunctionReturn(0); 762 } 763 /* restore a checkpoint */ 764 ierr = StackPop(stack,&e);CHKERRQ(ierr); 765 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 766 tjsch->recompute = PETSC_TRUE; 767 tjsch->skip_trajectory = PETSC_TRUE; 768 ierr = TurnForward(ts);CHKERRQ(ierr); 769 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 770 tjsch->skip_trajectory = PETSC_FALSE; 771 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 772 } else { 773 /* fill stack with info */ 774 if (localstepnum == 0 && tjsch->total_steps-stepnum >= laststridesize) { 775 id = stepnum/tjsch->stride; 776 if (tjsch->save_stack) { 777 ierr = StackLoadAll(tj,ts,stack,id);CHKERRQ(ierr); 778 } else { 779 ierr = LoadSingle(tj,ts,stack,id);CHKERRQ(ierr); 780 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 781 ierr = ElementSet(ts,stack,&e,(id-1)*tjsch->stride+1,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 782 ierr = StackPush(stack,e);CHKERRQ(ierr); 783 tjsch->recompute = PETSC_TRUE; 784 ierr = TurnForward(ts);CHKERRQ(ierr); 785 ierr = ReCompute(ts,tjsch,e->stepnum,id*tjsch->stride);CHKERRQ(ierr); 786 } 787 PetscFunctionReturn(0); 788 } 789 /* restore a checkpoint */ 790 ierr = StackPop(stack,&e);CHKERRQ(ierr); 791 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 792 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 793 } 794 PetscFunctionReturn(0); 795 } 796 797 #if defined(PETSC_HAVE_REVOLVE) 798 static PetscErrorCode printwhattodo(PetscViewer viewer,PetscInt whattodo,RevolveCTX *rctx,PetscInt shift) 799 { 800 PetscErrorCode ierr; 801 802 PetscFunctionBegin; 803 if (!viewer) PetscFunctionReturn(0); 804 805 switch(whattodo) { 806 case 1: 807 ierr = PetscViewerASCIIPrintf(viewer,"Advance from %D to %D\n",rctx->oldcapo+shift,rctx->capo+shift);CHKERRQ(ierr); 808 break; 809 case 2: 810 ierr = PetscViewerASCIIPrintf(viewer,"Store in checkpoint number %D (located in RAM)\n",rctx->check);CHKERRQ(ierr); 811 break; 812 case 3: 813 ierr = PetscViewerASCIIPrintf(viewer,"First turn: Initialize adjoints and reverse first step\n");CHKERRQ(ierr); 814 break; 815 case 4: 816 ierr = PetscViewerASCIIPrintf(viewer,"Forward and reverse one step\n");CHKERRQ(ierr); 817 break; 818 case 5: 819 ierr = PetscViewerASCIIPrintf(viewer,"Restore in checkpoint number %D (located in RAM)\n",rctx->check);CHKERRQ(ierr); 820 break; 821 case 7: 822 ierr = PetscViewerASCIIPrintf(viewer,"Store in checkpoint number %D (located on disk)\n",rctx->check);CHKERRQ(ierr); 823 break; 824 case 8: 825 ierr = PetscViewerASCIIPrintf(viewer,"Restore in checkpoint number %D (located on disk)\n",rctx->check);CHKERRQ(ierr); 826 break; 827 case -1: 828 ierr = PetscViewerASCIIPrintf(viewer,"Error!");CHKERRQ(ierr); 829 break; 830 } 831 PetscFunctionReturn(0); 832 } 833 834 static PetscErrorCode printwhattodo2(PetscViewer viewer,PetscInt whattodo,RevolveCTX *rctx,PetscInt shift) 835 { 836 PetscErrorCode ierr; 837 838 PetscFunctionBegin; 839 if (!viewer) PetscFunctionReturn(0); 840 841 switch(whattodo) { 842 case 1: 843 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Advance from stride %D to stride %D\n",rctx->oldcapo+shift,rctx->capo+shift);CHKERRQ(ierr); 844 break; 845 case 2: 846 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Store in checkpoint number %D\n",rctx->check);CHKERRQ(ierr); 847 break; 848 case 3: 849 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] First turn: Initialize adjoints and reverse first stride\n");CHKERRQ(ierr); 850 break; 851 case 4: 852 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Forward and reverse one stride\n");CHKERRQ(ierr); 853 break; 854 case 5: 855 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Restore in checkpoint number %D\n",rctx->check);CHKERRQ(ierr); 856 break; 857 case 7: 858 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Store in top-level checkpoint number %D\n",rctx->check);CHKERRQ(ierr); 859 break; 860 case 8: 861 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Restore in top-level checkpoint number %D\n",rctx->check);CHKERRQ(ierr); 862 break; 863 case -1: 864 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Error!");CHKERRQ(ierr); 865 break; 866 } 867 PetscFunctionReturn(0); 868 } 869 870 static PetscErrorCode InitRevolve(PetscInt fine,PetscInt snaps,RevolveCTX *rctx) 871 { 872 PetscFunctionBegin; 873 revolve_reset(); 874 revolve_create_offline(fine,snaps); 875 rctx->snaps_in = snaps; 876 rctx->fine = fine; 877 rctx->check = 0; 878 rctx->capo = 0; 879 rctx->reverseonestep = PETSC_FALSE; 880 /* check stepsleft? */ 881 PetscFunctionReturn(0); 882 } 883 884 static PetscErrorCode FastForwardRevolve(RevolveCTX *rctx) 885 { 886 PetscInt whattodo; 887 888 PetscFunctionBegin; 889 whattodo = 0; 890 while(whattodo!=3) { /* we have to fast forward revolve to the beginning of the backward sweep due to unfriendly revolve interface */ 891 whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 892 } 893 PetscFunctionReturn(0); 894 } 895 896 static PetscErrorCode ApplyRevolve(PetscViewer viewer,SchedulerType stype,RevolveCTX *rctx,PetscInt total_steps,PetscInt stepnum,PetscInt localstepnum,PetscBool toplevel,PetscInt *store) 897 { 898 PetscErrorCode ierr; 899 PetscInt shift,whattodo; 900 901 PetscFunctionBegin; 902 *store = 0; 903 if (rctx->stepsleft > 0) { /* advance the solution without checkpointing anything as Revolve requires */ 904 rctx->stepsleft--; 905 PetscFunctionReturn(0); 906 } 907 /* let Revolve determine what to do next */ 908 shift = stepnum-localstepnum; 909 rctx->oldcapo = rctx->capo; 910 rctx->capo = localstepnum; 911 912 if (!toplevel) whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 913 else whattodo = revolve2_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 914 if (stype == REVOLVE_ONLINE && whattodo == 8) whattodo = 5; 915 if (stype == REVOLVE_ONLINE && whattodo == 7) whattodo = 2; 916 if (!toplevel) {ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 917 else {ierr = printwhattodo2(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 918 if (whattodo == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in the Revolve library"); 919 if (whattodo == 1) { /* advance some time steps */ 920 if (stype == REVOLVE_ONLINE && rctx->capo >= total_steps-1) { 921 revolve_turn(total_steps,&rctx->capo,&rctx->fine); 922 if (!toplevel) {ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 923 else {ierr = printwhattodo2(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 924 } 925 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 926 } 927 if (whattodo == 3 || whattodo == 4) { /* ready for a reverse step */ 928 rctx->reverseonestep = PETSC_TRUE; 929 } 930 if (whattodo == 5) { /* restore a checkpoint and ask Revolve what to do next */ 931 rctx->oldcapo = rctx->capo; 932 if (!toplevel) whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); /* must return 1 or 3 or 4*/ 933 else whattodo = revolve2_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 934 if (!toplevel) {ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 935 else {ierr = printwhattodo2(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 936 if (whattodo == 3 || whattodo == 4) rctx->reverseonestep = PETSC_TRUE; 937 if (whattodo == 1) rctx->stepsleft = rctx->capo-rctx->oldcapo; 938 } 939 if (whattodo == 7) { /* save the checkpoint to disk */ 940 *store = 2; 941 rctx->oldcapo = rctx->capo; 942 whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); /* must return 1 */ 943 ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr); 944 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 945 } 946 if (whattodo == 2) { /* store a checkpoint to RAM and ask Revolve how many time steps to advance next */ 947 *store = 1; 948 rctx->oldcapo = rctx->capo; 949 if (!toplevel) whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); /* must return 1 */ 950 else whattodo = revolve2_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 951 if (!toplevel) {ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 952 else {ierr = printwhattodo2(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 953 if (stype == REVOLVE_ONLINE && rctx->capo >= total_steps-1) { 954 revolve_turn(total_steps,&rctx->capo,&rctx->fine); 955 ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr); 956 } 957 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 958 } 959 PetscFunctionReturn(0); 960 } 961 962 static PetscErrorCode SetTrajROF(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 963 { 964 Stack *stack = &tjsch->stack; 965 PetscInt store; 966 StackElement e; 967 PetscErrorCode ierr; 968 969 PetscFunctionBegin; 970 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 971 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 972 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,stepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 973 if (store == 1) { 974 if (stepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 975 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 976 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 977 ierr = StackPush(stack,e);CHKERRQ(ierr); 978 } 979 PetscFunctionReturn(0); 980 } 981 982 static PetscErrorCode GetTrajROF(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 983 { 984 Stack *stack = &tjsch->stack; 985 PetscInt whattodo,shift,store; 986 StackElement e; 987 PetscErrorCode ierr; 988 989 PetscFunctionBegin; 990 if (stepnum == 0 || stepnum == tjsch->total_steps) { 991 ierr = TurnBackward(ts);CHKERRQ(ierr); 992 tjsch->rctx->reverseonestep = PETSC_FALSE; 993 PetscFunctionReturn(0); 994 } 995 /* restore a checkpoint */ 996 ierr = StackTop(stack,&e);CHKERRQ(ierr); 997 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 998 if (stack->solution_only) { /* start with restoring a checkpoint */ 999 tjsch->rctx->capo = stepnum; 1000 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1001 shift = 0; 1002 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); 1003 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1004 } else { /* 2 revolve actions: restore a checkpoint and then advance */ 1005 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,stepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1006 if (tj->monitor) { 1007 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1008 ierr = PetscViewerASCIIPrintf(tj->monitor,"Skip the step from %D to %D (stage values already checkpointed)\n",tjsch->rctx->oldcapo,tjsch->rctx->oldcapo+1);CHKERRQ(ierr); 1009 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1010 } 1011 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1012 } 1013 if (stack->solution_only || (!stack->solution_only && e->stepnum < stepnum)) { 1014 tjsch->recompute = PETSC_TRUE; 1015 ierr = TurnForward(ts);CHKERRQ(ierr); 1016 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1017 } 1018 if ((stack->solution_only && e->stepnum+1 == stepnum) || (!stack->solution_only && e->stepnum == stepnum)) { 1019 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1020 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1021 } 1022 tjsch->rctx->reverseonestep = PETSC_FALSE; 1023 PetscFunctionReturn(0); 1024 } 1025 1026 static PetscErrorCode SetTrajRON(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 1027 { 1028 Stack *stack = &tjsch->stack; 1029 Vec *Y; 1030 PetscInt i,store; 1031 PetscReal timeprev; 1032 StackElement e; 1033 RevolveCTX *rctx = tjsch->rctx; 1034 PetscErrorCode ierr; 1035 1036 PetscFunctionBegin; 1037 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 1038 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 1039 ierr = ApplyRevolve(tj->monitor,tjsch->stype,rctx,tjsch->total_steps,stepnum,stepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1040 if (store == 1) { 1041 if (rctx->check != stack->top+1) { /* overwrite some non-top checkpoint in the stack */ 1042 ierr = StackFind(stack,&e,rctx->check);CHKERRQ(ierr); 1043 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 1044 if (stack->numY > 0 && !stack->solution_only) { 1045 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 1046 for (i=0;i<stack->numY;i++) { 1047 ierr = VecCopy(Y[i],e->Y[i]);CHKERRQ(ierr); 1048 } 1049 } 1050 e->stepnum = stepnum; 1051 e->time = time; 1052 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 1053 e->timeprev = timeprev; 1054 } else { 1055 if (stepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 1056 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1057 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 1058 ierr = StackPush(stack,e);CHKERRQ(ierr); 1059 } 1060 } 1061 PetscFunctionReturn(0); 1062 } 1063 1064 static PetscErrorCode GetTrajRON(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 1065 { 1066 Stack *stack = &tjsch->stack; 1067 PetscInt whattodo,shift; 1068 StackElement e; 1069 PetscErrorCode ierr; 1070 1071 PetscFunctionBegin; 1072 if (stepnum == 0 || stepnum == tjsch->total_steps) { 1073 ierr = TurnBackward(ts);CHKERRQ(ierr); 1074 tjsch->rctx->reverseonestep = PETSC_FALSE; 1075 PetscFunctionReturn(0); 1076 } 1077 tjsch->rctx->capo = stepnum; 1078 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1079 shift = 0; 1080 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); /* whattodo=restore */ 1081 if (whattodo == 8) whattodo = 5; 1082 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1083 /* restore a checkpoint */ 1084 ierr = StackFind(stack,&e,tjsch->rctx->check);CHKERRQ(ierr); 1085 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1086 if (!stack->solution_only) { /* whattodo must be 5 */ 1087 /* ask Revolve what to do next */ 1088 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1089 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); /* must return 1 or 3 or 4*/ 1090 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1091 if (whattodo == 3 || whattodo == 4) tjsch->rctx->reverseonestep = PETSC_TRUE; 1092 if (whattodo == 1) tjsch->rctx->stepsleft = tjsch->rctx->capo-tjsch->rctx->oldcapo; 1093 if (tj->monitor) { 1094 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1095 ierr = PetscViewerASCIIPrintf(tj->monitor,"Skip the step from %D to %D (stage values already checkpointed)\n",tjsch->rctx->oldcapo,tjsch->rctx->oldcapo+1);CHKERRQ(ierr); 1096 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1097 } 1098 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1099 } 1100 if (stack->solution_only || (!stack->solution_only && e->stepnum < stepnum)) { 1101 tjsch->recompute = PETSC_TRUE; 1102 ierr = TurnForward(ts);CHKERRQ(ierr); 1103 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1104 } 1105 tjsch->rctx->reverseonestep = PETSC_FALSE; 1106 PetscFunctionReturn(0); 1107 } 1108 1109 static PetscErrorCode SetTrajTLR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 1110 { 1111 Stack *stack = &tjsch->stack; 1112 PetscInt store,localstepnum,laststridesize; 1113 StackElement e; 1114 PetscBool done = PETSC_FALSE; 1115 PetscErrorCode ierr; 1116 1117 PetscFunctionBegin; 1118 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 1119 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 1120 1121 localstepnum = stepnum%tjsch->stride; 1122 laststridesize = tjsch->total_steps%tjsch->stride; 1123 if (!laststridesize) laststridesize = tjsch->stride; 1124 1125 if (!tjsch->recompute) { 1126 ierr = TopLevelStore(tj,ts,tjsch,stepnum,localstepnum,laststridesize,&done);CHKERRQ(ierr); 1127 /* revolve is needed for the last stride; different starting points for last stride between solutin_only and !solutin_only */ 1128 if (!stack->solution_only && !tjsch->save_stack && stepnum <= tjsch->total_steps-laststridesize) PetscFunctionReturn(0); 1129 if (stack->solution_only && !tjsch->save_stack && stepnum < tjsch->total_steps-laststridesize) PetscFunctionReturn(0); 1130 } 1131 if (tjsch->save_stack && done) { 1132 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1133 PetscFunctionReturn(0); 1134 } 1135 if (laststridesize < tjsch->stride) { 1136 if (stack->solution_only && stepnum == tjsch->total_steps-laststridesize && !tjsch->recompute) { /* step tjsch->total_steps-laststridesize-1 is skipped, but the next step is not */ 1137 ierr = InitRevolve(laststridesize,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1138 } 1139 if (!stack->solution_only && stepnum == tjsch->total_steps-laststridesize+1 && !tjsch->recompute) { /* step tjsch->total_steps-laststridesize is skipped, but the next step is not */ 1140 ierr = InitRevolve(laststridesize,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1141 } 1142 } 1143 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,localstepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1144 if (store == 1) { 1145 if (localstepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 1146 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1147 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 1148 ierr = StackPush(stack,e);CHKERRQ(ierr); 1149 } 1150 PetscFunctionReturn(0); 1151 } 1152 1153 static PetscErrorCode GetTrajTLR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 1154 { 1155 Stack *stack = &tjsch->stack; 1156 PetscInt whattodo,shift; 1157 PetscInt localstepnum,stridenum,laststridesize,store; 1158 StackElement e; 1159 PetscErrorCode ierr; 1160 1161 PetscFunctionBegin; 1162 localstepnum = stepnum%tjsch->stride; 1163 stridenum = stepnum/tjsch->stride; 1164 if (stepnum == tjsch->total_steps) { 1165 ierr = TurnBackward(ts);CHKERRQ(ierr); 1166 tjsch->rctx->reverseonestep = PETSC_FALSE; 1167 PetscFunctionReturn(0); 1168 } 1169 laststridesize = tjsch->total_steps%tjsch->stride; 1170 if (!laststridesize) laststridesize = tjsch->stride; 1171 if (stack->solution_only) { 1172 /* fill stack */ 1173 if (localstepnum == 0 && stepnum <= tjsch->total_steps-laststridesize) { 1174 if (tjsch->save_stack) { 1175 ierr = StackLoadAll(tj,ts,stack,stridenum);CHKERRQ(ierr); 1176 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1177 ierr = FastForwardRevolve(tjsch->rctx);CHKERRQ(ierr); 1178 tjsch->recompute = PETSC_TRUE; 1179 tjsch->skip_trajectory = PETSC_TRUE; 1180 ierr = TurnForward(ts);CHKERRQ(ierr); 1181 ierr = ReCompute(ts,tjsch,stridenum*tjsch->stride-1,stridenum*tjsch->stride);CHKERRQ(ierr); 1182 tjsch->skip_trajectory = PETSC_FALSE; 1183 } else { 1184 ierr = LoadSingle(tj,ts,stack,stridenum);CHKERRQ(ierr); 1185 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1186 tjsch->recompute = PETSC_TRUE; 1187 ierr = TurnForward(ts);CHKERRQ(ierr); 1188 ierr = ReCompute(ts,tjsch,(stridenum-1)*tjsch->stride,stridenum*tjsch->stride);CHKERRQ(ierr); 1189 } 1190 PetscFunctionReturn(0); 1191 } 1192 /* restore a checkpoint */ 1193 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1194 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1195 /* start with restoring a checkpoint */ 1196 tjsch->rctx->capo = stepnum; 1197 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1198 shift = stepnum-localstepnum; 1199 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); 1200 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1201 tjsch->recompute = PETSC_TRUE; 1202 ierr = TurnForward(ts);CHKERRQ(ierr); 1203 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1204 if (e->stepnum+1 == stepnum) { 1205 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1206 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1207 } 1208 } else { 1209 /* fill stack with info */ 1210 if (localstepnum == 0 && tjsch->total_steps-stepnum >= laststridesize) { 1211 if (tjsch->save_stack) { 1212 ierr = StackLoadAll(tj,ts,stack,stridenum);CHKERRQ(ierr); 1213 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1214 ierr = FastForwardRevolve(tjsch->rctx);CHKERRQ(ierr); 1215 } else { 1216 ierr = LoadSingle(tj,ts,stack,stridenum);CHKERRQ(ierr); 1217 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1218 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,(stridenum-1)*tjsch->stride+1,1,PETSC_FALSE,&store);CHKERRQ(ierr); 1219 if (tj->monitor) { 1220 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1221 ierr = PetscViewerASCIIPrintf(tj->monitor,"Skip the step from %D to %D (stage values already checkpointed)\n",(stridenum-1)*tjsch->stride+tjsch->rctx->oldcapo,(stridenum-1)*tjsch->stride+tjsch->rctx->oldcapo+1);CHKERRQ(ierr); 1222 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1223 } 1224 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1225 ierr = ElementSet(ts,stack,&e,(stridenum-1)*tjsch->stride+1,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1226 ierr = StackPush(stack,e);CHKERRQ(ierr); 1227 tjsch->recompute = PETSC_TRUE; 1228 ierr = TurnForward(ts);CHKERRQ(ierr); 1229 ierr = ReCompute(ts,tjsch,e->stepnum,stridenum*tjsch->stride);CHKERRQ(ierr); 1230 } 1231 PetscFunctionReturn(0); 1232 } 1233 /* restore a checkpoint */ 1234 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1235 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1236 /* 2 revolve actions: restore a checkpoint and then advance */ 1237 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,localstepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1238 if (tj->monitor) { 1239 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1240 ierr = PetscViewerASCIIPrintf(tj->monitor,"Skip the step from %D to %D (stage values already checkpointed)\n",stepnum-localstepnum+tjsch->rctx->oldcapo,stepnum-localstepnum+tjsch->rctx->oldcapo+1);CHKERRQ(ierr); 1241 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1242 } 1243 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1244 if (e->stepnum < stepnum) { 1245 tjsch->recompute = PETSC_TRUE; 1246 ierr = TurnForward(ts);CHKERRQ(ierr); 1247 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1248 } 1249 if (e->stepnum == stepnum) { 1250 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1251 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1252 } 1253 } 1254 tjsch->rctx->reverseonestep = PETSC_FALSE; 1255 PetscFunctionReturn(0); 1256 } 1257 1258 static PetscErrorCode SetTrajTLTR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 1259 { 1260 Stack *stack = &tjsch->stack; 1261 PetscInt store,localstepnum,stridenum,laststridesize; 1262 StackElement e; 1263 PetscBool done = PETSC_FALSE; 1264 PetscErrorCode ierr; 1265 1266 PetscFunctionBegin; 1267 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 1268 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 1269 1270 localstepnum = stepnum%tjsch->stride; /* index at the bottom level (inside a stride) */ 1271 stridenum = stepnum/tjsch->stride; /* index at the top level */ 1272 laststridesize = tjsch->total_steps%tjsch->stride; 1273 if (!laststridesize) laststridesize = tjsch->stride; 1274 if (stack->solution_only && localstepnum == 0 && !tjsch->rctx2->reverseonestep) { 1275 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx2,(tjsch->total_steps+tjsch->stride-1)/tjsch->stride,stridenum,stridenum,PETSC_TRUE,&tjsch->store_stride);CHKERRQ(ierr); 1276 if (laststridesize < tjsch->stride && stepnum == tjsch->total_steps-laststridesize) { 1277 ierr = InitRevolve(laststridesize,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1278 } 1279 } 1280 if (!stack->solution_only && localstepnum == 1 && !tjsch->rctx2->reverseonestep) { 1281 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx2,(tjsch->total_steps+tjsch->stride-1)/tjsch->stride,stridenum,stridenum,PETSC_TRUE,&tjsch->store_stride);CHKERRQ(ierr); 1282 if (laststridesize < tjsch->stride && stepnum == tjsch->total_steps-laststridesize+1) { 1283 ierr = InitRevolve(laststridesize,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1284 } 1285 } 1286 if (tjsch->store_stride) { 1287 ierr = TopLevelStore(tj,ts,tjsch,stepnum,localstepnum,laststridesize,&done);CHKERRQ(ierr); 1288 if (done) { 1289 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1290 PetscFunctionReturn(0); 1291 } 1292 } 1293 if (stepnum < tjsch->total_steps-laststridesize) { 1294 if (tjsch->save_stack && !tjsch->store_stride && !tjsch->rctx2->reverseonestep) PetscFunctionReturn(0); /* store or forward-and-reverse at top level trigger revolve at bottom level */ 1295 if (!tjsch->save_stack && !tjsch->rctx2->reverseonestep) PetscFunctionReturn(0); /* store operation does not require revolve be called at bottom level */ 1296 } 1297 /* Skipping stepnum=0 for !stack->only is enough for TLR, but not for TLTR. Here we skip the first step for each stride so that the top-level revolve is applied (always at localstepnum=1) ahead of the bottom-level revolve */ 1298 if (!stack->solution_only && localstepnum == 0 && stepnum != tjsch->total_steps && !tjsch->recompute) PetscFunctionReturn(0); 1299 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,localstepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1300 if (store == 1) { 1301 if (localstepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 1302 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1303 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 1304 ierr = StackPush(stack,e);CHKERRQ(ierr); 1305 } 1306 PetscFunctionReturn(0); 1307 } 1308 1309 static PetscErrorCode GetTrajTLTR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 1310 { 1311 Stack *stack = &tjsch->stack; 1312 DiskStack *diskstack = &tjsch->diskstack; 1313 PetscInt whattodo,shift; 1314 PetscInt localstepnum,stridenum,restoredstridenum,laststridesize,store; 1315 StackElement e; 1316 PetscErrorCode ierr; 1317 1318 PetscFunctionBegin; 1319 localstepnum = stepnum%tjsch->stride; 1320 stridenum = stepnum/tjsch->stride; 1321 if (stepnum == tjsch->total_steps) { 1322 ierr = TurnBackward(ts);CHKERRQ(ierr); 1323 tjsch->rctx->reverseonestep = PETSC_FALSE; 1324 PetscFunctionReturn(0); 1325 } 1326 laststridesize = tjsch->total_steps%tjsch->stride; 1327 if (!laststridesize) laststridesize = tjsch->stride; 1328 /* 1329 Last stride can be adjoined directly. All the other strides require that the stack in memory be ready before an adjoint step is taken (at the end of each stride). The following two cases need to be addressed differently: 1330 Case 1 (save_stack) 1331 Restore a disk checkpoint; update TS with the last element in the restored data; recompute to the current point. 1332 Case 2 (!save_stack) 1333 Restore a disk checkpoint; update TS with the restored point; recompute to the current point. 1334 */ 1335 if (localstepnum == 0 && stepnum <= tjsch->total_steps-laststridesize) { 1336 /* restore the top element in the stack for disk checkpoints */ 1337 restoredstridenum = diskstack->container[diskstack->top]; 1338 tjsch->rctx2->reverseonestep = PETSC_FALSE; 1339 /* top-level revolve must be applied before current step, just like the solution_only mode for single-level revolve */ 1340 if (!tjsch->save_stack && stack->solution_only) { /* start with restoring a checkpoint */ 1341 tjsch->rctx2->capo = stridenum; 1342 tjsch->rctx2->oldcapo = tjsch->rctx2->capo; 1343 shift = 0; 1344 whattodo = revolve2_action(&tjsch->rctx2->check,&tjsch->rctx2->capo,&tjsch->rctx2->fine,tjsch->rctx2->snaps_in,&tjsch->rctx2->info,&tjsch->rctx2->where); 1345 ierr = printwhattodo2(tj->monitor,whattodo,tjsch->rctx2,shift);CHKERRQ(ierr); 1346 } else { /* 2 revolve actions: restore a checkpoint and then advance */ 1347 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx2,(tjsch->total_steps+tjsch->stride-1)/tjsch->stride,stridenum,stridenum,PETSC_TRUE,&tjsch->store_stride);CHKERRQ(ierr); 1348 if (tj->monitor) { 1349 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1350 ierr = PetscViewerASCIIPrintf(tj->monitor,"[Top Level] Skip the stride from %D to %D (stage values already checkpointed)\n",tjsch->rctx2->oldcapo,tjsch->rctx2->oldcapo+1);CHKERRQ(ierr); 1351 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1352 } 1353 if (!tjsch->rctx2->reverseonestep && tjsch->rctx2->stepsleft > 0) tjsch->rctx2->stepsleft--; 1354 } 1355 /* fill stack */ 1356 if (stack->solution_only) { 1357 if (tjsch->save_stack) { 1358 if (restoredstridenum < stridenum) { 1359 ierr = StackLoadLast(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1360 } else { 1361 ierr = StackLoadAll(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1362 } 1363 /* recompute one step ahead */ 1364 tjsch->recompute = PETSC_TRUE; 1365 tjsch->skip_trajectory = PETSC_TRUE; 1366 ierr = TurnForward(ts);CHKERRQ(ierr); 1367 ierr = ReCompute(ts,tjsch,stridenum*tjsch->stride-1,stridenum*tjsch->stride);CHKERRQ(ierr); 1368 tjsch->skip_trajectory = PETSC_FALSE; 1369 if (restoredstridenum < stridenum) { 1370 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1371 tjsch->recompute = PETSC_TRUE; 1372 ierr = TurnForward(ts);CHKERRQ(ierr); 1373 ierr = ReCompute(ts,tjsch,restoredstridenum*tjsch->stride,stepnum);CHKERRQ(ierr); 1374 } else { /* stack ready, fast forward revolve status */ 1375 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1376 ierr = FastForwardRevolve(tjsch->rctx);CHKERRQ(ierr); 1377 } 1378 } else { 1379 ierr = LoadSingle(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1380 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1381 tjsch->recompute = PETSC_TRUE; 1382 ierr = TurnForward(ts);CHKERRQ(ierr); 1383 ierr = ReCompute(ts,tjsch,(restoredstridenum-1)*tjsch->stride,stepnum);CHKERRQ(ierr); 1384 } 1385 } else { 1386 if (tjsch->save_stack) { 1387 if (restoredstridenum < stridenum) { 1388 ierr = StackLoadLast(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1389 /* reset revolve */ 1390 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1391 tjsch->recompute = PETSC_TRUE; 1392 ierr = TurnForward(ts);CHKERRQ(ierr); 1393 ierr = ReCompute(ts,tjsch,restoredstridenum*tjsch->stride,stepnum);CHKERRQ(ierr); 1394 } else { /* stack ready, fast forward revolve status */ 1395 ierr = StackLoadAll(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1396 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1397 ierr = FastForwardRevolve(tjsch->rctx);CHKERRQ(ierr); 1398 } 1399 } else { 1400 ierr = LoadSingle(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1401 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1402 /* push first element to stack */ 1403 if (tjsch->store_stride || tjsch->rctx2->reverseonestep) { 1404 shift = (restoredstridenum-1)*tjsch->stride-localstepnum; 1405 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,(restoredstridenum-1)*tjsch->stride+1,1,PETSC_FALSE,&store);CHKERRQ(ierr); 1406 if (tj->monitor) { 1407 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1408 ierr = PetscViewerASCIIPrintf(tj->monitor,"Skip the step from %D to %D (stage values already checkpointed)\n",(restoredstridenum-1)*tjsch->stride,(restoredstridenum-1)*tjsch->stride+1);CHKERRQ(ierr); 1409 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1410 } 1411 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1412 ierr = ElementSet(ts,stack,&e,(restoredstridenum-1)*tjsch->stride+1,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1413 ierr = StackPush(stack,e);CHKERRQ(ierr); 1414 } 1415 tjsch->recompute = PETSC_TRUE; 1416 ierr = TurnForward(ts);CHKERRQ(ierr); 1417 ierr = ReCompute(ts,tjsch,(restoredstridenum-1)*tjsch->stride+1,stepnum);CHKERRQ(ierr); 1418 } 1419 } 1420 if (restoredstridenum == stridenum) diskstack->top--; 1421 tjsch->rctx->reverseonestep = PETSC_FALSE; 1422 PetscFunctionReturn(0); 1423 } 1424 1425 if (stack->solution_only) { 1426 /* restore a checkpoint */ 1427 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1428 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1429 /* start with restoring a checkpoint */ 1430 tjsch->rctx->capo = stepnum; 1431 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1432 shift = stepnum-localstepnum; 1433 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); 1434 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1435 tjsch->recompute = PETSC_TRUE; 1436 ierr = TurnForward(ts);CHKERRQ(ierr); 1437 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1438 if (e->stepnum+1 == stepnum) { 1439 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1440 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1441 } 1442 } else { 1443 /* restore a checkpoint */ 1444 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1445 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1446 /* 2 revolve actions: restore a checkpoint and then advance */ 1447 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,localstepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1448 if (tj->monitor) { 1449 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1450 ierr = PetscViewerASCIIPrintf(tj->monitor,"Skip the step from %D to %D (stage values already checkpointed)\n",stepnum-localstepnum+tjsch->rctx->oldcapo,stepnum-localstepnum+tjsch->rctx->oldcapo+1);CHKERRQ(ierr); 1451 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1452 } 1453 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1454 if (e->stepnum < stepnum) { 1455 tjsch->recompute = PETSC_TRUE; 1456 ierr = TurnForward(ts);CHKERRQ(ierr); 1457 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1458 } 1459 if (e->stepnum == stepnum) { 1460 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1461 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1462 } 1463 } 1464 tjsch->rctx->reverseonestep = PETSC_FALSE; 1465 PetscFunctionReturn(0); 1466 } 1467 1468 static PetscErrorCode SetTrajRMS(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 1469 { 1470 Stack *stack = &tjsch->stack; 1471 PetscInt store; 1472 StackElement e; 1473 PetscErrorCode ierr; 1474 1475 PetscFunctionBegin; 1476 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 1477 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 1478 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,stepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1479 if (store == 1){ 1480 if (stepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 1481 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1482 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 1483 ierr = StackPush(stack,e);CHKERRQ(ierr); 1484 } else if (store == 2) { 1485 ierr = DumpSingle(tj,ts,stack,tjsch->rctx->check+1);CHKERRQ(ierr); 1486 } 1487 PetscFunctionReturn(0); 1488 } 1489 1490 static PetscErrorCode GetTrajRMS(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 1491 { 1492 Stack *stack = &tjsch->stack; 1493 PetscInt whattodo,shift; 1494 PetscInt restart; 1495 PetscBool ondisk; 1496 StackElement e; 1497 PetscErrorCode ierr; 1498 1499 PetscFunctionBegin; 1500 if (stepnum == 0 || stepnum == tjsch->total_steps) { 1501 ierr = TurnBackward(ts);CHKERRQ(ierr); 1502 tjsch->rctx->reverseonestep = PETSC_FALSE; 1503 PetscFunctionReturn(0); 1504 } 1505 tjsch->rctx->capo = stepnum; 1506 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1507 shift = 0; 1508 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); /* whattodo=restore */ 1509 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1510 /* restore a checkpoint */ 1511 restart = tjsch->rctx->capo; 1512 if (!tjsch->rctx->where) { 1513 ondisk = PETSC_TRUE; 1514 ierr = LoadSingle(tj,ts,stack,tjsch->rctx->check+1);CHKERRQ(ierr); 1515 ierr = TurnBackward(ts);CHKERRQ(ierr); 1516 } else { 1517 ondisk = PETSC_FALSE; 1518 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1519 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1520 } 1521 if (!stack->solution_only) { /* whattodo must be 5 or 8 */ 1522 /* ask Revolve what to do next */ 1523 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1524 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); /* must return 1 or 3 or 4*/ 1525 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1526 if (whattodo == 3 || whattodo == 4) tjsch->rctx->reverseonestep = PETSC_TRUE; 1527 if (whattodo == 1) tjsch->rctx->stepsleft = tjsch->rctx->capo-tjsch->rctx->oldcapo; 1528 if (tj->monitor) { 1529 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1530 ierr = PetscViewerASCIIPrintf(tj->monitor,"Skip the step from %D to %D (stage values already checkpointed)\n",tjsch->rctx->oldcapo,tjsch->rctx->oldcapo+1);CHKERRQ(ierr); 1531 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1532 } 1533 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1534 restart++; /* skip one step */ 1535 } 1536 if (stack->solution_only || (!stack->solution_only && restart < stepnum)) { 1537 tjsch->recompute = PETSC_TRUE; 1538 ierr = TurnForward(ts);CHKERRQ(ierr); 1539 ierr = ReCompute(ts,tjsch,restart,stepnum);CHKERRQ(ierr); 1540 } 1541 if (!ondisk && ( (stack->solution_only && e->stepnum+1 == stepnum) || (!stack->solution_only && e->stepnum == stepnum) )) { 1542 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1543 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1544 } 1545 tjsch->rctx->reverseonestep = PETSC_FALSE; 1546 PetscFunctionReturn(0); 1547 } 1548 #endif 1549 1550 static PetscErrorCode TSTrajectorySet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X) 1551 { 1552 TJScheduler *tjsch = (TJScheduler*)tj->data; 1553 PetscErrorCode ierr; 1554 1555 PetscFunctionBegin; 1556 if (!tjsch->recompute) { /* use global stepnum in the forward sweep */ 1557 ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr); 1558 } 1559 /* for consistency */ 1560 if (!tjsch->recompute && stepnum == 0) ts->ptime_prev = ts->ptime-ts->time_step; 1561 switch (tjsch->stype) { 1562 case NONE: 1563 if (tj->adjoint_solve_mode) { 1564 ierr = SetTrajN(ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1565 } else { 1566 ierr = SetTrajN_2(ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1567 } 1568 break; 1569 case TWO_LEVEL_NOREVOLVE: 1570 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1571 ierr = SetTrajTLNR(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1572 break; 1573 #if defined(PETSC_HAVE_REVOLVE) 1574 case TWO_LEVEL_REVOLVE: 1575 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1576 ierr = SetTrajTLR(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1577 break; 1578 case TWO_LEVEL_TWO_REVOLVE: 1579 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1580 ierr = SetTrajTLTR(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1581 break; 1582 case REVOLVE_OFFLINE: 1583 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1584 ierr = SetTrajROF(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1585 break; 1586 case REVOLVE_ONLINE: 1587 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1588 ierr = SetTrajRON(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1589 break; 1590 case REVOLVE_MULTISTAGE: 1591 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1592 ierr = SetTrajRMS(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1593 break; 1594 #endif 1595 default: 1596 break; 1597 } 1598 PetscFunctionReturn(0); 1599 } 1600 1601 static PetscErrorCode TSTrajectoryGet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t) 1602 { 1603 TJScheduler *tjsch = (TJScheduler*)tj->data; 1604 PetscErrorCode ierr; 1605 1606 PetscFunctionBegin; 1607 if (tj->adjoint_solve_mode && stepnum == 0) PetscFunctionReturn(0); 1608 switch (tjsch->stype) { 1609 case NONE: 1610 if (tj->adjoint_solve_mode) { 1611 ierr = GetTrajN(ts,tjsch,stepnum);CHKERRQ(ierr); 1612 } else { 1613 ierr = GetTrajN_2(ts,tjsch,stepnum);CHKERRQ(ierr); 1614 } 1615 break; 1616 case TWO_LEVEL_NOREVOLVE: 1617 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1618 ierr = GetTrajTLNR(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1619 break; 1620 #if defined(PETSC_HAVE_REVOLVE) 1621 case TWO_LEVEL_REVOLVE: 1622 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1623 ierr = GetTrajTLR(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1624 break; 1625 case TWO_LEVEL_TWO_REVOLVE: 1626 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1627 ierr = GetTrajTLTR(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1628 break; 1629 case REVOLVE_OFFLINE: 1630 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1631 ierr = GetTrajROF(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1632 break; 1633 case REVOLVE_ONLINE: 1634 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1635 ierr = GetTrajRON(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1636 break; 1637 case REVOLVE_MULTISTAGE: 1638 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1639 ierr = GetTrajRMS(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1640 break; 1641 #endif 1642 default: 1643 break; 1644 } 1645 PetscFunctionReturn(0); 1646 } 1647 1648 PETSC_UNUSED static PetscErrorCode TSTrajectorySetStride_Memory(TSTrajectory tj,PetscInt stride) 1649 { 1650 TJScheduler *tjsch = (TJScheduler*)tj->data; 1651 1652 PetscFunctionBegin; 1653 tjsch->stride = stride; 1654 PetscFunctionReturn(0); 1655 } 1656 1657 PETSC_UNUSED static PetscErrorCode TSTrajectorySetMaxCpsRAM_Memory(TSTrajectory tj,PetscInt max_cps_ram) 1658 { 1659 TJScheduler *tjsch = (TJScheduler*)tj->data; 1660 1661 PetscFunctionBegin; 1662 tjsch->max_cps_ram = max_cps_ram; 1663 PetscFunctionReturn(0); 1664 } 1665 1666 PETSC_UNUSED static PetscErrorCode TSTrajectorySetMaxCpsDisk_Memory(TSTrajectory tj,PetscInt max_cps_disk) 1667 { 1668 TJScheduler *tjsch = (TJScheduler*)tj->data; 1669 1670 PetscFunctionBegin; 1671 tjsch->max_cps_disk = max_cps_disk; 1672 PetscFunctionReturn(0); 1673 } 1674 1675 #if defined(PETSC_HAVE_REVOLVE) 1676 PETSC_UNUSED static PetscErrorCode TSTrajectorySetRevolveOnline(TSTrajectory tj,PetscBool use_online) 1677 { 1678 TJScheduler *tjsch = (TJScheduler*)tj->data; 1679 1680 PetscFunctionBegin; 1681 tjsch->use_online = use_online; 1682 PetscFunctionReturn(0); 1683 } 1684 #endif 1685 1686 PETSC_UNUSED static PetscErrorCode TSTrajectorySetSaveStack(TSTrajectory tj,PetscBool save_stack) 1687 { 1688 TJScheduler *tjsch = (TJScheduler*)tj->data; 1689 1690 PetscFunctionBegin; 1691 tjsch->save_stack = save_stack; 1692 PetscFunctionReturn(0); 1693 } 1694 1695 PETSC_UNUSED static PetscErrorCode TSTrajectorySetUseDRAM(TSTrajectory tj,PetscBool use_dram) 1696 { 1697 TJScheduler *tjsch = (TJScheduler*)tj->data; 1698 1699 PetscFunctionBegin; 1700 tjsch->stack.use_dram = use_dram; 1701 PetscFunctionReturn(0); 1702 } 1703 1704 static PetscErrorCode TSTrajectorySetFromOptions_Memory(PetscOptionItems *PetscOptionsObject,TSTrajectory tj) 1705 { 1706 TJScheduler *tjsch = (TJScheduler*)tj->data; 1707 PetscErrorCode ierr; 1708 1709 PetscFunctionBegin; 1710 ierr = PetscOptionsHead(PetscOptionsObject,"Memory based TS trajectory options");CHKERRQ(ierr); 1711 { 1712 ierr = PetscOptionsInt("-ts_trajectory_max_cps_ram","Maximum number of checkpoints in RAM","TSTrajectorySetMaxCpsRAM_Memory",tjsch->max_cps_ram,&tjsch->max_cps_ram,NULL);CHKERRQ(ierr); 1713 ierr = PetscOptionsInt("-ts_trajectory_max_cps_disk","Maximum number of checkpoints on disk","TSTrajectorySetMaxCpsDisk_Memory",tjsch->max_cps_disk,&tjsch->max_cps_disk,NULL);CHKERRQ(ierr); 1714 ierr = PetscOptionsInt("-ts_trajectory_stride","Stride to save checkpoints to file","TSTrajectorySetStride_Memory",tjsch->stride,&tjsch->stride,NULL);CHKERRQ(ierr); 1715 #if defined(PETSC_HAVE_REVOLVE) 1716 ierr = PetscOptionsBool("-ts_trajectory_revolve_online","Trick TS trajectory into using online mode of revolve","TSTrajectorySetRevolveOnline",tjsch->use_online,&tjsch->use_online,NULL);CHKERRQ(ierr); 1717 #endif 1718 ierr = PetscOptionsBool("-ts_trajectory_save_stack","Save all stack to disk","TSTrajectorySetSaveStack",tjsch->save_stack,&tjsch->save_stack,NULL);CHKERRQ(ierr); 1719 ierr = PetscOptionsBool("-ts_trajectory_use_dram","Use DRAM for checkpointing","TSTrajectorySetUseDRAM",tjsch->stack.use_dram,&tjsch->stack.use_dram,NULL);CHKERRQ(ierr); 1720 } 1721 ierr = PetscOptionsTail();CHKERRQ(ierr); 1722 tjsch->stack.solution_only = tj->solution_only; 1723 PetscFunctionReturn(0); 1724 } 1725 1726 static PetscErrorCode TSTrajectorySetUp_Memory(TSTrajectory tj,TS ts) 1727 { 1728 TJScheduler *tjsch = (TJScheduler*)tj->data; 1729 Stack *stack = &tjsch->stack; 1730 #if defined(PETSC_HAVE_REVOLVE) 1731 RevolveCTX *rctx,*rctx2; 1732 DiskStack *diskstack = &tjsch->diskstack; 1733 PetscInt diskblocks; 1734 #endif 1735 PetscInt numY,total_steps; 1736 PetscBool fixedtimestep; 1737 PetscErrorCode ierr; 1738 1739 PetscFunctionBegin; 1740 if (ts->adapt) { 1741 ierr = PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&fixedtimestep);CHKERRQ(ierr); 1742 } else fixedtimestep = PETSC_TRUE; 1743 total_steps = (PetscInt)(PetscCeilReal((ts->max_time-ts->ptime)/ts->time_step)); 1744 total_steps = total_steps < 0 ? PETSC_MAX_INT : total_steps; 1745 if (fixedtimestep) tjsch->total_steps = PetscMin(ts->max_steps,total_steps); 1746 if (tjsch->max_cps_ram > 0) stack->stacksize = tjsch->max_cps_ram; 1747 1748 if (tjsch->stride > 1) { /* two level mode */ 1749 if (tjsch->save_stack && tjsch->max_cps_disk > 1 && tjsch->max_cps_disk <= tjsch->max_cps_ram) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"The specified disk capacity is not enough to store a full stack of RAM checkpoints. You might want to change the disk capacity or use single level checkpointing instead."); 1750 if (tjsch->max_cps_disk <= 1 && tjsch->max_cps_ram > 1 && tjsch->max_cps_ram <= tjsch->stride-1) tjsch->stype = TWO_LEVEL_REVOLVE; /* use revolve_offline for each stride */ 1751 if (tjsch->max_cps_disk > 1 && tjsch->max_cps_ram > 1 && tjsch->max_cps_ram <= tjsch->stride-1) tjsch->stype = TWO_LEVEL_TWO_REVOLVE; /* use revolve_offline for each stride */ 1752 if (tjsch->max_cps_disk <= 1 && (tjsch->max_cps_ram >= tjsch->stride || tjsch->max_cps_ram == -1)) tjsch->stype = TWO_LEVEL_NOREVOLVE; /* can also be handled by TWO_LEVEL_REVOLVE */ 1753 } else { /* single level mode */ 1754 if (fixedtimestep) { 1755 if (tjsch->max_cps_ram >= tjsch->total_steps-1 || tjsch->max_cps_ram < 1) tjsch->stype = NONE; /* checkpoint all */ 1756 else tjsch->stype = (tjsch->max_cps_disk>1) ? REVOLVE_MULTISTAGE : REVOLVE_OFFLINE; 1757 } else tjsch->stype = NONE; /* checkpoint all for adaptive time step */ 1758 #if defined(PETSC_HAVE_REVOLVE) 1759 if (tjsch->use_online) tjsch->stype = REVOLVE_ONLINE; /* trick into online (for testing purpose only) */ 1760 #endif 1761 } 1762 1763 if (tjsch->stype > TWO_LEVEL_NOREVOLVE) { 1764 #ifndef PETSC_HAVE_REVOLVE 1765 SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"revolve is needed when there is not enough memory to checkpoint all time steps according to the user's settings, please reconfigure with the additional option --download-revolve."); 1766 #else 1767 switch (tjsch->stype) { 1768 case TWO_LEVEL_REVOLVE: 1769 revolve_create_offline(tjsch->stride,tjsch->max_cps_ram); 1770 break; 1771 case TWO_LEVEL_TWO_REVOLVE: 1772 diskblocks = tjsch->save_stack ? tjsch->max_cps_disk/(tjsch->max_cps_ram+1) : tjsch->max_cps_disk; /* The block size depends on whether the stack is saved. */ 1773 diskstack->stacksize = diskblocks; 1774 revolve_create_offline(tjsch->stride,tjsch->max_cps_ram); 1775 revolve2_create_offline((tjsch->total_steps+tjsch->stride-1)/tjsch->stride,diskblocks); 1776 ierr = PetscCalloc1(1,&rctx2);CHKERRQ(ierr); 1777 rctx2->snaps_in = diskblocks; 1778 rctx2->reverseonestep = PETSC_FALSE; 1779 rctx2->check = 0; 1780 rctx2->oldcapo = 0; 1781 rctx2->capo = 0; 1782 rctx2->info = 2; 1783 rctx2->fine = (tjsch->total_steps+tjsch->stride-1)/tjsch->stride; 1784 tjsch->rctx2 = rctx2; 1785 diskstack->top = -1; 1786 ierr = PetscMalloc1(diskstack->stacksize*sizeof(PetscInt),&diskstack->container);CHKERRQ(ierr); 1787 break; 1788 case REVOLVE_OFFLINE: 1789 revolve_create_offline(tjsch->total_steps,tjsch->max_cps_ram); 1790 break; 1791 case REVOLVE_ONLINE: 1792 stack->stacksize = tjsch->max_cps_ram; 1793 revolve_create_online(tjsch->max_cps_ram); 1794 break; 1795 case REVOLVE_MULTISTAGE: 1796 revolve_create_multistage(tjsch->total_steps,tjsch->max_cps_ram+tjsch->max_cps_disk,tjsch->max_cps_ram); 1797 break; 1798 default: 1799 break; 1800 } 1801 ierr = PetscCalloc1(1,&rctx);CHKERRQ(ierr); 1802 rctx->snaps_in = tjsch->max_cps_ram; /* for theta methods snaps_in=2*max_cps_ram */ 1803 rctx->reverseonestep = PETSC_FALSE; 1804 rctx->check = 0; 1805 rctx->oldcapo = 0; 1806 rctx->capo = 0; 1807 rctx->info = 2; 1808 rctx->fine = (tjsch->stride > 1) ? tjsch->stride : tjsch->total_steps; 1809 tjsch->rctx = rctx; 1810 if (tjsch->stype == REVOLVE_ONLINE) rctx->fine = -1; 1811 #endif 1812 } else { 1813 if (tjsch->stype == TWO_LEVEL_NOREVOLVE) stack->stacksize = tjsch->stride-1; /* need tjsch->stride-1 at most */ 1814 if (tjsch->stype == NONE) { 1815 if (fixedtimestep) stack->stacksize = stack->solution_only ? tjsch->total_steps : tjsch->total_steps-1; 1816 else { /* adaptive time step */ 1817 /* if max_cps_ram is not specified, use maximal allowed number of steps for stack size */ 1818 if(tjsch->max_cps_ram == -1) stack->stacksize = ts->max_steps < PETSC_MAX_INT ? ts->max_steps : 10000; 1819 tjsch->total_steps = stack->solution_only ? stack->stacksize : stack->stacksize+1; /* will be updated as time integration advances */ 1820 } 1821 } 1822 } 1823 1824 tjsch->recompute = PETSC_FALSE; 1825 ierr = TSGetStages(ts,&numY,PETSC_IGNORE);CHKERRQ(ierr); 1826 ierr = StackCreate(stack,stack->stacksize,numY);CHKERRQ(ierr); 1827 PetscFunctionReturn(0); 1828 } 1829 1830 static PetscErrorCode TSTrajectoryReset_Memory(TSTrajectory tj) 1831 { 1832 TJScheduler *tjsch = (TJScheduler*)tj->data; 1833 PetscErrorCode ierr; 1834 1835 PetscFunctionBegin; 1836 if (tjsch->stype > TWO_LEVEL_NOREVOLVE) { 1837 #if defined(PETSC_HAVE_REVOLVE) 1838 revolve_reset(); 1839 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) { 1840 revolve2_reset(); 1841 ierr = PetscFree(tjsch->diskstack.container);CHKERRQ(ierr); 1842 } 1843 #endif 1844 } 1845 ierr = StackDestroy(&tjsch->stack);CHKERRQ(ierr); 1846 #if defined(PETSC_HAVE_REVOLVE) 1847 if (tjsch->stype > TWO_LEVEL_NOREVOLVE) { 1848 ierr = PetscFree(tjsch->rctx);CHKERRQ(ierr); 1849 ierr = PetscFree(tjsch->rctx2);CHKERRQ(ierr); 1850 } 1851 #endif 1852 PetscFunctionReturn(0); 1853 } 1854 1855 static PetscErrorCode TSTrajectoryDestroy_Memory(TSTrajectory tj) 1856 { 1857 TJScheduler *tjsch = (TJScheduler*)tj->data; 1858 PetscErrorCode ierr; 1859 1860 PetscFunctionBegin; 1861 ierr = PetscFree(tjsch);CHKERRQ(ierr); 1862 PetscFunctionReturn(0); 1863 } 1864 1865 /*MC 1866 TSTRAJECTORYMEMORY - Stores each solution of the ODE/ADE in memory 1867 1868 Level: intermediate 1869 1870 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType() 1871 1872 M*/ 1873 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory tj,TS ts) 1874 { 1875 TJScheduler *tjsch; 1876 PetscErrorCode ierr; 1877 1878 PetscFunctionBegin; 1879 tj->ops->set = TSTrajectorySet_Memory; 1880 tj->ops->get = TSTrajectoryGet_Memory; 1881 tj->ops->setup = TSTrajectorySetUp_Memory; 1882 tj->ops->setfromoptions = TSTrajectorySetFromOptions_Memory; 1883 tj->ops->reset = TSTrajectoryReset_Memory; 1884 tj->ops->destroy = TSTrajectoryDestroy_Memory; 1885 1886 ierr = PetscCalloc1(1,&tjsch);CHKERRQ(ierr); 1887 tjsch->stype = NONE; 1888 tjsch->max_cps_ram = -1; /* -1 indicates that it is not set */ 1889 tjsch->max_cps_disk = -1; /* -1 indicates that it is not set */ 1890 tjsch->stride = 0; /* if not zero, two-level checkpointing will be used */ 1891 #if defined(PETSC_HAVE_REVOLVE) 1892 tjsch->use_online = PETSC_FALSE; 1893 #endif 1894 tjsch->save_stack = PETSC_TRUE; 1895 1896 tjsch->stack.solution_only = tj->solution_only; 1897 1898 tj->data = tjsch; 1899 PetscFunctionReturn(0); 1900 } 1901