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