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