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