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