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