1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 2 #include <petscsys.h> 3 4 #ifdef PETSC_HAVE_REVOLVE 5 #include <revolve_c.h> 6 #endif 7 8 PetscLogEvent Disk_Write, Disk_Read; 9 10 typedef enum {NONE,TWO_LEVEL_NOREVOLVE,TWO_LEVEL_REVOLVE,REVOLVE_OFFLINE,REVOLVE_ONLINE,REVOLVE_MULTISTAGE} SchedulerType; 11 12 typedef struct _StackElement { 13 PetscInt stepnum; 14 Vec X; 15 Vec *Y; 16 PetscReal time; 17 PetscReal timeprev; 18 } *StackElement; 19 20 typedef struct _RevolveCTX { 21 PetscBool reverseonestep; 22 PetscInt where; 23 PetscInt snaps_in; 24 PetscInt stepsleft; 25 PetscInt check; 26 PetscInt oldcapo; 27 PetscInt capo; 28 PetscInt fine; 29 PetscInt info; 30 } RevolveCTX; 31 32 typedef struct _Stack { 33 SchedulerType stype; 34 PetscBool use_online; 35 PetscBool recompute; 36 PetscBool solution_only; 37 PetscBool save_stack; 38 MPI_Comm comm; 39 RevolveCTX *rctx; 40 PetscInt max_cps_ram; /* maximum checkpoints in RAM */ 41 PetscInt max_cps_disk; /* maximum checkpoints on disk */ 42 PetscInt stride; 43 PetscInt total_steps; /* total number of steps */ 44 PetscInt numY; 45 PetscInt stacksize; 46 PetscInt top; /* top of the stack */ 47 StackElement *stack; /* container */ 48 } Stack; 49 50 static PetscErrorCode StackCreate(MPI_Comm,Stack *,PetscInt,PetscInt); 51 static PetscErrorCode StackDestroy(Stack*); 52 static PetscErrorCode StackPush(Stack*,StackElement); 53 static PetscErrorCode StackPop(Stack*,StackElement*); 54 static PetscErrorCode StackTop(Stack*,StackElement*); 55 static PetscErrorCode StackDumpAll(TS,Stack*,PetscInt); 56 static PetscErrorCode StackLoadAll(TS,Stack*,PetscInt); 57 58 #ifdef PETSC_HAVE_REVOLVE 59 static void printwhattodo(PetscInt whattodo,RevolveCTX *rctx,PetscInt shift) 60 { 61 switch(whattodo) { 62 case 1: 63 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mAdvance from %D to %D\033[0m\n",rctx->oldcapo+shift,rctx->capo+shift); 64 break; 65 case 2: 66 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mStore in checkpoint number %D (located in RAM)\033[0m\n",rctx->check); 67 break; 68 case 3: 69 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mFirst turn: Initialize adjoints and reverse first step\033[0m\n"); 70 break; 71 case 4: 72 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mForward and reverse one step\033[0m\n"); 73 break; 74 case 5: 75 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mRestore in checkpoint number %D (located in RAM)\033[0m\n",rctx->check); 76 break; 77 case 7: 78 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mStore in checkpoint number %D (located on disk)\033[0m\n",rctx->check); 79 break; 80 case 8: 81 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mRestore in checkpoint number %D (located on disk)\033[0m\n",rctx->check); 82 break; 83 case -1: 84 PetscPrintf(PETSC_COMM_WORLD,"\x1B[35mError!"); 85 break; 86 } 87 } 88 #endif 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "StackCreate" 92 static PetscErrorCode StackCreate(MPI_Comm comm,Stack *s,PetscInt size,PetscInt ny) 93 { 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 s->top = -1; 98 s->comm = comm; 99 s->numY = ny; 100 101 ierr = PetscMalloc1(size*sizeof(StackElement),&s->stack);CHKERRQ(ierr); 102 PetscFunctionReturn(0); 103 } 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "StackDestroy" 107 static PetscErrorCode StackDestroy(Stack *s) 108 { 109 PetscInt i; 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 if (s->top>-1) { 114 for (i=0;i<=s->top;i++) { 115 ierr = VecDestroy(&s->stack[i]->X);CHKERRQ(ierr); 116 if (!s->solution_only) { 117 ierr = VecDestroyVecs(s->numY,&s->stack[i]->Y);CHKERRQ(ierr); 118 } 119 ierr = PetscFree(s->stack[i]);CHKERRQ(ierr); 120 } 121 } 122 ierr = PetscFree(s->stack);CHKERRQ(ierr); 123 if (s->stype) { 124 ierr = PetscFree(s->rctx);CHKERRQ(ierr); 125 } 126 ierr = PetscFree(s);CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "StackPush" 132 static PetscErrorCode StackPush(Stack *s,StackElement e) 133 { 134 PetscFunctionBegin; 135 if (s->top+1 >= s->stacksize) SETERRQ1(s->comm,PETSC_ERR_MEMC,"Maximum stack size (%D) exceeded",s->stacksize); 136 s->stack[++s->top] = e; 137 PetscFunctionReturn(0); 138 } 139 140 #undef __FUNCT__ 141 #define __FUNCT__ "StackPop" 142 static PetscErrorCode StackPop(Stack *s,StackElement *e) 143 { 144 PetscFunctionBegin; 145 if (s->top == -1) SETERRQ(s->comm,PETSC_ERR_MEMC,"Empty stack"); 146 *e = s->stack[s->top--]; 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "StackTop" 152 static PetscErrorCode StackTop(Stack *s,StackElement *e) 153 { 154 PetscFunctionBegin; 155 *e = s->stack[s->top]; 156 PetscFunctionReturn(0); 157 } 158 159 #undef __FUNCT__ 160 #define __FUNCT__ "StackFind" 161 static PetscErrorCode StackFind(Stack *s,StackElement *e,PetscInt index) 162 { 163 PetscFunctionBegin; 164 *e = s->stack[index]; 165 PetscFunctionReturn(0); 166 } 167 168 #undef __FUNCT__ 169 #define __FUNCT__ "OutputBIN" 170 static PetscErrorCode OutputBIN(const char *filename,PetscViewer *viewer) 171 { 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 ierr = PetscViewerCreate(PETSC_COMM_WORLD,viewer);CHKERRQ(ierr); 176 ierr = PetscViewerSetType(*viewer,PETSCVIEWERBINARY);CHKERRQ(ierr); 177 ierr = PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);CHKERRQ(ierr); 178 ierr = PetscViewerFileSetName(*viewer,filename);CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "WriteToDisk" 184 static PetscErrorCode WriteToDisk(PetscInt stepnum,PetscReal time,PetscReal timeprev,Vec X,Vec *Y,PetscInt numY,PetscBool solution_only,PetscViewer viewer) 185 { 186 PetscInt i; 187 PetscErrorCode ierr; 188 189 PetscFunctionBegin; 190 ierr = PetscViewerBinaryWrite(viewer,&stepnum,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 191 ierr = VecView(X,viewer);CHKERRQ(ierr); 192 for (i=0;!solution_only && i<numY;i++) { 193 ierr = VecView(Y[i],viewer);CHKERRQ(ierr); 194 } 195 ierr = PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 196 ierr = PetscViewerBinaryWrite(viewer,&timeprev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "ReadFromDisk" 202 static PetscErrorCode ReadFromDisk(PetscInt *stepnum,PetscReal *time,PetscReal *timeprev,Vec X,Vec *Y,PetscInt numY,PetscBool solution_only,PetscViewer viewer) 203 { 204 PetscInt i; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 ierr = PetscViewerBinaryRead(viewer,stepnum,1,NULL,PETSC_INT);CHKERRQ(ierr); 209 ierr = VecLoad(X,viewer);CHKERRQ(ierr); 210 for (i=0;!solution_only && i<numY;i++) { 211 ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr); 212 } 213 ierr = PetscViewerBinaryRead(viewer,time,1,NULL,PETSC_REAL);CHKERRQ(ierr); 214 ierr = PetscViewerBinaryRead(viewer,timeprev,1,NULL,PETSC_REAL);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "StackDumpAll" 220 static PetscErrorCode StackDumpAll(TS ts,Stack *s,PetscInt id) 221 { 222 Vec *Y; 223 PetscInt i; 224 StackElement e; 225 PetscViewer viewer; 226 char filename[PETSC_MAX_PATH_LEN]; 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 if (id == 1) { 231 PetscMPIInt rank; 232 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);CHKERRQ(ierr); 233 if (!rank) { 234 ierr = PetscRMTree("SA-data");CHKERRQ(ierr); 235 ierr = PetscMkdir("SA-data");CHKERRQ(ierr); 236 } 237 } 238 ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-STACK%06d.bin",id);CHKERRQ(ierr); 239 ierr = OutputBIN(filename,&viewer);CHKERRQ(ierr); 240 for (i=0;i<s->stacksize;i++) { 241 e = s->stack[i]; 242 ierr = WriteToDisk(e->stepnum,e->time,e->timeprev,e->X,e->Y,s->numY,s->solution_only,viewer);CHKERRQ(ierr); 243 } 244 /* 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 */ 245 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 246 ierr = WriteToDisk(ts->total_steps,ts->ptime,ts->ptime_prev,ts->vec_sol,Y,s->numY,s->solution_only,viewer);CHKERRQ(ierr); 247 for (i=0;i<s->stacksize;i++) { 248 ierr = StackPop(s,&e);CHKERRQ(ierr); 249 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 250 if (!s->solution_only) { 251 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 252 } 253 ierr = PetscFree(e);CHKERRQ(ierr); 254 } 255 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 256 PetscFunctionReturn(0); 257 } 258 259 #undef __FUNCT__ 260 #define __FUNCT__ "StackLoadAll" 261 static PetscErrorCode StackLoadAll(TS ts,Stack *s,PetscInt id) 262 { 263 Vec *Y; 264 PetscInt i; 265 StackElement e; 266 PetscViewer viewer; 267 char filename[PETSC_MAX_PATH_LEN]; 268 PetscErrorCode ierr; 269 270 PetscFunctionBegin; 271 ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-STACK%06d.bin",id);CHKERRQ(ierr); 272 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 273 for (i=0;i<s->stacksize;i++) { 274 ierr = PetscCalloc1(1,&e); 275 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 276 ierr = VecDuplicate(Y[0],&e->X);CHKERRQ(ierr); 277 if (!s->solution_only && s->numY>0) { 278 ierr = VecDuplicateVecs(Y[0],s->numY,&e->Y);CHKERRQ(ierr); 279 } 280 ierr = StackPush(s,e);CHKERRQ(ierr); 281 } 282 for (i=0;i<s->stacksize;i++) { 283 e = s->stack[i]; 284 ierr = ReadFromDisk(&e->stepnum,&e->time,&e->timeprev,e->X,e->Y,s->numY,s->solution_only,viewer);CHKERRQ(ierr); 285 } 286 /* load the last step into TS */ 287 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 288 ierr = ReadFromDisk(&ts->total_steps,&ts->ptime,&ts->ptime_prev,ts->vec_sol,Y,s->numY,s->solution_only,viewer);CHKERRQ(ierr); 289 ierr = TSSetTimeStep(ts,ts->ptime-ts->ptime_prev);CHKERRQ(ierr); 290 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "DumpSingle" 296 static PetscErrorCode DumpSingle(TS ts,Stack *s,PetscInt id) 297 { 298 Vec *Y; 299 PetscInt stepnum; 300 PetscViewer viewer; 301 char filename[PETSC_MAX_PATH_LEN]; 302 PetscErrorCode ierr; 303 304 PetscFunctionBegin; 305 ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr); 306 if (id == 0) { 307 PetscMPIInt rank; 308 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);CHKERRQ(ierr); 309 if (!rank) { 310 ierr = PetscRMTree("SA-data");CHKERRQ(ierr); 311 ierr = PetscMkdir("SA-data");CHKERRQ(ierr); 312 } 313 } 314 ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-CPS%06d.bin",id);CHKERRQ(ierr); 315 ierr = OutputBIN(filename,&viewer);CHKERRQ(ierr); 316 317 ierr = PetscLogEventBegin(Disk_Write,ts,0,0,0);CHKERRQ(ierr); 318 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 319 ierr = WriteToDisk(stepnum,ts->ptime,ts->ptime_prev,ts->vec_sol,Y,s->numY,s->solution_only,viewer);CHKERRQ(ierr); 320 ierr = PetscLogEventEnd(Disk_Write,ts,0,0,0);CHKERRQ(ierr); 321 322 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 323 PetscFunctionReturn(0); 324 } 325 326 #undef __FUNCT__ 327 #define __FUNCT__ "LoadSingle" 328 static PetscErrorCode LoadSingle(TS ts,Stack *s,PetscInt id) 329 { 330 Vec *Y; 331 PetscViewer viewer; 332 char filename[PETSC_MAX_PATH_LEN]; 333 PetscErrorCode ierr; 334 335 PetscFunctionBegin; 336 ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-CPS%06d.bin",id);CHKERRQ(ierr); 337 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 338 339 ierr = PetscLogEventBegin(Disk_Read,ts,0,0,0);CHKERRQ(ierr); 340 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 341 ierr = ReadFromDisk(&ts->total_steps,&ts->ptime,&ts->ptime_prev,ts->vec_sol,Y,s->numY,s->solution_only,viewer);CHKERRQ(ierr); 342 ierr = PetscLogEventEnd(Disk_Read,ts,0,0,0);CHKERRQ(ierr); 343 344 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 345 PetscFunctionReturn(0); 346 } 347 348 #undef __FUNCT__ 349 #define __FUNCT__ "TSTrajectorySetStride_Memory" 350 PetscErrorCode TSTrajectorySetStride_Memory(TSTrajectory tj,TS ts,PetscInt stride) 351 { 352 Stack *s = (Stack*)tj->data; 353 354 PetscFunctionBegin; 355 s->stride = stride; 356 PetscFunctionReturn(0); 357 } 358 359 #undef __FUNCT__ 360 #define __FUNCT__ "TSTrajectorySetMaxCpsRAM_Memory" 361 PetscErrorCode TSTrajectorySetMaxCpsRAM_Memory(TSTrajectory tj,TS ts,PetscInt max_cps_ram) 362 { 363 Stack *s = (Stack*)tj->data; 364 365 PetscFunctionBegin; 366 s->max_cps_ram = max_cps_ram; 367 PetscFunctionReturn(0); 368 } 369 370 #undef __FUNCT__ 371 #define __FUNCT__ "TSTrajectorySetMaxCpsDisk_Memory" 372 PetscErrorCode TSTrajectorySetMaxCpsDisk_Memory(TSTrajectory tj,TS ts,PetscInt max_cps_disk) 373 { 374 Stack *s = (Stack*)tj->data; 375 376 PetscFunctionBegin; 377 s->max_cps_disk = max_cps_disk; 378 PetscFunctionReturn(0); 379 } 380 381 #undef __FUNCT__ 382 #define __FUNCT__ "TSTrajectorySetRevolveOnline" 383 PetscErrorCode TSTrajectorySetRevolveOnline(TSTrajectory tj,PetscBool use_online) 384 { 385 Stack *s = (Stack*)tj->data; 386 387 PetscFunctionBegin; 388 s->use_online = use_online; 389 PetscFunctionReturn(0); 390 } 391 392 #undef __FUNCT__ 393 #define __FUNCT__ "TSTrajectorySetFromOptions_Memory" 394 PetscErrorCode TSTrajectorySetFromOptions_Memory(PetscOptions *PetscOptionsObject,TSTrajectory tj) 395 { 396 Stack *s = (Stack*)tj->data; 397 PetscErrorCode ierr; 398 399 PetscFunctionBegin; 400 ierr = PetscOptionsHead(PetscOptionsObject,"Memory based TS trajectory options");CHKERRQ(ierr); 401 { 402 ierr = PetscOptionsInt("-tstrajectory_max_cps_ram","Maximum number of checkpoints in RAM","TSTrajectorySetMaxCpsRAM_Memory",s->max_cps_ram,&s->max_cps_ram,NULL);CHKERRQ(ierr); 403 ierr = PetscOptionsInt("-tstrajectory_max_cps_disk","Maximum number of checkpoints on disk","TSTrajectorySetMaxCpsDisk_Memory",s->max_cps_disk,&s->max_cps_disk,NULL);CHKERRQ(ierr); 404 ierr = PetscOptionsInt("-tstrajectory_stride","Stride to save checkpoints to file","TSTrajectorySetStride_Memory",s->stride,&s->stride,NULL);CHKERRQ(ierr); 405 ierr = PetscOptionsBool("-tstrajectory_revolve_online","Trick TS trajectory into using online mode of revolve","TSTrajectorySetRevolveOnline",s->use_online,&s->use_online,NULL);CHKERRQ(ierr); 406 } 407 ierr = PetscOptionsTail();CHKERRQ(ierr); 408 PetscFunctionReturn(0); 409 } 410 411 #undef __FUNCT__ 412 #define __FUNCT__ "TSTrajectorySetUp_Memory" 413 PetscErrorCode TSTrajectorySetUp_Memory(TSTrajectory tj,TS ts) 414 { 415 Stack *s = (Stack*)tj->data; 416 RevolveCTX *rctx; 417 PetscInt numY; 418 PetscBool flg; 419 PetscErrorCode ierr; 420 421 PetscFunctionBegin; 422 PetscStrcmp(((PetscObject)ts->adapt)->type_name,TSADAPTNONE,&flg); 423 if (flg) { /* fixed time step */ 424 s->total_steps = PetscMin(ts->max_steps,(PetscInt)(ceil(ts->max_time/ts->time_step))); 425 } 426 if (s->max_cps_ram > 1) s->stacksize = s->max_cps_ram; 427 if (s->stride > 1) { /* two level mode works for both fixed time step and adaptive time step */ 428 if (s->max_cps_ram > 1 && s->max_cps_ram < s->stride-1) { /* use revolve_offline for each stride */ 429 s->stype = TWO_LEVEL_REVOLVE; 430 }else { /* checkpoint all for each stride */ 431 s->stype = TWO_LEVEL_NOREVOLVE; 432 s->stacksize = s->solution_only ? s->stride : s->stride-1; 433 } 434 } else { 435 if (flg) { /* fixed time step */ 436 if (s->max_cps_ram >= s->total_steps-1 || s->max_cps_ram < 1) { /* checkpoint all */ 437 s->stype = NONE; 438 s->stacksize = s->solution_only ? s->total_steps : s->total_steps-1; 439 } else { 440 if (s->max_cps_disk > 1) { /* disk can be used */ 441 s->stype = REVOLVE_MULTISTAGE; 442 } else { /* memory only */ 443 s->stype = REVOLVE_OFFLINE; 444 } 445 } 446 } else { /* adaptive time step */ 447 s->stype = REVOLVE_ONLINE; 448 } 449 if (s->use_online) { /* trick into online */ 450 s->stype = REVOLVE_ONLINE; 451 s->stacksize = s->max_cps_ram; 452 } 453 } 454 455 if (s->stype > TWO_LEVEL_NOREVOLVE) { 456 #ifndef PETSC_HAVE_REVOLVE 457 SETERRQ(s->comm,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."); 458 #else 459 if (s->stype == TWO_LEVEL_REVOLVE) revolve_create_offline(s->stride,s->max_cps_ram); 460 else if (s->stype == REVOLVE_OFFLINE) revolve_create_offline(s->total_steps,s->max_cps_ram); 461 else if (s->stype == REVOLVE_ONLINE) revolve_create_online(s->max_cps_ram); 462 else if (s->stype ==REVOLVE_MULTISTAGE) revolve_create_multistage(s->total_steps,s->max_cps_ram+s->max_cps_disk,s->max_cps_ram); 463 464 ierr = PetscCalloc1(1,&rctx);CHKERRQ(ierr); 465 rctx->snaps_in = s->max_cps_ram; /* for theta methods snaps_in=2*max_cps_ram */ 466 rctx->reverseonestep = PETSC_FALSE; 467 rctx->check = 0; 468 rctx->oldcapo = 0; 469 rctx->capo = 0; 470 rctx->info = 2; 471 rctx->fine = (s->stride > 1) ? s->stride : s->total_steps; 472 473 s->rctx = rctx; 474 if (s->stype == REVOLVE_ONLINE) rctx->fine = -1; 475 #endif 476 } 477 478 s->recompute = PETSC_FALSE; 479 480 ierr = TSGetStages(ts,&numY,PETSC_IGNORE);CHKERRQ(ierr); 481 ierr = StackCreate(PetscObjectComm((PetscObject)ts),s,s->stacksize,numY);CHKERRQ(ierr); 482 PetscFunctionReturn(0); 483 } 484 485 #undef __FUNCT__ 486 #define __FUNCT__ "TSTrajectorySet_Memory" 487 PetscErrorCode TSTrajectorySet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X) 488 { 489 PetscInt i; 490 Vec *Y; 491 PetscReal timeprev; 492 StackElement e; 493 Stack *s = (Stack*)tj->data; 494 PetscInt localstepnum,id; 495 RevolveCTX *rctx; 496 #ifdef PETSC_HAVE_REVOLVE 497 PetscInt whattodo,shift; 498 #endif 499 PetscErrorCode ierr; 500 501 PetscFunctionBegin; 502 if (!s->recompute) { /* use global stepnum in the forward sweep */ 503 ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr); 504 } 505 /* for consistency */ 506 if (!s->recompute && stepnum == 0) ts->ptime_prev = ts->ptime-ts->time_step; 507 508 if (s->stype == TWO_LEVEL_REVOLVE) { /* two-level mode */ 509 localstepnum = stepnum%s->stride; 510 if (stepnum!=s->total_steps && localstepnum==0 && !s->recompute) { 511 #ifdef PETSC_HAVE_REVOLVE 512 PetscPrintf(PETSC_COMM_WORLD,"\x1B[33mDump stack to file\033[0m\n"); 513 #endif 514 id = stepnum/s->stride; 515 if (s->save_stack) { 516 ierr = StackDumpAll(ts,s,id);CHKERRQ(ierr); 517 s->top = -1; /* reset top */ 518 #ifdef PETSC_HAVE_REVOLVE 519 if (s->stype == TWO_LEVEL_REVOLVE) { 520 revolve_reset(); 521 revolve_create_offline(s->stride,s->max_cps_ram); 522 rctx = s->rctx; 523 rctx->check = 0; 524 rctx->capo = 0; 525 rctx->fine = s->stride; 526 } 527 #endif 528 } else { 529 ierr = DumpSingle(ts,s,id);CHKERRQ(ierr); 530 } 531 } 532 /* first forward sweep only checkpoints once in each stride */ 533 if (!s->recompute) PetscFunctionReturn(0); 534 } else if (s->stype == TWO_LEVEL_NOREVOLVE) { 535 localstepnum = stepnum%s->stride; 536 if (stepnum != s->total_steps && stepnum != 0 && localstepnum==0 && !s->recompute) { 537 id = stepnum/s->stride; 538 if (s->save_stack) { 539 ierr = StackDumpAll(ts,s,id);CHKERRQ(ierr); 540 s->top = -1; /* reset top */ 541 } else { 542 ierr = DumpSingle(ts,s,id);CHKERRQ(ierr); 543 } 544 } 545 } else { 546 localstepnum = stepnum; 547 } 548 549 if (s->stype > TWO_LEVEL_NOREVOLVE) { 550 #ifdef PETSC_HAVE_REVOLVE 551 rctx = s->rctx; 552 if (rctx->reverseonestep && stepnum==s->total_steps) { /* intermediate information is ready inside TS, this happens at last time step */ 553 rctx->reverseonestep = PETSC_FALSE; //? 554 PetscFunctionReturn(0); 555 } 556 if (rctx->stepsleft != 0) { /* advance the solution without checkpointing anything as Revolve requires */ 557 rctx->stepsleft--; 558 PetscFunctionReturn(0); 559 } 560 561 /* let Revolve determine what to do next */ 562 shift = stepnum-localstepnum; 563 rctx->capo = localstepnum; 564 rctx->oldcapo = rctx->capo; 565 whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); 566 if (s->stype == REVOLVE_ONLINE && whattodo ==7) whattodo = 2; 567 printwhattodo(whattodo,rctx,shift); 568 if (whattodo == -1) SETERRQ(s->comm,PETSC_ERR_LIB,"Error in the Revolve library"); 569 if (whattodo == 1) { /* advance some time steps */ 570 if (s->stype == REVOLVE_ONLINE && rctx->capo >= s->total_steps-1) { 571 revolve_turn(s->total_steps,&rctx->capo,&rctx->fine); 572 printwhattodo(whattodo,rctx,shift); 573 } 574 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 575 } 576 if (whattodo == 3 || whattodo == 4) { /* ready for a reverse step */ 577 rctx->reverseonestep = PETSC_TRUE; 578 } 579 if (whattodo == 5) { /* restore a checkpoint and ask Revolve what to do next */ 580 rctx->oldcapo = rctx->capo; 581 whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); /* must return 1*/ 582 printwhattodo(whattodo,rctx,shift); 583 rctx->stepsleft = rctx->capo-rctx->oldcapo; 584 } 585 if (whattodo == 7) { /* save the checkpoint to disk */ 586 ierr = DumpSingle(ts,s,rctx->check);CHKERRQ(ierr); 587 rctx->oldcapo = rctx->capo; 588 whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); /* must return 1*/ 589 printwhattodo(whattodo,rctx,shift); 590 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 591 } 592 if (whattodo != 2) { 593 PetscFunctionReturn(0); 594 } else { /* store a checkpoint to RAM and ask Revolve how many time steps to advance next */ 595 rctx->oldcapo = rctx->capo; 596 whattodo = revolve_action(&rctx->check,&rctx->capo,&rctx->fine,rctx->snaps_in,&rctx->info,&rctx->where); /* must return 1*/ 597 printwhattodo(whattodo,rctx,shift); 598 if (s->stype == REVOLVE_ONLINE && rctx->capo >= s->total_steps-1) { 599 revolve_turn(s->total_steps,&rctx->capo,&rctx->fine); 600 printwhattodo(whattodo,rctx,shift); 601 } 602 rctx->stepsleft = rctx->capo-rctx->oldcapo-1; 603 } 604 #endif 605 } else { /* Revolve is not used */ 606 if (s->solution_only) { 607 /* skip the last two steps of each stride or the whole interval */ 608 if (stepnum >= s->total_steps-1 || s->recompute) PetscFunctionReturn(0); //? 609 } else { 610 /* skip the first and the last steps of each stride or the whole interval */ 611 if (localstepnum == 0 || stepnum == s->total_steps) PetscFunctionReturn(0); 612 } 613 } 614 615 /* checkpoint to memory, rctx->check gives the index in the stack */ 616 if (s->stype == REVOLVE_ONLINE && rctx->check != s->top+1) { /* overwrite some non-top checkpoint in the stack*/ 617 ierr = StackFind(s,&e,rctx->check);CHKERRQ(ierr); 618 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 619 e->stepnum = stepnum; 620 e->time = time; 621 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 622 e->timeprev = timeprev; 623 } else { /* push to stack */ 624 if (localstepnum < s->top) SETERRQ(s->comm,PETSC_ERR_MEMC,"Illegal modification of a non-top stack element"); 625 ierr = PetscCalloc1(1,&e);CHKERRQ(ierr); 626 ierr = VecDuplicate(X,&e->X);CHKERRQ(ierr); 627 ierr = VecCopy(X,e->X);CHKERRQ(ierr); 628 if (s->numY > 0 && !s->solution_only) { 629 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 630 ierr = VecDuplicateVecs(Y[0],s->numY,&e->Y);CHKERRQ(ierr); 631 for (i=0;i<s->numY;i++) { 632 ierr = VecCopy(Y[i],e->Y[i]);CHKERRQ(ierr); 633 } 634 } 635 e->stepnum = stepnum; 636 e->time = time; 637 /* for consistency */ 638 if (stepnum == 0) { 639 e->timeprev = e->time - ts->time_step; 640 } else { 641 ierr = TSGetPrevTime(ts,&timeprev);CHKERRQ(ierr); 642 e->timeprev = timeprev; 643 } 644 ierr = StackPush(s,e);CHKERRQ(ierr); 645 } 646 647 PetscFunctionReturn(0); 648 } 649 650 #undef __FUNCT__ 651 #define __FUNCT__ "TSTrajectoryGet_Memory" 652 PetscErrorCode TSTrajectoryGet_Memory(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t) 653 { 654 Vec *Y; 655 PetscInt i; 656 StackElement e = NULL; 657 Stack *s = (Stack*)tj->data; 658 PetscReal stepsize; 659 PetscInt localstepnum,id,steps; 660 #ifdef PETSC_HAVE_REVOLVE 661 PetscInt whattodo,shift; 662 #endif 663 PetscErrorCode ierr; 664 665 PetscFunctionBegin; 666 ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr); /* stepnum should be larger than 0 */ 667 if (s->stype == TWO_LEVEL_REVOLVE) { /* two-level mode */ 668 localstepnum = stepnum%s->stride; 669 if (localstepnum == 0) { 670 #ifdef PETSC_HAVE_REVOLVE 671 PetscPrintf(PETSC_COMM_WORLD,"\x1B[33mLoad stack from file\033[0m\n"); 672 #endif 673 if (s->save_stack) { 674 id = stepnum/s->stride; 675 ierr = StackLoadAll(ts,s,id);CHKERRQ(ierr); 676 s->top = s->stacksize-1; 677 } else { 678 /* save ts context */ 679 steps = ts->steps; 680 id = stepnum/s->stride-1; 681 ierr = LoadSingle(ts,s,id);CHKERRQ(ierr); 682 ts->steps = ts->total_steps; 683 s->recompute = PETSC_TRUE; 684 } 685 686 #ifdef PETSC_HAVE_REVOLVE 687 revolve_reset(); 688 revolve_create_offline(s->stride,s->max_cps_ram); 689 s->rctx->check = 0; 690 s->rctx->capo = 0; 691 s->rctx->fine = s->stride; 692 if (!s->solution_only) { 693 whattodo = 0; 694 while(whattodo!=3) { /* stupid revolve */ 695 whattodo = revolve_action(&s->rctx->check,&s->rctx->capo,&s->rctx->fine,s->rctx->snaps_in,&s->rctx->info,&s->rctx->where); 696 } 697 } 698 #endif 699 for (i=ts->steps;i<stepnum;i++) { /* assume fixed step size */ 700 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 701 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 702 ierr = TSStep(ts);CHKERRQ(ierr); 703 if (ts->event) { 704 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 705 } 706 if (!ts->steprollback) { 707 ierr = TSPostStep(ts);CHKERRQ(ierr); 708 } 709 } 710 ts->steps = steps; 711 ts->total_steps = stepnum; 712 } 713 } else if (s->stype == TWO_LEVEL_NOREVOLVE) { 714 localstepnum = stepnum%s->stride; 715 if (stepnum != s->total_steps && localstepnum==0) { 716 id = stepnum/s->stride; 717 if (s->save_stack) { 718 ierr = StackLoadAll(ts,s,id);CHKERRQ(ierr); 719 } else { 720 ierr = LoadSingle(ts,s,id);CHKERRQ(ierr); 721 } 722 } 723 } else { 724 localstepnum = stepnum; 725 } 726 #ifdef PETSC_HAVE_REVOLVE 727 if (s->stype > TWO_LEVEL_NOREVOLVE && s->rctx->reverseonestep) { /* ready for the reverse step */ 728 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 729 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); 730 s->rctx->reverseonestep = PETSC_FALSE; 731 PetscFunctionReturn(0); 732 } 733 #endif 734 if ((!s->solution_only && localstepnum == 0) || stepnum == s->total_steps) { 735 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 736 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); 737 s->rctx->reverseonestep = PETSC_FALSE;//? 738 PetscFunctionReturn(0); 739 } 740 741 /* recomputation is needed */ 742 steps = ts->steps; /* save the TS context */ 743 s->recompute = PETSC_TRUE; 744 745 if (s->stype > TWO_LEVEL_NOREVOLVE) { 746 #ifdef PETSC_HAVE_REVOLVE 747 s->rctx->capo = localstepnum; 748 shift = stepnum-localstepnum; 749 whattodo = revolve_action(&s->rctx->check,&s->rctx->capo,&s->rctx->fine,s->rctx->snaps_in,&s->rctx->info,&s->rctx->where); 750 if (s->stype == REVOLVE_ONLINE && whattodo == 8) whattodo = 5; 751 printwhattodo(whattodo,s->rctx,shift); 752 #endif 753 if (s->stype == REVOLVE_MULTISTAGE && !s->rctx->where) { 754 ierr = LoadSingle(ts,s,s->rctx->check);CHKERRQ(ierr); 755 ts->steps = ts->total_steps; 756 } 757 } 758 /* restore a checkpoint */ 759 if (!(s->stype == REVOLVE_MULTISTAGE && !s->rctx->where)) { 760 if (s->stype == REVOLVE_ONLINE) { 761 ierr = StackFind(s,&e,s->rctx->check);CHKERRQ(ierr); 762 } else { 763 ierr = StackTop(s,&e);CHKERRQ(ierr); 764 } 765 ierr = VecCopy(e->X,ts->vec_sol);CHKERRQ(ierr); 766 if (!s->solution_only) { 767 ierr = TSGetStages(ts,&s->numY,&Y);CHKERRQ(ierr); 768 for (i=0;i<s->numY;i++) { 769 ierr = VecCopy(e->Y[i],Y[i]);CHKERRQ(ierr); 770 } 771 } 772 *t = e->time; 773 ierr = TSSetTimeStep(ts,(*t)-e->timeprev);CHKERRQ(ierr); 774 /* reset ts context */ 775 ts->steps = e->stepnum; /* global stepnum */ 776 ts->ptime = e->time; 777 ts->ptime_prev = e->timeprev; 778 } 779 for (i=ts->steps;i<stepnum;i++) { /* assume fixed step size */ 780 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 781 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 782 ierr = TSStep(ts);CHKERRQ(ierr); 783 if (ts->event) { 784 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 785 } 786 if (!ts->steprollback) { 787 ierr = TSPostStep(ts);CHKERRQ(ierr); 788 } 789 } 790 /* reverseonestep must be true after the for loop */ 791 ts->steps = steps; 792 ts->total_steps = stepnum; 793 ierr = TSGetTimeStep(ts,&stepsize);CHKERRQ(ierr); 794 ierr = TSSetTimeStep(ts,-stepsize);CHKERRQ(ierr); 795 if (e && e->stepnum >= stepnum) { 796 SETERRQ2(s->comm,PETSC_ERR_ARG_OUTOFRANGE,"The current step no. is %D, but the step number at top of the stack is %D",stepnum,e->stepnum); 797 } 798 if (s->stype != REVOLVE_ONLINE && e && stepnum-e->stepnum==1) { /* the restored checkpoint can be deleted now */ 799 ierr = StackPop(s,&e);CHKERRQ(ierr); 800 ierr = VecDestroy(&e->X);CHKERRQ(ierr); 801 if (!s->solution_only) { 802 ierr = VecDestroyVecs(s->numY,&e->Y);CHKERRQ(ierr); 803 } 804 ierr = PetscFree(e);CHKERRQ(ierr); 805 } 806 if (s->stype > TWO_LEVEL_NOREVOLVE) { 807 #ifdef PETSC_HAVE_REVOLVE 808 s->rctx->reverseonestep = PETSC_FALSE; 809 #endif 810 } 811 PetscFunctionReturn(0); 812 } 813 814 #undef __FUNCT__ 815 #define __FUNCT__ "TSTrajectoryDestroy_Memory" 816 PETSC_EXTERN PetscErrorCode TSTrajectoryDestroy_Memory(TSTrajectory tj) 817 { 818 Stack *s = (Stack*)tj->data; 819 PetscErrorCode ierr; 820 821 PetscFunctionBegin; 822 if (s->stype > TWO_LEVEL_NOREVOLVE) { 823 #ifdef PETSC_HAVE_REVOLVE 824 revolve_reset(); 825 #endif 826 } 827 ierr = StackDestroy(s);CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830 831 /*MC 832 TSTRAJECTORYMEMORY - Stores each solution of the ODE/ADE in memory 833 834 Level: intermediate 835 836 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType() 837 838 M*/ 839 #undef __FUNCT__ 840 #define __FUNCT__ "TSTrajectoryCreate_Memory" 841 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory tj,TS ts) 842 { 843 Stack *s; 844 PetscErrorCode ierr; 845 846 PetscFunctionBegin; 847 tj->ops->set = TSTrajectorySet_Memory; 848 tj->ops->get = TSTrajectoryGet_Memory; 849 tj->ops->setup = TSTrajectorySetUp_Memory; 850 tj->ops->destroy = TSTrajectoryDestroy_Memory; 851 tj->ops->setfromoptions = TSTrajectorySetFromOptions_Memory; 852 853 ierr = PetscCalloc1(1,&s);CHKERRQ(ierr); 854 s->stype = NONE; 855 s->max_cps_ram = -1; /* -1 indicates that it is not set */ 856 s->max_cps_disk = -1; /* -1 indicates that it is not set */ 857 s->stride = 0; /* if not zero, two-level checkpointing will be used */ 858 s->use_online = PETSC_FALSE; 859 s->save_stack = PETSC_FALSE; 860 s->solution_only= PETSC_TRUE; 861 tj->data = s; 862 PetscFunctionReturn(0); 863 } 864