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