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