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