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