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 tjsch->recompute = PETSC_TRUE; /* hints TSTrajectorySet() that it is in recompute mode */ 528 ierr = TSSetStepNumber(ts,stepnumbegin);CHKERRQ(ierr);/* global step number */ 529 for (i=stepnumbegin;i<stepnumend;i++) { /* assume fixed step size */ 530 if (stack->solution_only && !tjsch->skip_trajectory) { /* revolve online need this */ 531 /* don't use the public interface as it will update the TSHistory: this need a better fix */ 532 ierr = TSTrajectorySet_Memory(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 533 } 534 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 535 ierr = TSStep(ts);CHKERRQ(ierr); 536 if (!stack->solution_only && !tjsch->skip_trajectory) { 537 /* don't use the public interface as it will update the TSHistory: this need a better fix */ 538 ierr = TSTrajectorySet_Memory(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 539 } 540 ierr = TSEventHandler(ts);CHKERRQ(ierr); 541 if (!ts->steprollback) { 542 ierr = TSPostStep(ts);CHKERRQ(ierr); 543 } 544 } 545 ierr = TurnBackward(ts);CHKERRQ(ierr); 546 ts->trajectory->recomps += stepnumend-stepnumbegin; /* recomputation counter */ 547 ierr = TSSetStepNumber(ts,stepnumend);CHKERRQ(ierr); 548 tjsch->recompute = PETSC_FALSE; /* reset the flag for recompute mode */ 549 PetscFunctionReturn(0); 550 } 551 552 static PetscErrorCode TopLevelStore(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscInt localstepnum,PetscInt laststridesize,PetscBool *done) 553 { 554 Stack *stack = &tjsch->stack; 555 DiskStack *diskstack = &tjsch->diskstack; 556 PetscInt stridenum; 557 PetscErrorCode ierr; 558 559 PetscFunctionBegin; 560 *done = PETSC_FALSE; 561 stridenum = stepnum/tjsch->stride; 562 /* make sure saved checkpoint id starts from 1 563 skip last stride when using stridenum+1 564 skip first stride when using stridenum */ 565 if (stack->solution_only) { 566 if (tjsch->save_stack) { 567 if (localstepnum == tjsch->stride-1 && stepnum < tjsch->total_steps-laststridesize) { /* current step will be saved without going through stack */ 568 ierr = StackDumpAll(tj,ts,stack,stridenum+1);CHKERRQ(ierr); 569 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) diskstack->container[++diskstack->top] = stridenum+1; 570 *done = PETSC_TRUE; 571 } 572 } else { 573 if (localstepnum == 0 && stepnum < tjsch->total_steps-laststridesize) { 574 ierr = DumpSingle(tj,ts,stack,stridenum+1);CHKERRQ(ierr); 575 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) diskstack->container[++diskstack->top] = stridenum+1; 576 *done = PETSC_TRUE; 577 } 578 } 579 } else { 580 if (tjsch->save_stack) { 581 if (localstepnum == 0 && stepnum < tjsch->total_steps && stepnum != 0) { /* skip the first stride */ 582 ierr = StackDumpAll(tj,ts,stack,stridenum);CHKERRQ(ierr); 583 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) diskstack->container[++diskstack->top] = stridenum; 584 *done = PETSC_TRUE; 585 } 586 } else { 587 if (localstepnum == 1 && stepnum < tjsch->total_steps-laststridesize) { 588 ierr = DumpSingle(tj,ts,stack,stridenum+1);CHKERRQ(ierr); 589 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) diskstack->container[++diskstack->top] = stridenum+1; 590 *done = PETSC_TRUE; 591 } 592 } 593 } 594 PetscFunctionReturn(0); 595 } 596 597 static PetscErrorCode SetTrajN(TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 598 { 599 Stack *stack = &tjsch->stack; 600 StackElement e; 601 PetscErrorCode ierr; 602 603 PetscFunctionBegin; 604 /* skip the last step */ 605 if (ts->reason) { /* only affect the forward run */ 606 /* update total_steps in the end of forward run */ 607 if (stepnum != tjsch->total_steps) tjsch->total_steps = stepnum; 608 if (stack->solution_only) { 609 /* get rid of the solution at second last step */ 610 ierr = StackPop(stack,&e);CHKERRQ(ierr); 611 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 612 } 613 PetscFunctionReturn(0); 614 } 615 /* do not save trajectory at the recompute stage for solution_only mode */ 616 if (stack->solution_only && tjsch->recompute) PetscFunctionReturn(0); 617 /* skip the first step for no_solution_only mode */ 618 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 619 620 /* resize the stack */ 621 if (stack->top+1 == stack->stacksize) { 622 ierr = StackResize(stack,2*stack->stacksize);CHKERRQ(ierr); 623 } 624 /* update timenext for the previous step; necessary for step adaptivity */ 625 if (stack->top > -1) { 626 ierr = StackTop(stack,&e);CHKERRQ(ierr); 627 e->timenext = ts->ptime; 628 } 629 if (stepnum < stack->top) { 630 SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 631 } 632 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 633 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 634 ierr = StackPush(stack,e);CHKERRQ(ierr); 635 PetscFunctionReturn(0); 636 } 637 638 static PetscErrorCode SetTrajN_2(TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 639 { 640 Stack *stack = &tjsch->stack; 641 StackElement e; 642 PetscErrorCode ierr; 643 644 PetscFunctionBegin; 645 if (stack->top+1 == stack->stacksize) { 646 ierr = StackResize(stack,2*stack->stacksize);CHKERRQ(ierr); 647 } 648 /* update timenext for the previous step; necessary for step adaptivity */ 649 if (stack->top > -1) { 650 ierr = StackTop(stack,&e);CHKERRQ(ierr); 651 e->timenext = ts->ptime; 652 } 653 if (stepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 654 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 655 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 656 ierr = StackPush(stack,e);CHKERRQ(ierr); 657 PetscFunctionReturn(0); 658 } 659 660 static PetscErrorCode GetTrajN(TS ts,TJScheduler *tjsch,PetscInt stepnum) 661 { 662 Stack *stack = &tjsch->stack; 663 StackElement e; 664 PetscInt ns; 665 PetscErrorCode ierr; 666 667 PetscFunctionBegin; 668 if (stepnum == tjsch->total_steps) { 669 ierr = TurnBackward(ts);CHKERRQ(ierr); 670 PetscFunctionReturn(0); 671 } 672 /* restore a checkpoint */ 673 ierr = StackTop(stack,&e);CHKERRQ(ierr); 674 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 675 ierr = TSGetStages(ts,&ns,NULL);CHKERRQ(ierr); 676 if (stack->solution_only && ns) { /* recompute one step */ 677 ierr = TurnForwardWithStepsize(ts,e->timenext-e->time);CHKERRQ(ierr); 678 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 679 } 680 ierr = StackPop(stack,&e);CHKERRQ(ierr); 681 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 682 PetscFunctionReturn(0); 683 } 684 685 static PetscErrorCode GetTrajN_2(TS ts,TJScheduler *tjsch,PetscInt stepnum) 686 { 687 Stack *stack = &tjsch->stack; 688 StackElement e = NULL; 689 PetscErrorCode ierr; 690 691 PetscFunctionBegin; 692 ierr = StackFind(stack,&e,stepnum);CHKERRQ(ierr); 693 if (stepnum != e->stepnum) SETERRQ2(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Inconsistent steps! %D != %D",stepnum,e->stepnum); 694 ierr = UpdateTS(ts,stack,e,PETSC_FALSE);CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 698 static PetscErrorCode SetTrajTLNR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 699 { 700 Stack *stack = &tjsch->stack; 701 PetscInt localstepnum,laststridesize; 702 StackElement e; 703 PetscBool done; 704 PetscErrorCode ierr; 705 706 PetscFunctionBegin; 707 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 708 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 709 if (tjsch->save_stack && tjsch->recompute) PetscFunctionReturn(0); 710 711 localstepnum = stepnum%tjsch->stride; 712 /* (stridesize-1) checkpoints are saved in each stride; an extra point is added by StackDumpAll() */ 713 laststridesize = tjsch->total_steps%tjsch->stride; 714 if (!laststridesize) laststridesize = tjsch->stride; 715 716 if (!tjsch->recompute) { 717 ierr = TopLevelStore(tj,ts,tjsch,stepnum,localstepnum,laststridesize,&done);CHKERRQ(ierr); 718 if (!tjsch->save_stack && stepnum < tjsch->total_steps-laststridesize) PetscFunctionReturn(0); 719 } 720 if (!stack->solution_only && localstepnum == 0) PetscFunctionReturn(0); /* skip last point in each stride at recompute stage or last stride */ 721 if (stack->solution_only && localstepnum == tjsch->stride-1) PetscFunctionReturn(0); /* skip last step in each stride at recompute stage or last stride */ 722 723 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 724 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 725 ierr = StackPush(stack,e);CHKERRQ(ierr); 726 PetscFunctionReturn(0); 727 } 728 729 static PetscErrorCode GetTrajTLNR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 730 { 731 Stack *stack = &tjsch->stack; 732 PetscInt id,localstepnum,laststridesize; 733 StackElement e; 734 PetscErrorCode ierr; 735 736 PetscFunctionBegin; 737 if (stepnum == tjsch->total_steps) { 738 ierr = TurnBackward(ts);CHKERRQ(ierr); 739 PetscFunctionReturn(0); 740 } 741 742 localstepnum = stepnum%tjsch->stride; 743 laststridesize = tjsch->total_steps%tjsch->stride; 744 if (!laststridesize) laststridesize = tjsch->stride; 745 if (stack->solution_only) { 746 /* fill stack with info */ 747 if (localstepnum == 0 && tjsch->total_steps-stepnum >= laststridesize) { 748 id = stepnum/tjsch->stride; 749 if (tjsch->save_stack) { 750 ierr = StackLoadAll(tj,ts,stack,id);CHKERRQ(ierr); 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 ierr = TurnForward(ts);CHKERRQ(ierr); 758 ierr = ReCompute(ts,tjsch,(id-1)*tjsch->stride,id*tjsch->stride);CHKERRQ(ierr); 759 } 760 PetscFunctionReturn(0); 761 } 762 /* restore a checkpoint */ 763 ierr = StackPop(stack,&e);CHKERRQ(ierr); 764 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 765 tjsch->skip_trajectory = PETSC_TRUE; 766 ierr = TurnForward(ts);CHKERRQ(ierr); 767 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 768 tjsch->skip_trajectory = PETSC_FALSE; 769 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 770 } else { 771 /* fill stack with info */ 772 if (localstepnum == 0 && tjsch->total_steps-stepnum >= laststridesize) { 773 id = stepnum/tjsch->stride; 774 if (tjsch->save_stack) { 775 ierr = StackLoadAll(tj,ts,stack,id);CHKERRQ(ierr); 776 } else { 777 ierr = LoadSingle(tj,ts,stack,id);CHKERRQ(ierr); 778 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 779 ierr = ElementSet(ts,stack,&e,(id-1)*tjsch->stride+1,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 780 ierr = StackPush(stack,e);CHKERRQ(ierr); 781 ierr = TurnForward(ts);CHKERRQ(ierr); 782 ierr = ReCompute(ts,tjsch,e->stepnum,id*tjsch->stride);CHKERRQ(ierr); 783 } 784 PetscFunctionReturn(0); 785 } 786 /* restore a checkpoint */ 787 ierr = StackPop(stack,&e);CHKERRQ(ierr); 788 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 789 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 790 } 791 PetscFunctionReturn(0); 792 } 793 794 #if defined(PETSC_HAVE_REVOLVE) 795 static PetscErrorCode printwhattodo(PetscViewer viewer,PetscInt whattodo,RevolveCTX *rctx,PetscInt shift) 796 { 797 PetscErrorCode ierr; 798 799 PetscFunctionBegin; 800 if (!viewer) PetscFunctionReturn(0); 801 802 switch(whattodo) { 803 case 1: 804 ierr = PetscViewerASCIIPrintf(viewer,"Advance from %D to %D\n",rctx->oldcapo+shift,rctx->capo+shift);CHKERRQ(ierr); 805 break; 806 case 2: 807 ierr = PetscViewerASCIIPrintf(viewer,"Store in checkpoint number %D (located in RAM)\n",rctx->check);CHKERRQ(ierr); 808 break; 809 case 3: 810 ierr = PetscViewerASCIIPrintf(viewer,"First turn: Initialize adjoints and reverse first step\n");CHKERRQ(ierr); 811 break; 812 case 4: 813 ierr = PetscViewerASCIIPrintf(viewer,"Forward and reverse one step\n");CHKERRQ(ierr); 814 break; 815 case 5: 816 ierr = PetscViewerASCIIPrintf(viewer,"Restore in checkpoint number %D (located in RAM)\n",rctx->check);CHKERRQ(ierr); 817 break; 818 case 7: 819 ierr = PetscViewerASCIIPrintf(viewer,"Store in checkpoint number %D (located on disk)\n",rctx->check);CHKERRQ(ierr); 820 break; 821 case 8: 822 ierr = PetscViewerASCIIPrintf(viewer,"Restore in checkpoint number %D (located on disk)\n",rctx->check);CHKERRQ(ierr); 823 break; 824 case -1: 825 ierr = PetscViewerASCIIPrintf(viewer,"Error!");CHKERRQ(ierr); 826 break; 827 } 828 PetscFunctionReturn(0); 829 } 830 831 static PetscErrorCode printwhattodo2(PetscViewer viewer,PetscInt whattodo,RevolveCTX *rctx,PetscInt shift) 832 { 833 PetscErrorCode ierr; 834 835 PetscFunctionBegin; 836 if (!viewer) PetscFunctionReturn(0); 837 838 switch(whattodo) { 839 case 1: 840 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Advance from stride %D to stride %D\n",rctx->oldcapo+shift,rctx->capo+shift);CHKERRQ(ierr); 841 break; 842 case 2: 843 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Store in checkpoint number %D\n",rctx->check);CHKERRQ(ierr); 844 break; 845 case 3: 846 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] First turn: Initialize adjoints and reverse first stride\n");CHKERRQ(ierr); 847 break; 848 case 4: 849 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Forward and reverse one stride\n");CHKERRQ(ierr); 850 break; 851 case 5: 852 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Restore in checkpoint number %D\n",rctx->check);CHKERRQ(ierr); 853 break; 854 case 7: 855 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Store in top-level checkpoint number %D\n",rctx->check);CHKERRQ(ierr); 856 break; 857 case 8: 858 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Restore in top-level checkpoint number %D\n",rctx->check);CHKERRQ(ierr); 859 break; 860 case -1: 861 ierr = PetscViewerASCIIPrintf(viewer,"[Top Level] Error!");CHKERRQ(ierr); 862 break; 863 } 864 PetscFunctionReturn(0); 865 } 866 867 static PetscErrorCode InitRevolve(PetscInt fine,PetscInt snaps,RevolveCTX *rctx) 868 { 869 PetscFunctionBegin; 870 revolve_reset(); 871 revolve_create_offline(fine,snaps); 872 rctx->snaps_in = snaps; 873 rctx->fine = fine; 874 rctx->check = 0; 875 rctx->capo = 0; 876 rctx->reverseonestep = PETSC_FALSE; 877 /* check stepsleft? */ 878 PetscFunctionReturn(0); 879 } 880 881 static PetscErrorCode FastForwardRevolve(RevolveCTX *rctx) 882 { 883 PetscInt whattodo; 884 885 PetscFunctionBegin; 886 whattodo = 0; 887 while(whattodo!=3) { /* we have to fast forward revolve to the beginning of the backward sweep due to unfriendly revolve interface */ 888 whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 889 } 890 PetscFunctionReturn(0); 891 } 892 893 static PetscErrorCode ApplyRevolve(PetscViewer viewer,SchedulerType stype,RevolveCTX *rctx,PetscInt total_steps,PetscInt stepnum,PetscInt localstepnum,PetscBool toplevel,PetscInt *store) 894 { 895 PetscErrorCode ierr; 896 PetscInt shift,whattodo; 897 898 PetscFunctionBegin; 899 *store = 0; 900 if (rctx->stepsleft > 0) { /* advance the solution without checkpointing anything as Revolve requires */ 901 rctx->stepsleft--; 902 PetscFunctionReturn(0); 903 } 904 /* let Revolve determine what to do next */ 905 shift = stepnum-localstepnum; 906 rctx->oldcapo = rctx->capo; 907 rctx->capo = localstepnum; 908 909 if (!toplevel) whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 910 else whattodo = revolve2_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 911 if (stype == REVOLVE_ONLINE && whattodo == 8) whattodo = 5; 912 if (stype == REVOLVE_ONLINE && whattodo == 7) whattodo = 2; 913 if (!toplevel) {ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 914 else {ierr = printwhattodo2(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 915 if (whattodo == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in the Revolve library"); 916 if (whattodo == 1) { /* advance some time steps */ 917 if (stype == REVOLVE_ONLINE && rctx->capo >= total_steps-1) { 918 revolve_turn(total_steps,&rctx->capo,&rctx->fine); 919 if (!toplevel) {ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 920 else {ierr = printwhattodo2(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 921 } 922 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 923 } 924 if (whattodo == 3 || whattodo == 4) { /* ready for a reverse step */ 925 rctx->reverseonestep = PETSC_TRUE; 926 } 927 if (whattodo == 5) { /* restore a checkpoint and ask Revolve what to do next */ 928 rctx->oldcapo = rctx->capo; 929 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*/ 930 else whattodo = revolve2_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 931 if (!toplevel) {ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 932 else {ierr = printwhattodo2(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 933 if (whattodo == 3 || whattodo == 4) rctx->reverseonestep = PETSC_TRUE; 934 if (whattodo == 1) rctx->stepsleft = rctx->capo-rctx->oldcapo; 935 } 936 if (whattodo == 7) { /* save the checkpoint to disk */ 937 *store = 2; 938 rctx->oldcapo = rctx->capo; 939 whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); /* must return 1 */ 940 ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr); 941 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 942 } 943 if (whattodo == 2) { /* store a checkpoint to RAM and ask Revolve how many time steps to advance next */ 944 *store = 1; 945 rctx->oldcapo = rctx->capo; 946 if (!toplevel) whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); /* must return 1 */ 947 else whattodo = revolve2_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 948 if (!toplevel) {ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 949 else {ierr = printwhattodo2(viewer,whattodo,rctx,shift);CHKERRQ(ierr);} 950 if (stype == REVOLVE_ONLINE && rctx->capo >= total_steps-1) { 951 revolve_turn(total_steps,&rctx->capo,&rctx->fine); 952 ierr = printwhattodo(viewer,whattodo,rctx,shift);CHKERRQ(ierr); 953 } 954 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 955 } 956 PetscFunctionReturn(0); 957 } 958 959 static PetscErrorCode SetTrajROF(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 960 { 961 Stack *stack = &tjsch->stack; 962 PetscInt store; 963 StackElement e; 964 PetscErrorCode ierr; 965 966 PetscFunctionBegin; 967 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 968 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 969 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,stepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 970 if (store == 1) { 971 if (stepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 972 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 973 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 974 ierr = StackPush(stack,e);CHKERRQ(ierr); 975 } 976 PetscFunctionReturn(0); 977 } 978 979 static PetscErrorCode GetTrajROF(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 980 { 981 Stack *stack = &tjsch->stack; 982 PetscInt whattodo,shift,store; 983 StackElement e; 984 PetscErrorCode ierr; 985 986 PetscFunctionBegin; 987 if (stepnum == 0 || stepnum == tjsch->total_steps) { 988 ierr = TurnBackward(ts);CHKERRQ(ierr); 989 tjsch->rctx->reverseonestep = PETSC_FALSE; 990 PetscFunctionReturn(0); 991 } 992 /* restore a checkpoint */ 993 ierr = StackTop(stack,&e);CHKERRQ(ierr); 994 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 995 if (stack->solution_only) { /* start with restoring a checkpoint */ 996 tjsch->rctx->capo = stepnum; 997 tjsch->rctx->oldcapo = tjsch->rctx->capo; 998 shift = 0; 999 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); 1000 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1001 } else { /* 2 revolve actions: restore a checkpoint and then advance */ 1002 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,stepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1003 if (tj->monitor) { 1004 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1005 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); 1006 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1007 } 1008 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1009 } 1010 if (stack->solution_only || (!stack->solution_only && e->stepnum < stepnum)) { 1011 ierr = TurnForward(ts);CHKERRQ(ierr); 1012 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1013 } 1014 if ((stack->solution_only && e->stepnum+1 == stepnum) || (!stack->solution_only && e->stepnum == stepnum)) { 1015 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1016 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1017 } 1018 tjsch->rctx->reverseonestep = PETSC_FALSE; 1019 PetscFunctionReturn(0); 1020 } 1021 1022 static PetscErrorCode SetTrajRON(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 1023 { 1024 Stack *stack = &tjsch->stack; 1025 Vec *Y; 1026 PetscInt i,store; 1027 PetscReal timeprev; 1028 StackElement e; 1029 RevolveCTX *rctx = tjsch->rctx; 1030 PetscErrorCode ierr; 1031 1032 PetscFunctionBegin; 1033 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 1034 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 1035 ierr = ApplyRevolve(tj->monitor,tjsch->stype,rctx,tjsch->total_steps,stepnum,stepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1036 if (store == 1) { 1037 if (rctx->check != stack->top+1) { /* overwrite some non-top checkpoint in the stack */ 1038 ierr = StackFind(stack,&e,rctx->check);CHKERRQ(ierr); 1039 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 1040 if (stack->numY > 0 && !stack->solution_only) { 1041 ierr = TSGetStages(ts,&stack->numY,&Y);CHKERRQ(ierr); 1042 for (i=0;i<stack->numY;i++) { 1043 ierr = VecCopy(Y[i],e->Y[i]);CHKERRQ(ierr); 1044 } 1045 } 1046 e->stepnum = stepnum; 1047 e->time = time; 1048 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 1049 e->timeprev = timeprev; 1050 } else { 1051 if (stepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 1052 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1053 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 1054 ierr = StackPush(stack,e);CHKERRQ(ierr); 1055 } 1056 } 1057 PetscFunctionReturn(0); 1058 } 1059 1060 static PetscErrorCode GetTrajRON(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 1061 { 1062 Stack *stack = &tjsch->stack; 1063 PetscInt whattodo,shift; 1064 StackElement e; 1065 PetscErrorCode ierr; 1066 1067 PetscFunctionBegin; 1068 if (stepnum == 0 || stepnum == tjsch->total_steps) { 1069 ierr = TurnBackward(ts);CHKERRQ(ierr); 1070 tjsch->rctx->reverseonestep = PETSC_FALSE; 1071 PetscFunctionReturn(0); 1072 } 1073 tjsch->rctx->capo = stepnum; 1074 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1075 shift = 0; 1076 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); /* whattodo=restore */ 1077 if (whattodo == 8) whattodo = 5; 1078 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1079 /* restore a checkpoint */ 1080 ierr = StackFind(stack,&e,tjsch->rctx->check);CHKERRQ(ierr); 1081 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1082 if (!stack->solution_only) { /* whattodo must be 5 */ 1083 /* ask Revolve what to do next */ 1084 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1085 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*/ 1086 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1087 if (whattodo == 3 || whattodo == 4) tjsch->rctx->reverseonestep = PETSC_TRUE; 1088 if (whattodo == 1) tjsch->rctx->stepsleft = tjsch->rctx->capo-tjsch->rctx->oldcapo; 1089 if (tj->monitor) { 1090 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1091 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); 1092 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1093 } 1094 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1095 } 1096 if (stack->solution_only || (!stack->solution_only && e->stepnum < stepnum)) { 1097 ierr = TurnForward(ts);CHKERRQ(ierr); 1098 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1099 } 1100 tjsch->rctx->reverseonestep = PETSC_FALSE; 1101 PetscFunctionReturn(0); 1102 } 1103 1104 static PetscErrorCode SetTrajTLR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 1105 { 1106 Stack *stack = &tjsch->stack; 1107 PetscInt store,localstepnum,laststridesize; 1108 StackElement e; 1109 PetscBool done = PETSC_FALSE; 1110 PetscErrorCode ierr; 1111 1112 PetscFunctionBegin; 1113 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 1114 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 1115 1116 localstepnum = stepnum%tjsch->stride; 1117 laststridesize = tjsch->total_steps%tjsch->stride; 1118 if (!laststridesize) laststridesize = tjsch->stride; 1119 1120 if (!tjsch->recompute) { 1121 ierr = TopLevelStore(tj,ts,tjsch,stepnum,localstepnum,laststridesize,&done);CHKERRQ(ierr); 1122 /* revolve is needed for the last stride; different starting points for last stride between solutin_only and !solutin_only */ 1123 if (!stack->solution_only && !tjsch->save_stack && stepnum <= tjsch->total_steps-laststridesize) PetscFunctionReturn(0); 1124 if (stack->solution_only && !tjsch->save_stack && stepnum < tjsch->total_steps-laststridesize) PetscFunctionReturn(0); 1125 } 1126 if (tjsch->save_stack && done) { 1127 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1128 PetscFunctionReturn(0); 1129 } 1130 if (laststridesize < tjsch->stride) { 1131 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 */ 1132 ierr = InitRevolve(laststridesize,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1133 } 1134 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 */ 1135 ierr = InitRevolve(laststridesize,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1136 } 1137 } 1138 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,localstepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1139 if (store == 1) { 1140 if (localstepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 1141 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1142 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 1143 ierr = StackPush(stack,e);CHKERRQ(ierr); 1144 } 1145 PetscFunctionReturn(0); 1146 } 1147 1148 static PetscErrorCode GetTrajTLR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 1149 { 1150 Stack *stack = &tjsch->stack; 1151 PetscInt whattodo,shift; 1152 PetscInt localstepnum,stridenum,laststridesize,store; 1153 StackElement e; 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 localstepnum = stepnum%tjsch->stride; 1158 stridenum = stepnum/tjsch->stride; 1159 if (stepnum == tjsch->total_steps) { 1160 ierr = TurnBackward(ts);CHKERRQ(ierr); 1161 tjsch->rctx->reverseonestep = PETSC_FALSE; 1162 PetscFunctionReturn(0); 1163 } 1164 laststridesize = tjsch->total_steps%tjsch->stride; 1165 if (!laststridesize) laststridesize = tjsch->stride; 1166 if (stack->solution_only) { 1167 /* fill stack */ 1168 if (localstepnum == 0 && stepnum <= tjsch->total_steps-laststridesize) { 1169 if (tjsch->save_stack) { 1170 ierr = StackLoadAll(tj,ts,stack,stridenum);CHKERRQ(ierr); 1171 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1172 ierr = FastForwardRevolve(tjsch->rctx);CHKERRQ(ierr); 1173 tjsch->skip_trajectory = PETSC_TRUE; 1174 ierr = TurnForward(ts);CHKERRQ(ierr); 1175 ierr = ReCompute(ts,tjsch,stridenum*tjsch->stride-1,stridenum*tjsch->stride);CHKERRQ(ierr); 1176 tjsch->skip_trajectory = PETSC_FALSE; 1177 } else { 1178 ierr = LoadSingle(tj,ts,stack,stridenum);CHKERRQ(ierr); 1179 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1180 ierr = TurnForward(ts);CHKERRQ(ierr); 1181 ierr = ReCompute(ts,tjsch,(stridenum-1)*tjsch->stride,stridenum*tjsch->stride);CHKERRQ(ierr); 1182 } 1183 PetscFunctionReturn(0); 1184 } 1185 /* restore a checkpoint */ 1186 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1187 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1188 /* start with restoring a checkpoint */ 1189 tjsch->rctx->capo = stepnum; 1190 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1191 shift = stepnum-localstepnum; 1192 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); 1193 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1194 ierr = TurnForward(ts);CHKERRQ(ierr); 1195 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1196 if (e->stepnum+1 == stepnum) { 1197 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1198 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1199 } 1200 } else { 1201 /* fill stack with info */ 1202 if (localstepnum == 0 && tjsch->total_steps-stepnum >= laststridesize) { 1203 if (tjsch->save_stack) { 1204 ierr = StackLoadAll(tj,ts,stack,stridenum);CHKERRQ(ierr); 1205 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1206 ierr = FastForwardRevolve(tjsch->rctx);CHKERRQ(ierr); 1207 } else { 1208 ierr = LoadSingle(tj,ts,stack,stridenum);CHKERRQ(ierr); 1209 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1210 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,(stridenum-1)*tjsch->stride+1,1,PETSC_FALSE,&store);CHKERRQ(ierr); 1211 if (tj->monitor) { 1212 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1213 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); 1214 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1215 } 1216 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1217 ierr = ElementSet(ts,stack,&e,(stridenum-1)*tjsch->stride+1,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1218 ierr = StackPush(stack,e);CHKERRQ(ierr); 1219 ierr = TurnForward(ts);CHKERRQ(ierr); 1220 ierr = ReCompute(ts,tjsch,e->stepnum,stridenum*tjsch->stride);CHKERRQ(ierr); 1221 } 1222 PetscFunctionReturn(0); 1223 } 1224 /* restore a checkpoint */ 1225 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1226 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1227 /* 2 revolve actions: restore a checkpoint and then advance */ 1228 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,localstepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1229 if (tj->monitor) { 1230 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1231 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); 1232 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1233 } 1234 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1235 if (e->stepnum < stepnum) { 1236 ierr = TurnForward(ts);CHKERRQ(ierr); 1237 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1238 } 1239 if (e->stepnum == stepnum) { 1240 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1241 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1242 } 1243 } 1244 tjsch->rctx->reverseonestep = PETSC_FALSE; 1245 PetscFunctionReturn(0); 1246 } 1247 1248 static PetscErrorCode SetTrajTLTR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 1249 { 1250 Stack *stack = &tjsch->stack; 1251 PetscInt store,localstepnum,stridenum,laststridesize; 1252 StackElement e; 1253 PetscBool done = PETSC_FALSE; 1254 PetscErrorCode ierr; 1255 1256 PetscFunctionBegin; 1257 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 1258 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 1259 1260 localstepnum = stepnum%tjsch->stride; /* index at the bottom level (inside a stride) */ 1261 stridenum = stepnum/tjsch->stride; /* index at the top level */ 1262 laststridesize = tjsch->total_steps%tjsch->stride; 1263 if (!laststridesize) laststridesize = tjsch->stride; 1264 if (stack->solution_only && localstepnum == 0 && !tjsch->rctx2->reverseonestep) { 1265 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); 1266 if (laststridesize < tjsch->stride && stepnum == tjsch->total_steps-laststridesize) { 1267 ierr = InitRevolve(laststridesize,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1268 } 1269 } 1270 if (!stack->solution_only && localstepnum == 1 && !tjsch->rctx2->reverseonestep) { 1271 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); 1272 if (laststridesize < tjsch->stride && stepnum == tjsch->total_steps-laststridesize+1) { 1273 ierr = InitRevolve(laststridesize,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1274 } 1275 } 1276 if (tjsch->store_stride) { 1277 ierr = TopLevelStore(tj,ts,tjsch,stepnum,localstepnum,laststridesize,&done);CHKERRQ(ierr); 1278 if (done) { 1279 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1280 PetscFunctionReturn(0); 1281 } 1282 } 1283 if (stepnum < tjsch->total_steps-laststridesize) { 1284 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 */ 1285 if (!tjsch->save_stack && !tjsch->rctx2->reverseonestep) PetscFunctionReturn(0); /* store operation does not require revolve be called at bottom level */ 1286 } 1287 /* 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 */ 1288 if (!stack->solution_only && localstepnum == 0 && stepnum != tjsch->total_steps && !tjsch->recompute) PetscFunctionReturn(0); 1289 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,localstepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1290 if (store == 1) { 1291 if (localstepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 1292 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1293 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 1294 ierr = StackPush(stack,e);CHKERRQ(ierr); 1295 } 1296 PetscFunctionReturn(0); 1297 } 1298 1299 static PetscErrorCode GetTrajTLTR(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 1300 { 1301 Stack *stack = &tjsch->stack; 1302 DiskStack *diskstack = &tjsch->diskstack; 1303 PetscInt whattodo,shift; 1304 PetscInt localstepnum,stridenum,restoredstridenum,laststridesize,store; 1305 StackElement e; 1306 PetscErrorCode ierr; 1307 1308 PetscFunctionBegin; 1309 localstepnum = stepnum%tjsch->stride; 1310 stridenum = stepnum/tjsch->stride; 1311 if (stepnum == tjsch->total_steps) { 1312 ierr = TurnBackward(ts);CHKERRQ(ierr); 1313 tjsch->rctx->reverseonestep = PETSC_FALSE; 1314 PetscFunctionReturn(0); 1315 } 1316 laststridesize = tjsch->total_steps%tjsch->stride; 1317 if (!laststridesize) laststridesize = tjsch->stride; 1318 /* 1319 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: 1320 Case 1 (save_stack) 1321 Restore a disk checkpoint; update TS with the last element in the restored data; recompute to the current point. 1322 Case 2 (!save_stack) 1323 Restore a disk checkpoint; update TS with the restored point; recompute to the current point. 1324 */ 1325 if (localstepnum == 0 && stepnum <= tjsch->total_steps-laststridesize) { 1326 /* restore the top element in the stack for disk checkpoints */ 1327 restoredstridenum = diskstack->container[diskstack->top]; 1328 tjsch->rctx2->reverseonestep = PETSC_FALSE; 1329 /* top-level revolve must be applied before current step, just like the solution_only mode for single-level revolve */ 1330 if (!tjsch->save_stack && stack->solution_only) { /* start with restoring a checkpoint */ 1331 tjsch->rctx2->capo = stridenum; 1332 tjsch->rctx2->oldcapo = tjsch->rctx2->capo; 1333 shift = 0; 1334 whattodo = revolve2_action(&tjsch->rctx2->check,&tjsch->rctx2->capo,&tjsch->rctx2->fine,tjsch->rctx2->snaps_in,&tjsch->rctx2->info,&tjsch->rctx2->where); 1335 ierr = printwhattodo2(tj->monitor,whattodo,tjsch->rctx2,shift);CHKERRQ(ierr); 1336 } else { /* 2 revolve actions: restore a checkpoint and then advance */ 1337 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); 1338 if (tj->monitor) { 1339 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1340 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); 1341 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1342 } 1343 if (!tjsch->rctx2->reverseonestep && tjsch->rctx2->stepsleft > 0) tjsch->rctx2->stepsleft--; 1344 } 1345 /* fill stack */ 1346 if (stack->solution_only) { 1347 if (tjsch->save_stack) { 1348 if (restoredstridenum < stridenum) { 1349 ierr = StackLoadLast(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1350 } else { 1351 ierr = StackLoadAll(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1352 } 1353 /* recompute one step ahead */ 1354 tjsch->skip_trajectory = PETSC_TRUE; 1355 ierr = TurnForward(ts);CHKERRQ(ierr); 1356 ierr = ReCompute(ts,tjsch,stridenum*tjsch->stride-1,stridenum*tjsch->stride);CHKERRQ(ierr); 1357 tjsch->skip_trajectory = PETSC_FALSE; 1358 if (restoredstridenum < stridenum) { 1359 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1360 ierr = TurnForward(ts);CHKERRQ(ierr); 1361 ierr = ReCompute(ts,tjsch,restoredstridenum*tjsch->stride,stepnum);CHKERRQ(ierr); 1362 } else { /* stack ready, fast forward revolve status */ 1363 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1364 ierr = FastForwardRevolve(tjsch->rctx);CHKERRQ(ierr); 1365 } 1366 } else { 1367 ierr = LoadSingle(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1368 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1369 ierr = TurnForward(ts);CHKERRQ(ierr); 1370 ierr = ReCompute(ts,tjsch,(restoredstridenum-1)*tjsch->stride,stepnum);CHKERRQ(ierr); 1371 } 1372 } else { 1373 if (tjsch->save_stack) { 1374 if (restoredstridenum < stridenum) { 1375 ierr = StackLoadLast(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1376 /* reset revolve */ 1377 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1378 ierr = TurnForward(ts);CHKERRQ(ierr); 1379 ierr = ReCompute(ts,tjsch,restoredstridenum*tjsch->stride,stepnum);CHKERRQ(ierr); 1380 } else { /* stack ready, fast forward revolve status */ 1381 ierr = StackLoadAll(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1382 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1383 ierr = FastForwardRevolve(tjsch->rctx);CHKERRQ(ierr); 1384 } 1385 } else { 1386 ierr = LoadSingle(tj,ts,stack,restoredstridenum);CHKERRQ(ierr); 1387 ierr = InitRevolve(tjsch->stride,tjsch->max_cps_ram,tjsch->rctx);CHKERRQ(ierr); 1388 /* push first element to stack */ 1389 if (tjsch->store_stride || tjsch->rctx2->reverseonestep) { 1390 shift = (restoredstridenum-1)*tjsch->stride-localstepnum; 1391 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,(restoredstridenum-1)*tjsch->stride+1,1,PETSC_FALSE,&store);CHKERRQ(ierr); 1392 if (tj->monitor) { 1393 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1394 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); 1395 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1396 } 1397 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1398 ierr = ElementSet(ts,stack,&e,(restoredstridenum-1)*tjsch->stride+1,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 1399 ierr = StackPush(stack,e);CHKERRQ(ierr); 1400 } 1401 ierr = TurnForward(ts);CHKERRQ(ierr); 1402 ierr = ReCompute(ts,tjsch,(restoredstridenum-1)*tjsch->stride+1,stepnum);CHKERRQ(ierr); 1403 } 1404 } 1405 if (restoredstridenum == stridenum) diskstack->top--; 1406 tjsch->rctx->reverseonestep = PETSC_FALSE; 1407 PetscFunctionReturn(0); 1408 } 1409 1410 if (stack->solution_only) { 1411 /* restore a checkpoint */ 1412 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1413 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1414 /* start with restoring a checkpoint */ 1415 tjsch->rctx->capo = stepnum; 1416 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1417 shift = stepnum-localstepnum; 1418 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); 1419 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1420 ierr = TurnForward(ts);CHKERRQ(ierr); 1421 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1422 if (e->stepnum+1 == stepnum) { 1423 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1424 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1425 } 1426 } else { 1427 /* restore a checkpoint */ 1428 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1429 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1430 /* 2 revolve actions: restore a checkpoint and then advance */ 1431 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,localstepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1432 if (tj->monitor) { 1433 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1434 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); 1435 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1436 } 1437 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1438 if (e->stepnum < stepnum) { 1439 ierr = TurnForward(ts);CHKERRQ(ierr); 1440 ierr = ReCompute(ts,tjsch,e->stepnum,stepnum);CHKERRQ(ierr); 1441 } 1442 if (e->stepnum == stepnum) { 1443 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1444 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1445 } 1446 } 1447 tjsch->rctx->reverseonestep = PETSC_FALSE; 1448 PetscFunctionReturn(0); 1449 } 1450 1451 static PetscErrorCode SetTrajRMS(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum,PetscReal time,Vec X) 1452 { 1453 Stack *stack = &tjsch->stack; 1454 PetscInt store; 1455 StackElement e; 1456 PetscErrorCode ierr; 1457 1458 PetscFunctionBegin; 1459 if (!stack->solution_only && stepnum == 0) PetscFunctionReturn(0); 1460 if (stack->solution_only && stepnum == tjsch->total_steps) PetscFunctionReturn(0); 1461 ierr = ApplyRevolve(tj->monitor,tjsch->stype,tjsch->rctx,tjsch->total_steps,stepnum,stepnum,PETSC_FALSE,&store);CHKERRQ(ierr); 1462 if (store == 1){ 1463 if (stepnum < stack->top) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 1464 ierr = ElementCreate(ts,stack,&e);CHKERRQ(ierr); 1465 ierr = ElementSet(ts,stack,&e,stepnum,time,X);CHKERRQ(ierr); 1466 ierr = StackPush(stack,e);CHKERRQ(ierr); 1467 } else if (store == 2) { 1468 ierr = DumpSingle(tj,ts,stack,tjsch->rctx->check+1);CHKERRQ(ierr); 1469 } 1470 PetscFunctionReturn(0); 1471 } 1472 1473 static PetscErrorCode GetTrajRMS(TSTrajectory tj,TS ts,TJScheduler *tjsch,PetscInt stepnum) 1474 { 1475 Stack *stack = &tjsch->stack; 1476 PetscInt whattodo,shift; 1477 PetscInt restart; 1478 PetscBool ondisk; 1479 StackElement e; 1480 PetscErrorCode ierr; 1481 1482 PetscFunctionBegin; 1483 if (stepnum == 0 || stepnum == tjsch->total_steps) { 1484 ierr = TurnBackward(ts);CHKERRQ(ierr); 1485 tjsch->rctx->reverseonestep = PETSC_FALSE; 1486 PetscFunctionReturn(0); 1487 } 1488 tjsch->rctx->capo = stepnum; 1489 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1490 shift = 0; 1491 whattodo = revolve_action(&tjsch->rctx->check,&tjsch->rctx->capo,&tjsch->rctx->fine,tjsch->rctx->snaps_in,&tjsch->rctx->info,&tjsch->rctx->where); /* whattodo=restore */ 1492 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1493 /* restore a checkpoint */ 1494 restart = tjsch->rctx->capo; 1495 if (!tjsch->rctx->where) { 1496 ondisk = PETSC_TRUE; 1497 ierr = LoadSingle(tj,ts,stack,tjsch->rctx->check+1);CHKERRQ(ierr); 1498 ierr = TurnBackward(ts);CHKERRQ(ierr); 1499 } else { 1500 ondisk = PETSC_FALSE; 1501 ierr = StackTop(stack,&e);CHKERRQ(ierr); 1502 ierr = UpdateTS(ts,stack,e,PETSC_TRUE);CHKERRQ(ierr); 1503 } 1504 if (!stack->solution_only) { /* whattodo must be 5 or 8 */ 1505 /* ask Revolve what to do next */ 1506 tjsch->rctx->oldcapo = tjsch->rctx->capo; 1507 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*/ 1508 ierr = printwhattodo(tj->monitor,whattodo,tjsch->rctx,shift);CHKERRQ(ierr); 1509 if (whattodo == 3 || whattodo == 4) tjsch->rctx->reverseonestep = PETSC_TRUE; 1510 if (whattodo == 1) tjsch->rctx->stepsleft = tjsch->rctx->capo-tjsch->rctx->oldcapo; 1511 if (tj->monitor) { 1512 ierr = PetscViewerASCIIAddTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1513 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); 1514 ierr = PetscViewerASCIISubtractTab(tj->monitor,((PetscObject)tj)->tablevel);CHKERRQ(ierr); 1515 } 1516 if (!tjsch->rctx->reverseonestep && tjsch->rctx->stepsleft > 0) tjsch->rctx->stepsleft--; 1517 restart++; /* skip one step */ 1518 } 1519 if (stack->solution_only || (!stack->solution_only && restart < stepnum)) { 1520 ierr = TurnForward(ts);CHKERRQ(ierr); 1521 ierr = ReCompute(ts,tjsch,restart,stepnum);CHKERRQ(ierr); 1522 } 1523 if (!ondisk && ( (stack->solution_only && e->stepnum+1 == stepnum) || (!stack->solution_only && e->stepnum == stepnum) )) { 1524 ierr = StackPop(stack,&e);CHKERRQ(ierr); 1525 ierr = ElementDestroy(stack,e);CHKERRQ(ierr); 1526 } 1527 tjsch->rctx->reverseonestep = PETSC_FALSE; 1528 PetscFunctionReturn(0); 1529 } 1530 #endif 1531 1532 static PetscErrorCode TSTrajectorySet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X) 1533 { 1534 TJScheduler *tjsch = (TJScheduler*)tj->data; 1535 PetscErrorCode ierr; 1536 1537 PetscFunctionBegin; 1538 if (!tjsch->recompute) { /* use global stepnum in the forward sweep */ 1539 ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr); 1540 } 1541 /* for consistency */ 1542 if (!tjsch->recompute && stepnum == 0) ts->ptime_prev = ts->ptime-ts->time_step; 1543 switch (tjsch->stype) { 1544 case NONE: 1545 if (tj->adjoint_solve_mode) { 1546 ierr = SetTrajN(ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1547 } else { 1548 ierr = SetTrajN_2(ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1549 } 1550 break; 1551 case TWO_LEVEL_NOREVOLVE: 1552 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1553 ierr = SetTrajTLNR(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1554 break; 1555 #if defined(PETSC_HAVE_REVOLVE) 1556 case TWO_LEVEL_REVOLVE: 1557 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1558 ierr = SetTrajTLR(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1559 break; 1560 case TWO_LEVEL_TWO_REVOLVE: 1561 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1562 ierr = SetTrajTLTR(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1563 break; 1564 case REVOLVE_OFFLINE: 1565 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1566 ierr = SetTrajROF(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1567 break; 1568 case REVOLVE_ONLINE: 1569 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1570 ierr = SetTrajRON(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1571 break; 1572 case REVOLVE_MULTISTAGE: 1573 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1574 ierr = SetTrajRMS(tj,ts,tjsch,stepnum,time,X);CHKERRQ(ierr); 1575 break; 1576 #endif 1577 default: 1578 break; 1579 } 1580 PetscFunctionReturn(0); 1581 } 1582 1583 static PetscErrorCode TSTrajectoryGet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t) 1584 { 1585 TJScheduler *tjsch = (TJScheduler*)tj->data; 1586 PetscErrorCode ierr; 1587 1588 PetscFunctionBegin; 1589 if (tj->adjoint_solve_mode && stepnum == 0) PetscFunctionReturn(0); 1590 switch (tjsch->stype) { 1591 case NONE: 1592 if (tj->adjoint_solve_mode) { 1593 ierr = GetTrajN(ts,tjsch,stepnum);CHKERRQ(ierr); 1594 } else { 1595 ierr = GetTrajN_2(ts,tjsch,stepnum);CHKERRQ(ierr); 1596 } 1597 break; 1598 case TWO_LEVEL_NOREVOLVE: 1599 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1600 ierr = GetTrajTLNR(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1601 break; 1602 #if defined(PETSC_HAVE_REVOLVE) 1603 case TWO_LEVEL_REVOLVE: 1604 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1605 ierr = GetTrajTLR(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1606 break; 1607 case TWO_LEVEL_TWO_REVOLVE: 1608 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1609 ierr = GetTrajTLTR(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1610 break; 1611 case REVOLVE_OFFLINE: 1612 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1613 ierr = GetTrajROF(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1614 break; 1615 case REVOLVE_ONLINE: 1616 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1617 ierr = GetTrajRON(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1618 break; 1619 case REVOLVE_MULTISTAGE: 1620 if (!tj->adjoint_solve_mode) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_SUP,"Not implemented"); 1621 ierr = GetTrajRMS(tj,ts,tjsch,stepnum);CHKERRQ(ierr); 1622 break; 1623 #endif 1624 default: 1625 break; 1626 } 1627 PetscFunctionReturn(0); 1628 } 1629 1630 PETSC_UNUSED static PetscErrorCode TSTrajectorySetStride_Memory(TSTrajectory tj,PetscInt stride) 1631 { 1632 TJScheduler *tjsch = (TJScheduler*)tj->data; 1633 1634 PetscFunctionBegin; 1635 tjsch->stride = stride; 1636 PetscFunctionReturn(0); 1637 } 1638 1639 PETSC_UNUSED static PetscErrorCode TSTrajectorySetMaxCpsRAM_Memory(TSTrajectory tj,PetscInt max_cps_ram) 1640 { 1641 TJScheduler *tjsch = (TJScheduler*)tj->data; 1642 1643 PetscFunctionBegin; 1644 tjsch->max_cps_ram = max_cps_ram; 1645 PetscFunctionReturn(0); 1646 } 1647 1648 PETSC_UNUSED static PetscErrorCode TSTrajectorySetMaxCpsDisk_Memory(TSTrajectory tj,PetscInt max_cps_disk) 1649 { 1650 TJScheduler *tjsch = (TJScheduler*)tj->data; 1651 1652 PetscFunctionBegin; 1653 tjsch->max_cps_disk = max_cps_disk; 1654 PetscFunctionReturn(0); 1655 } 1656 1657 #if defined(PETSC_HAVE_REVOLVE) 1658 PETSC_UNUSED static PetscErrorCode TSTrajectorySetRevolveOnline(TSTrajectory tj,PetscBool use_online) 1659 { 1660 TJScheduler *tjsch = (TJScheduler*)tj->data; 1661 1662 PetscFunctionBegin; 1663 tjsch->use_online = use_online; 1664 PetscFunctionReturn(0); 1665 } 1666 #endif 1667 1668 PETSC_UNUSED static PetscErrorCode TSTrajectorySetSaveStack(TSTrajectory tj,PetscBool save_stack) 1669 { 1670 TJScheduler *tjsch = (TJScheduler*)tj->data; 1671 1672 PetscFunctionBegin; 1673 tjsch->save_stack = save_stack; 1674 PetscFunctionReturn(0); 1675 } 1676 1677 PETSC_UNUSED static PetscErrorCode TSTrajectorySetUseDRAM(TSTrajectory tj,PetscBool use_dram) 1678 { 1679 TJScheduler *tjsch = (TJScheduler*)tj->data; 1680 1681 PetscFunctionBegin; 1682 tjsch->stack.use_dram = use_dram; 1683 PetscFunctionReturn(0); 1684 } 1685 1686 static PetscErrorCode TSTrajectorySetFromOptions_Memory(PetscOptionItems *PetscOptionsObject,TSTrajectory tj) 1687 { 1688 TJScheduler *tjsch = (TJScheduler*)tj->data; 1689 PetscErrorCode ierr; 1690 1691 PetscFunctionBegin; 1692 ierr = PetscOptionsHead(PetscOptionsObject,"Memory based TS trajectory options");CHKERRQ(ierr); 1693 { 1694 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); 1695 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); 1696 ierr = PetscOptionsInt("-ts_trajectory_stride","Stride to save checkpoints to file","TSTrajectorySetStride_Memory",tjsch->stride,&tjsch->stride,NULL);CHKERRQ(ierr); 1697 #if defined(PETSC_HAVE_REVOLVE) 1698 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); 1699 #endif 1700 ierr = PetscOptionsBool("-ts_trajectory_save_stack","Save all stack to disk","TSTrajectorySetSaveStack",tjsch->save_stack,&tjsch->save_stack,NULL);CHKERRQ(ierr); 1701 ierr = PetscOptionsBool("-ts_trajectory_use_dram","Use DRAM for checkpointing","TSTrajectorySetUseDRAM",tjsch->stack.use_dram,&tjsch->stack.use_dram,NULL);CHKERRQ(ierr); 1702 } 1703 ierr = PetscOptionsTail();CHKERRQ(ierr); 1704 tjsch->stack.solution_only = tj->solution_only; 1705 PetscFunctionReturn(0); 1706 } 1707 1708 static PetscErrorCode TSTrajectorySetUp_Memory(TSTrajectory tj,TS ts) 1709 { 1710 TJScheduler *tjsch = (TJScheduler*)tj->data; 1711 Stack *stack = &tjsch->stack; 1712 #if defined(PETSC_HAVE_REVOLVE) 1713 RevolveCTX *rctx,*rctx2; 1714 DiskStack *diskstack = &tjsch->diskstack; 1715 PetscInt diskblocks; 1716 #endif 1717 PetscInt numY,total_steps; 1718 PetscBool fixedtimestep; 1719 PetscErrorCode ierr; 1720 1721 PetscFunctionBegin; 1722 if (ts->adapt) { 1723 ierr = PetscObjectTypeCompare((PetscObject)ts->adapt,TSADAPTNONE,&fixedtimestep);CHKERRQ(ierr); 1724 } else fixedtimestep = PETSC_TRUE; 1725 total_steps = (PetscInt)(PetscCeilReal((ts->max_time-ts->ptime)/ts->time_step)); 1726 total_steps = total_steps < 0 ? PETSC_MAX_INT : total_steps; 1727 if (fixedtimestep) tjsch->total_steps = PetscMin(ts->max_steps,total_steps); 1728 if (tjsch->max_cps_ram > 0) stack->stacksize = tjsch->max_cps_ram; 1729 1730 if (tjsch->stride > 1) { /* two level mode */ 1731 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."); 1732 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 */ 1733 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 */ 1734 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 */ 1735 } else { /* single level mode */ 1736 if (fixedtimestep) { 1737 if (tjsch->max_cps_ram >= tjsch->total_steps-1 || tjsch->max_cps_ram < 1) tjsch->stype = NONE; /* checkpoint all */ 1738 else tjsch->stype = (tjsch->max_cps_disk>1) ? REVOLVE_MULTISTAGE : REVOLVE_OFFLINE; 1739 } else tjsch->stype = NONE; /* checkpoint all for adaptive time step */ 1740 #if defined(PETSC_HAVE_REVOLVE) 1741 if (tjsch->use_online) tjsch->stype = REVOLVE_ONLINE; /* trick into online (for testing purpose only) */ 1742 #endif 1743 } 1744 1745 if (tjsch->stype > TWO_LEVEL_NOREVOLVE) { 1746 #ifndef PETSC_HAVE_REVOLVE 1747 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."); 1748 #else 1749 switch (tjsch->stype) { 1750 case TWO_LEVEL_REVOLVE: 1751 revolve_create_offline(tjsch->stride,tjsch->max_cps_ram); 1752 break; 1753 case TWO_LEVEL_TWO_REVOLVE: 1754 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. */ 1755 diskstack->stacksize = diskblocks; 1756 revolve_create_offline(tjsch->stride,tjsch->max_cps_ram); 1757 revolve2_create_offline((tjsch->total_steps+tjsch->stride-1)/tjsch->stride,diskblocks); 1758 ierr = PetscCalloc1(1,&rctx2);CHKERRQ(ierr); 1759 rctx2->snaps_in = diskblocks; 1760 rctx2->reverseonestep = PETSC_FALSE; 1761 rctx2->check = 0; 1762 rctx2->oldcapo = 0; 1763 rctx2->capo = 0; 1764 rctx2->info = 2; 1765 rctx2->fine = (tjsch->total_steps+tjsch->stride-1)/tjsch->stride; 1766 tjsch->rctx2 = rctx2; 1767 diskstack->top = -1; 1768 ierr = PetscMalloc1(diskstack->stacksize*sizeof(PetscInt),&diskstack->container);CHKERRQ(ierr); 1769 break; 1770 case REVOLVE_OFFLINE: 1771 revolve_create_offline(tjsch->total_steps,tjsch->max_cps_ram); 1772 break; 1773 case REVOLVE_ONLINE: 1774 stack->stacksize = tjsch->max_cps_ram; 1775 revolve_create_online(tjsch->max_cps_ram); 1776 break; 1777 case REVOLVE_MULTISTAGE: 1778 revolve_create_multistage(tjsch->total_steps,tjsch->max_cps_ram+tjsch->max_cps_disk,tjsch->max_cps_ram); 1779 break; 1780 default: 1781 break; 1782 } 1783 ierr = PetscCalloc1(1,&rctx);CHKERRQ(ierr); 1784 rctx->snaps_in = tjsch->max_cps_ram; /* for theta methods snaps_in=2*max_cps_ram */ 1785 rctx->reverseonestep = PETSC_FALSE; 1786 rctx->check = 0; 1787 rctx->oldcapo = 0; 1788 rctx->capo = 0; 1789 rctx->info = 2; 1790 rctx->fine = (tjsch->stride > 1) ? tjsch->stride : tjsch->total_steps; 1791 tjsch->rctx = rctx; 1792 if (tjsch->stype == REVOLVE_ONLINE) rctx->fine = -1; 1793 #endif 1794 } else { 1795 if (tjsch->stype == TWO_LEVEL_NOREVOLVE) stack->stacksize = tjsch->stride-1; /* need tjsch->stride-1 at most */ 1796 if (tjsch->stype == NONE) { 1797 if (fixedtimestep) stack->stacksize = stack->solution_only ? tjsch->total_steps : tjsch->total_steps-1; 1798 else { /* adaptive time step */ 1799 /* if max_cps_ram is not specified, use maximal allowed number of steps for stack size */ 1800 if(tjsch->max_cps_ram == -1) stack->stacksize = ts->max_steps < PETSC_MAX_INT ? ts->max_steps : 10000; 1801 tjsch->total_steps = stack->solution_only ? stack->stacksize : stack->stacksize+1; /* will be updated as time integration advances */ 1802 } 1803 } 1804 } 1805 1806 tjsch->recompute = PETSC_FALSE; 1807 ierr = TSGetStages(ts,&numY,PETSC_IGNORE);CHKERRQ(ierr); 1808 ierr = StackCreate(stack,stack->stacksize,numY);CHKERRQ(ierr); 1809 PetscFunctionReturn(0); 1810 } 1811 1812 static PetscErrorCode TSTrajectoryReset_Memory(TSTrajectory tj) 1813 { 1814 TJScheduler *tjsch = (TJScheduler*)tj->data; 1815 PetscErrorCode ierr; 1816 1817 PetscFunctionBegin; 1818 if (tjsch->stype > TWO_LEVEL_NOREVOLVE) { 1819 #if defined(PETSC_HAVE_REVOLVE) 1820 revolve_reset(); 1821 if (tjsch->stype == TWO_LEVEL_TWO_REVOLVE) { 1822 revolve2_reset(); 1823 ierr = PetscFree(tjsch->diskstack.container);CHKERRQ(ierr); 1824 } 1825 #endif 1826 } 1827 ierr = StackDestroy(&tjsch->stack);CHKERRQ(ierr); 1828 #if defined(PETSC_HAVE_REVOLVE) 1829 if (tjsch->stype > TWO_LEVEL_NOREVOLVE) { 1830 ierr = PetscFree(tjsch->rctx);CHKERRQ(ierr); 1831 ierr = PetscFree(tjsch->rctx2);CHKERRQ(ierr); 1832 } 1833 #endif 1834 PetscFunctionReturn(0); 1835 } 1836 1837 static PetscErrorCode TSTrajectoryDestroy_Memory(TSTrajectory tj) 1838 { 1839 TJScheduler *tjsch = (TJScheduler*)tj->data; 1840 PetscErrorCode ierr; 1841 1842 PetscFunctionBegin; 1843 ierr = PetscFree(tjsch);CHKERRQ(ierr); 1844 PetscFunctionReturn(0); 1845 } 1846 1847 /*MC 1848 TSTRAJECTORYMEMORY - Stores each solution of the ODE/ADE in memory 1849 1850 Level: intermediate 1851 1852 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType() 1853 1854 M*/ 1855 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory tj,TS ts) 1856 { 1857 TJScheduler *tjsch; 1858 PetscErrorCode ierr; 1859 1860 PetscFunctionBegin; 1861 tj->ops->set = TSTrajectorySet_Memory; 1862 tj->ops->get = TSTrajectoryGet_Memory; 1863 tj->ops->setup = TSTrajectorySetUp_Memory; 1864 tj->ops->setfromoptions = TSTrajectorySetFromOptions_Memory; 1865 tj->ops->reset = TSTrajectoryReset_Memory; 1866 tj->ops->destroy = TSTrajectoryDestroy_Memory; 1867 1868 ierr = PetscCalloc1(1,&tjsch);CHKERRQ(ierr); 1869 tjsch->stype = NONE; 1870 tjsch->max_cps_ram = -1; /* -1 indicates that it is not set */ 1871 tjsch->max_cps_disk = -1; /* -1 indicates that it is not set */ 1872 tjsch->stride = 0; /* if not zero, two-level checkpointing will be used */ 1873 #if defined(PETSC_HAVE_REVOLVE) 1874 tjsch->use_online = PETSC_FALSE; 1875 #endif 1876 tjsch->save_stack = PETSC_TRUE; 1877 1878 tjsch->stack.solution_only = tj->solution_only; 1879 1880 tj->data = tjsch; 1881 PetscFunctionReturn(0); 1882 } 1883