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