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