Actual source code: traj.c

  1: #include <petsc/private/tsimpl.h>
  2: #include <petsc/private/tshistoryimpl.h>
  3: #include <petscdm.h>

  5: PetscFunctionList TSTrajectoryList              = NULL;
  6: PetscBool         TSTrajectoryRegisterAllCalled = PETSC_FALSE;
  7: PetscClassId      TSTRAJECTORY_CLASSID;
  8: PetscLogEvent     TSTrajectory_Set, TSTrajectory_Get, TSTrajectory_GetVecs, TSTrajectory_SetUp;

 10: /*@C
 11:   TSTrajectoryRegister - Adds a way of storing trajectories to the `TS` package

 13:   Not Collective

 15:   Input Parameters:
 16: + sname    - the name of a new user-defined creation routine
 17: - function - the creation routine itself

 19:   Level: developer

 21:   Note:
 22:   `TSTrajectoryRegister()` may be called multiple times to add several user-defined tses.

 24: .seealso: [](ch_ts), `TSTrajectoryRegisterAll()`
 25: @*/
 26: PetscErrorCode TSTrajectoryRegister(const char sname[], PetscErrorCode (*function)(TSTrajectory, TS))
 27: {
 28:   PetscFunctionBegin;
 29:   PetscCall(PetscFunctionListAdd(&TSTrajectoryList, sname, function));
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: /*@
 34:   TSTrajectorySet - Sets a vector of state in the trajectory object

 36:   Collective

 38:   Input Parameters:
 39: + tj      - the trajectory object
 40: . ts      - the time stepper object (optional)
 41: . stepnum - the step number
 42: . time    - the current time
 43: - X       - the current solution

 45:   Level: developer

 47:   Note:
 48:   Usually one does not call this routine, it is called automatically during `TSSolve()`

 50: .seealso: [](ch_ts), `TSTrajectorySetUp()`, `TSTrajectoryDestroy()`, `TSTrajectorySetType()`, `TSTrajectorySetVariableNames()`, `TSGetTrajectory()`, `TSTrajectoryGet()`, `TSTrajectoryGetVecs()`
 51: @*/
 52: PetscErrorCode TSTrajectorySet(TSTrajectory tj, TS ts, PetscInt stepnum, PetscReal time, Vec X)
 53: {
 54:   PetscFunctionBegin;
 55:   if (!tj) PetscFunctionReturn(PETSC_SUCCESS);
 61:   PetscCheck(tj->setupcalled, PetscObjectComm((PetscObject)tj), PETSC_ERR_ORDER, "TSTrajectorySetUp should be called first");
 62:   if (tj->monitor) PetscCall(PetscViewerASCIIPrintf(tj->monitor, "TSTrajectorySet: stepnum %" PetscInt_FMT ", time %g (stages %" PetscInt_FMT ")\n", stepnum, (double)time, (PetscInt)!tj->solution_only));
 63:   PetscCall(PetscLogEventBegin(TSTrajectory_Set, tj, ts, 0, 0));
 64:   PetscUseTypeMethod(tj, set, ts, stepnum, time, X);
 65:   PetscCall(PetscLogEventEnd(TSTrajectory_Set, tj, ts, 0, 0));
 66:   if (tj->usehistory) PetscCall(TSHistoryUpdate(tj->tsh, stepnum, time));
 67:   if (tj->lag.caching) tj->lag.Udotcached.time = PETSC_MIN_REAL;
 68:   PetscFunctionReturn(PETSC_SUCCESS);
 69: }

 71: /*@
 72:   TSTrajectoryGetNumSteps - Return the number of steps registered in the `TSTrajectory` via `TSTrajectorySet()`.

 74:   Not Collective.

 76:   Input Parameter:
 77: . tj - the trajectory object

 79:   Output Parameter:
 80: . steps - the number of steps

 82:   Level: developer

 84: .seealso: [](ch_ts), `TS`, `TSTrajectorySet()`
 85: @*/
 86: PetscErrorCode TSTrajectoryGetNumSteps(TSTrajectory tj, PetscInt *steps)
 87: {
 88:   PetscFunctionBegin;
 90:   PetscAssertPointer(steps, 2);
 91:   PetscCall(TSHistoryGetNumSteps(tj->tsh, steps));
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: /*@
 96:   TSTrajectoryGet - Updates the solution vector of a time stepper object by querying the `TSTrajectory`

 98:   Collective

100:   Input Parameters:
101: + tj      - the trajectory object
102: . ts      - the time stepper object
103: - stepnum - the step number

105:   Output Parameter:
106: . time - the time associated with the step number

108:   Level: developer

110:   Note:
111:   Usually one does not call this routine, it is called automatically during `TSSolve()`

113: .seealso: [](ch_ts), `TS`, `TSSolve()`, `TSTrajectorySetUp()`, `TSTrajectoryDestroy()`, `TSTrajectorySetType()`, `TSTrajectorySetVariableNames()`, `TSGetTrajectory()`, `TSTrajectorySet()`, `TSTrajectoryGetVecs()`, `TSGetSolution()`
114: @*/
115: PetscErrorCode TSTrajectoryGet(TSTrajectory tj, TS ts, PetscInt stepnum, PetscReal *time)
116: {
117:   PetscFunctionBegin;
118:   PetscCheck(tj, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TS solver did not save trajectory");
122:   PetscAssertPointer(time, 4);
123:   PetscCheck(tj->setupcalled, PetscObjectComm((PetscObject)tj), PETSC_ERR_ORDER, "TSTrajectorySetUp should be called first");
124:   PetscCheck(stepnum >= 0, PetscObjectComm((PetscObject)tj), PETSC_ERR_PLIB, "Requesting negative step number");
125:   if (tj->monitor) {
126:     PetscCall(PetscViewerASCIIPrintf(tj->monitor, "TSTrajectoryGet: stepnum %" PetscInt_FMT ", stages %" PetscInt_FMT "\n", stepnum, (PetscInt)!tj->solution_only));
127:     PetscCall(PetscViewerFlush(tj->monitor));
128:   }
129:   PetscCall(PetscLogEventBegin(TSTrajectory_Get, tj, ts, 0, 0));
130:   PetscUseTypeMethod(tj, get, ts, stepnum, time);
131:   PetscCall(PetscLogEventEnd(TSTrajectory_Get, tj, ts, 0, 0));
132:   PetscFunctionReturn(PETSC_SUCCESS);
133: }

135: /*@
136:   TSTrajectoryGetVecs - Reconstructs the vector of state and its time derivative using information from the `TSTrajectory` and, possibly, from the `TS`

138:   Collective

140:   Input Parameters:
141: + tj      - the trajectory object
142: . ts      - the time stepper object (optional)
143: - stepnum - the requested step number

145:   Output Parameters:
146: + time - On input time for the step if step number is `PETSC_DECIDE`, on output the time associated with the step number
147: . U    - state vector (can be `NULL`)
148: - Udot - time derivative of state vector (can be `NULL`)

150:   Level: developer

152:   Notes:
153:   If the step number is `PETSC_DECIDE`, the time argument is used to inquire the trajectory.
154:   If the requested time does not match any in the trajectory, Lagrangian interpolations are returned.

156: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSTrajectorySetUp()`, `TSTrajectoryDestroy()`, `TSTrajectorySetType()`, `TSTrajectorySetVariableNames()`, `TSGetTrajectory()`, `TSTrajectorySet()`, `TSTrajectoryGet()`
157: @*/
158: PetscErrorCode TSTrajectoryGetVecs(TSTrajectory tj, TS ts, PetscInt stepnum, PetscReal *time, Vec U, Vec Udot)
159: {
160:   PetscFunctionBegin;
161:   PetscCheck(tj, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "TS solver did not save trajectory");
165:   PetscAssertPointer(time, 4);
168:   if (!U && !Udot) PetscFunctionReturn(PETSC_SUCCESS);
169:   PetscCheck(tj->setupcalled, PetscObjectComm((PetscObject)tj), PETSC_ERR_ORDER, "TSTrajectorySetUp should be called first");
170:   PetscCall(PetscLogEventBegin(TSTrajectory_GetVecs, tj, ts, 0, 0));
171:   if (tj->monitor) {
172:     PetscInt pU, pUdot;
173:     pU    = U ? 1 : 0;
174:     pUdot = Udot ? 1 : 0;
175:     PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Requested by GetVecs %" PetscInt_FMT " %" PetscInt_FMT ": stepnum %" PetscInt_FMT ", time %g\n", pU, pUdot, stepnum, (double)*time));
176:     PetscCall(PetscViewerFlush(tj->monitor));
177:   }
178:   if (U && tj->lag.caching) {
179:     PetscObjectId    id;
180:     PetscObjectState state;

182:     PetscCall(PetscObjectStateGet((PetscObject)U, &state));
183:     PetscCall(PetscObjectGetId((PetscObject)U, &id));
184:     if (stepnum == PETSC_DECIDE) {
185:       if (id == tj->lag.Ucached.id && *time == tj->lag.Ucached.time && state == tj->lag.Ucached.state) U = NULL;
186:     } else {
187:       if (id == tj->lag.Ucached.id && stepnum == tj->lag.Ucached.step && state == tj->lag.Ucached.state) U = NULL;
188:     }
189:     if (tj->monitor && !U) {
190:       PetscCall(PetscViewerASCIIPushTab(tj->monitor));
191:       PetscCall(PetscViewerASCIIPrintf(tj->monitor, "State vector cached\n"));
192:       PetscCall(PetscViewerASCIIPopTab(tj->monitor));
193:       PetscCall(PetscViewerFlush(tj->monitor));
194:     }
195:   }
196:   if (Udot && tj->lag.caching) {
197:     PetscObjectId    id;
198:     PetscObjectState state;

200:     PetscCall(PetscObjectStateGet((PetscObject)Udot, &state));
201:     PetscCall(PetscObjectGetId((PetscObject)Udot, &id));
202:     if (stepnum == PETSC_DECIDE) {
203:       if (id == tj->lag.Udotcached.id && *time == tj->lag.Udotcached.time && state == tj->lag.Udotcached.state) Udot = NULL;
204:     } else {
205:       if (id == tj->lag.Udotcached.id && stepnum == tj->lag.Udotcached.step && state == tj->lag.Udotcached.state) Udot = NULL;
206:     }
207:     if (tj->monitor && !Udot) {
208:       PetscCall(PetscViewerASCIIPushTab(tj->monitor));
209:       PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Derivative vector cached\n"));
210:       PetscCall(PetscViewerASCIIPopTab(tj->monitor));
211:       PetscCall(PetscViewerFlush(tj->monitor));
212:     }
213:   }
214:   if (!U && !Udot) {
215:     PetscCall(PetscLogEventEnd(TSTrajectory_GetVecs, tj, ts, 0, 0));
216:     PetscFunctionReturn(PETSC_SUCCESS);
217:   }

219:   if (stepnum == PETSC_DECIDE || Udot) { /* reverse search for requested time in TSHistory */
220:     if (tj->monitor) PetscCall(PetscViewerASCIIPushTab(tj->monitor));
221:     /* cached states will be updated in the function */
222:     PetscCall(TSTrajectoryReconstruct_Private(tj, ts, *time, U, Udot));
223:     if (tj->monitor) {
224:       PetscCall(PetscViewerASCIIPopTab(tj->monitor));
225:       PetscCall(PetscViewerFlush(tj->monitor));
226:     }
227:   } else if (U) { /* we were asked to load from stepnum, use TSTrajectoryGet */
228:     TS  fakets = ts;
229:     Vec U2;

231:     /* use a fake TS if ts is missing */
232:     if (!ts) {
233:       PetscCall(PetscObjectQuery((PetscObject)tj, "__fake_ts", (PetscObject *)&fakets));
234:       if (!fakets) {
235:         PetscCall(TSCreate(PetscObjectComm((PetscObject)tj), &fakets));
236:         PetscCall(PetscObjectCompose((PetscObject)tj, "__fake_ts", (PetscObject)fakets));
237:         PetscCall(PetscObjectDereference((PetscObject)fakets));
238:         PetscCall(VecDuplicate(U, &U2));
239:         PetscCall(TSSetSolution(fakets, U2));
240:         PetscCall(PetscObjectDereference((PetscObject)U2));
241:       }
242:     }
243:     PetscCall(TSTrajectoryGet(tj, fakets, stepnum, time));
244:     PetscCall(TSGetSolution(fakets, &U2));
245:     PetscCall(VecCopy(U2, U));
246:     PetscCall(PetscObjectStateGet((PetscObject)U, &tj->lag.Ucached.state));
247:     PetscCall(PetscObjectGetId((PetscObject)U, &tj->lag.Ucached.id));
248:     tj->lag.Ucached.time = *time;
249:     tj->lag.Ucached.step = stepnum;
250:   }
251:   PetscCall(PetscLogEventEnd(TSTrajectory_GetVecs, tj, ts, 0, 0));
252:   PetscFunctionReturn(PETSC_SUCCESS);
253: }

255: /*@C
256:   TSTrajectoryViewFromOptions - View a `TSTrajectory` based on values in the options database

258:   Collective

260:   Input Parameters:
261: + A    - the `TSTrajectory` context
262: . obj  - Optional object that provides prefix used for option name
263: - name - command line option

265:   Level: intermediate

267: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryView`, `PetscObjectViewFromOptions()`, `TSTrajectoryCreate()`
268: @*/
269: PetscErrorCode TSTrajectoryViewFromOptions(TSTrajectory A, PetscObject obj, const char name[])
270: {
271:   PetscFunctionBegin;
273:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: /*@C
278:   TSTrajectoryView - Prints information about the trajectory object

280:   Collective

282:   Input Parameters:
283: + tj     - the `TSTrajectory` context obtained from `TSTrajectoryCreate()`
284: - viewer - visualization context

286:   Options Database Key:
287: . -ts_trajectory_view - calls `TSTrajectoryView()` at end of `TSAdjointStep()`

289:   Level: developer

291:   Notes:
292:   The available visualization contexts include
293: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
294: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
295:   output where only the first processor opens
296:   the file.  All other processors send their
297:   data to the first processor to print.

299:   The user can open an alternative visualization context with
300:   `PetscViewerASCIIOpen()` - output to a specified file.

302: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `PetscViewer`, `PetscViewerASCIIOpen()`
303: @*/
304: PetscErrorCode TSTrajectoryView(TSTrajectory tj, PetscViewer viewer)
305: {
306:   PetscBool iascii;

308:   PetscFunctionBegin;
310:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tj), &viewer));
312:   PetscCheckSameComm(tj, 1, viewer, 2);

314:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
315:   if (iascii) {
316:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tj, viewer));
317:     PetscCall(PetscViewerASCIIPrintf(viewer, "  total number of recomputations for adjoint calculation = %" PetscInt_FMT "\n", tj->recomps));
318:     PetscCall(PetscViewerASCIIPrintf(viewer, "  disk checkpoint reads = %" PetscInt_FMT "\n", tj->diskreads));
319:     PetscCall(PetscViewerASCIIPrintf(viewer, "  disk checkpoint writes = %" PetscInt_FMT "\n", tj->diskwrites));
320:     PetscCall(PetscViewerASCIIPushTab(viewer));
321:     PetscTryTypeMethod(tj, view, viewer);
322:     PetscCall(PetscViewerASCIIPopTab(viewer));
323:   }
324:   PetscFunctionReturn(PETSC_SUCCESS);
325: }

327: /*@C
328:   TSTrajectorySetVariableNames - Sets the name of each component in the solution vector so that it may be saved with the trajectory

330:   Collective

332:   Input Parameters:
333: + ctx   - the trajectory context
334: - names - the names of the components, final string must be `NULL`

336:   Level: intermediate

338:   Fortran Notes:
339:   Fortran interface is not possible because of the string array argument

341: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`
342: @*/
343: PetscErrorCode TSTrajectorySetVariableNames(TSTrajectory ctx, const char *const *names)
344: {
345:   PetscFunctionBegin;
347:   PetscAssertPointer(names, 2);
348:   PetscCall(PetscStrArrayDestroy(&ctx->names));
349:   PetscCall(PetscStrArrayallocpy(names, &ctx->names));
350:   PetscFunctionReturn(PETSC_SUCCESS);
351: }

353: /*@C
354:   TSTrajectorySetTransform - Solution vector will be transformed by provided function before being saved to disk

356:   Collective

358:   Input Parameters:
359: + tj        - the `TSTrajectory` context
360: . transform - the transform function
361: . destroy   - function to destroy the optional context
362: - tctx      - optional context used by transform function

364:   Level: intermediate

366: .seealso: [](ch_ts), `TSTrajectorySetVariableNames()`, `TSTrajectory`, `TSMonitorLGSetTransform()`
367: @*/
368: PetscErrorCode TSTrajectorySetTransform(TSTrajectory tj, PetscErrorCode (*transform)(void *, Vec, Vec *), PetscErrorCode (*destroy)(void *), void *tctx)
369: {
370:   PetscFunctionBegin;
372:   tj->transform        = transform;
373:   tj->transformdestroy = destroy;
374:   tj->transformctx     = tctx;
375:   PetscFunctionReturn(PETSC_SUCCESS);
376: }

378: /*@
379:   TSTrajectoryCreate - This function creates an empty trajectory object used to store the time dependent solution of an ODE/DAE

381:   Collective

383:   Input Parameter:
384: . comm - the communicator

386:   Output Parameter:
387: . tj - the trajectory object

389:   Level: developer

391:   Notes:
392:   Usually one does not call this routine, it is called automatically when one calls `TSSetSaveTrajectory()`.

394: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSTrajectorySetUp()`, `TSTrajectoryDestroy()`, `TSTrajectorySetType()`, `TSTrajectorySetVariableNames()`, `TSGetTrajectory()`, `TSTrajectorySetKeepFiles()`
395: @*/
396: PetscErrorCode TSTrajectoryCreate(MPI_Comm comm, TSTrajectory *tj)
397: {
398:   TSTrajectory t;

400:   PetscFunctionBegin;
401:   PetscAssertPointer(tj, 2);
402:   *tj = NULL;
403:   PetscCall(TSInitializePackage());

405:   PetscCall(PetscHeaderCreate(t, TSTRAJECTORY_CLASSID, "TSTrajectory", "Time stepping", "TS", comm, TSTrajectoryDestroy, TSTrajectoryView));
406:   t->setupcalled = PETSC_FALSE;
407:   PetscCall(TSHistoryCreate(comm, &t->tsh));

409:   t->lag.order            = 1;
410:   t->lag.L                = NULL;
411:   t->lag.T                = NULL;
412:   t->lag.W                = NULL;
413:   t->lag.WW               = NULL;
414:   t->lag.TW               = NULL;
415:   t->lag.TT               = NULL;
416:   t->lag.caching          = PETSC_TRUE;
417:   t->lag.Ucached.id       = 0;
418:   t->lag.Ucached.state    = -1;
419:   t->lag.Ucached.time     = PETSC_MIN_REAL;
420:   t->lag.Ucached.step     = PETSC_MAX_INT;
421:   t->lag.Udotcached.id    = 0;
422:   t->lag.Udotcached.state = -1;
423:   t->lag.Udotcached.time  = PETSC_MIN_REAL;
424:   t->lag.Udotcached.step  = PETSC_MAX_INT;
425:   t->adjoint_solve_mode   = PETSC_TRUE;
426:   t->solution_only        = PETSC_FALSE;
427:   t->keepfiles            = PETSC_FALSE;
428:   t->usehistory           = PETSC_TRUE;
429:   *tj                     = t;
430:   PetscCall(TSTrajectorySetFiletemplate(t, "TS-%06" PetscInt_FMT ".bin"));
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

434: /*@C
435:   TSTrajectorySetType - Sets the storage method to be used as in a trajectory

437:   Collective

439:   Input Parameters:
440: + tj   - the `TSTrajectory` context
441: . ts   - the `TS` context
442: - type - a known method

444:   Options Database Key:
445: . -ts_trajectory_type <type> - Sets the method; use -help for a list of available methods (for instance, basic)

447:   Level: developer

449:   Developer Notes:
450:   Why does this option require access to the `TS`

452: .seealso: [](ch_ts), `TSTrajectory`, `TS`, `TSTrajectoryCreate()`, `TSTrajectorySetFromOptions()`, `TSTrajectoryDestroy()`, `TSTrajectoryGetType()`
453: @*/
454: PetscErrorCode TSTrajectorySetType(TSTrajectory tj, TS ts, TSTrajectoryType type)
455: {
456:   PetscErrorCode (*r)(TSTrajectory, TS);
457:   PetscBool match;

459:   PetscFunctionBegin;
461:   PetscCall(PetscObjectTypeCompare((PetscObject)tj, type, &match));
462:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

464:   PetscCall(PetscFunctionListFind(TSTrajectoryList, type, &r));
465:   PetscCheck(r, PetscObjectComm((PetscObject)tj), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSTrajectory type: %s", type);
466:   PetscTryTypeMethod(tj, destroy);
467:   tj->ops->destroy = NULL;
468:   PetscCall(PetscMemzero(tj->ops, sizeof(*tj->ops)));

470:   PetscCall(PetscObjectChangeTypeName((PetscObject)tj, type));
471:   PetscCall((*r)(tj, ts));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: /*@C
476:   TSTrajectoryGetType - Gets the trajectory type

478:   Collective

480:   Input Parameters:
481: + tj - the `TSTrajectory` context
482: - ts - the `TS` context

484:   Output Parameter:
485: . type - a known method

487:   Level: developer

489: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectorySetFromOptions()`, `TSTrajectoryDestroy()`, `TSTrajectorySetType()`
490: @*/
491: PetscErrorCode TSTrajectoryGetType(TSTrajectory tj, TS ts, TSTrajectoryType *type)
492: {
493:   PetscFunctionBegin;
495:   if (type) *type = ((PetscObject)tj)->type_name;
496:   PetscFunctionReturn(PETSC_SUCCESS);
497: }

499: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory, TS);
500: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Singlefile(TSTrajectory, TS);
501: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory, TS);
502: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Visualization(TSTrajectory, TS);

504: /*@C
505:   TSTrajectoryRegisterAll - Registers all of the `TSTrajectory` storage schecmes in the `TS` package.

507:   Not Collective

509:   Level: developer

511: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryRegister()`
512: @*/
513: PetscErrorCode TSTrajectoryRegisterAll(void)
514: {
515:   PetscFunctionBegin;
516:   if (TSTrajectoryRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
517:   TSTrajectoryRegisterAllCalled = PETSC_TRUE;

519:   PetscCall(TSTrajectoryRegister(TSTRAJECTORYBASIC, TSTrajectoryCreate_Basic));
520:   PetscCall(TSTrajectoryRegister(TSTRAJECTORYSINGLEFILE, TSTrajectoryCreate_Singlefile));
521:   PetscCall(TSTrajectoryRegister(TSTRAJECTORYMEMORY, TSTrajectoryCreate_Memory));
522:   PetscCall(TSTrajectoryRegister(TSTRAJECTORYVISUALIZATION, TSTrajectoryCreate_Visualization));
523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

526: /*@
527:   TSTrajectoryReset - Resets a trajectory context

529:   Collective

531:   Input Parameter:
532: . tj - the `TSTrajectory` context obtained from `TSGetTrajectory()`

534:   Level: developer

536: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectorySetUp()`
537: @*/
538: PetscErrorCode TSTrajectoryReset(TSTrajectory tj)
539: {
540:   PetscFunctionBegin;
541:   if (!tj) PetscFunctionReturn(PETSC_SUCCESS);
543:   PetscTryTypeMethod(tj, reset);
544:   PetscCall(PetscFree(tj->dirfiletemplate));
545:   PetscCall(TSHistoryDestroy(&tj->tsh));
546:   PetscCall(TSHistoryCreate(PetscObjectComm((PetscObject)tj), &tj->tsh));
547:   tj->setupcalled = PETSC_FALSE;
548:   PetscFunctionReturn(PETSC_SUCCESS);
549: }

551: /*@
552:   TSTrajectoryDestroy - Destroys a trajectory context

554:   Collective

556:   Input Parameter:
557: . tj - the `TSTrajectory` context obtained from `TSTrajectoryCreate()`

559:   Level: developer

561: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectorySetUp()`
562: @*/
563: PetscErrorCode TSTrajectoryDestroy(TSTrajectory *tj)
564: {
565:   PetscFunctionBegin;
566:   if (!*tj) PetscFunctionReturn(PETSC_SUCCESS);
568:   if (--((PetscObject)*tj)->refct > 0) {
569:     *tj = NULL;
570:     PetscFunctionReturn(PETSC_SUCCESS);
571:   }

573:   PetscCall(TSTrajectoryReset(*tj));
574:   PetscCall(TSHistoryDestroy(&(*tj)->tsh));
575:   PetscCall(VecDestroyVecs((*tj)->lag.order + 1, &(*tj)->lag.W));
576:   PetscCall(PetscFree5((*tj)->lag.L, (*tj)->lag.T, (*tj)->lag.WW, (*tj)->lag.TT, (*tj)->lag.TW));
577:   PetscCall(VecDestroy(&(*tj)->U));
578:   PetscCall(VecDestroy(&(*tj)->Udot));

580:   if ((*tj)->transformdestroy) PetscCall((*(*tj)->transformdestroy)((*tj)->transformctx));
581:   PetscTryTypeMethod(*tj, destroy);
582:   if (!((*tj)->keepfiles)) {
583:     PetscMPIInt rank;
584:     MPI_Comm    comm;

586:     PetscCall(PetscObjectGetComm((PetscObject)*tj, &comm));
587:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
588:     if (rank == 0 && (*tj)->dirname) { /* we own the directory, so we run PetscRMTree on it */
589:       PetscCall(PetscRMTree((*tj)->dirname));
590:     }
591:   }
592:   PetscCall(PetscStrArrayDestroy(&(*tj)->names));
593:   PetscCall(PetscFree((*tj)->dirname));
594:   PetscCall(PetscFree((*tj)->filetemplate));
595:   PetscCall(PetscHeaderDestroy(tj));
596:   PetscFunctionReturn(PETSC_SUCCESS);
597: }

599: /*
600:   TSTrajectorySetTypeFromOptions_Private - Sets the type of `TSTrajectory` from user options.

602:   Collective

604:   Input Parameter:
605: + tj - the `TSTrajectory` context
606: - ts - the TS context

608:   Options Database Key:
609: . -ts_trajectory_type <type> - TSTRAJECTORYBASIC, TSTRAJECTORYMEMORY, TSTRAJECTORYSINGLEFILE, TSTRAJECTORYVISUALIZATION

611:   Level: developer

613: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryType`, `TSTrajectorySetFromOptions()`, `TSTrajectorySetType()`
614: */
615: static PetscErrorCode TSTrajectorySetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject, TSTrajectory tj, TS ts)
616: {
617:   PetscBool   opt;
618:   const char *defaultType;
619:   char        typeName[256];

621:   PetscFunctionBegin;
622:   if (((PetscObject)tj)->type_name) defaultType = ((PetscObject)tj)->type_name;
623:   else defaultType = TSTRAJECTORYBASIC;

625:   PetscCall(TSTrajectoryRegisterAll());
626:   PetscCall(PetscOptionsFList("-ts_trajectory_type", "TSTrajectory method", "TSTrajectorySetType", TSTrajectoryList, defaultType, typeName, 256, &opt));
627:   if (opt) {
628:     PetscCall(TSTrajectorySetType(tj, ts, typeName));
629:   } else {
630:     PetscCall(TSTrajectorySetType(tj, ts, defaultType));
631:   }
632:   PetscFunctionReturn(PETSC_SUCCESS);
633: }

635: /*@
636:   TSTrajectorySetUseHistory - Use `TSHistory` in `TSTrajectory`

638:   Collective

640:   Input Parameters:
641: + tj  - the `TSTrajectory` context
642: - flg - `PETSC_TRUE` to save, `PETSC_FALSE` to disable

644:   Options Database Key:
645: . -ts_trajectory_use_history - have it use `TSHistory`

647:   Level: advanced

649: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectorySetUp()`
650: @*/
651: PetscErrorCode TSTrajectorySetUseHistory(TSTrajectory tj, PetscBool flg)
652: {
653:   PetscFunctionBegin;
656:   tj->usehistory = flg;
657:   PetscFunctionReturn(PETSC_SUCCESS);
658: }

660: /*@
661:   TSTrajectorySetMonitor - Monitor the schedules generated by the `TSTrajectory` checkpointing controller

663:   Collective

665:   Input Parameters:
666: + tj  - the `TSTrajectory` context
667: - flg - `PETSC_TRUE` to active a monitor, `PETSC_FALSE` to disable

669:   Options Database Key:
670: . -ts_trajectory_monitor - print `TSTrajectory` information

672:   Level: developer

674: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectorySetUp()`
675: @*/
676: PetscErrorCode TSTrajectorySetMonitor(TSTrajectory tj, PetscBool flg)
677: {
678:   PetscFunctionBegin;
681:   if (flg) tj->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)tj));
682:   else tj->monitor = NULL;
683:   PetscFunctionReturn(PETSC_SUCCESS);
684: }

686: /*@
687:   TSTrajectorySetKeepFiles - Keep the files generated by the `TSTrajectory` once the program is done

689:   Collective

691:   Input Parameters:
692: + tj  - the `TSTrajectory` context
693: - flg - `PETSC_TRUE` to save, `PETSC_FALSE` to disable

695:   Options Database Key:
696: . -ts_trajectory_keep_files - have it keep the files

698:   Level: advanced

700:   Note:
701:   By default the `TSTrajectory` used for adjoint computations, `TSTRAJECTORYBASIC`, removes the files it generates at the end of the run. This causes the files to be kept.

703: .seealso: [](ch_ts), `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectorySetUp()`, `TSTrajectorySetMonitor()`
704: @*/
705: PetscErrorCode TSTrajectorySetKeepFiles(TSTrajectory tj, PetscBool flg)
706: {
707:   PetscFunctionBegin;
710:   tj->keepfiles = flg;
711:   PetscFunctionReturn(PETSC_SUCCESS);
712: }

714: /*@C
715:   TSTrajectorySetDirname - Specify the name of the directory where `TSTrajectory` disk checkpoints are stored.

717:   Collective

719:   Input Parameters:
720: + tj      - the `TSTrajectory` context
721: - dirname - the directory name

723:   Options Database Key:
724: . -ts_trajectory_dirname - set the directory name

726:   Level: developer

728:   Notes:
729:   The final location of the files is determined by dirname/filetemplate where filetemplate was provided by `TSTrajectorySetFiletemplate()`

731:   If this is not called `TSTrajectory` selects a unique new name for the directory

733: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectorySetFiletemplate()`, `TSTrajectorySetUp()`
734: @*/
735: PetscErrorCode TSTrajectorySetDirname(TSTrajectory tj, const char dirname[])
736: {
737:   PetscBool flg;

739:   PetscFunctionBegin;
741:   PetscCall(PetscStrcmp(tj->dirname, dirname, &flg));
742:   PetscCheck(flg || !tj->dirfiletemplate, PetscObjectComm((PetscObject)tj), PETSC_ERR_ARG_WRONGSTATE, "Cannot set directoryname after TSTrajectory has been setup");
743:   PetscCall(PetscFree(tj->dirname));
744:   PetscCall(PetscStrallocpy(dirname, &tj->dirname));
745:   PetscFunctionReturn(PETSC_SUCCESS);
746: }

748: /*@C
749:   TSTrajectorySetFiletemplate - Specify the name template for the files storing `TSTrajectory` checkpoints.

751:   Collective

753:   Input Parameters:
754: + tj           - the `TSTrajectory` context
755: - filetemplate - the template

757:   Options Database Key:
758: . -ts_trajectory_file_template - set the file name template

760:   Level: developer

762:   Notes:
763:   The name template should be of the form, for example filename-%06" PetscInt_FMT ".bin It should not begin with a leading /

765:   The final location of the files is determined by dirname/filetemplate where dirname was provided by `TSTrajectorySetDirname()`. The %06" PetscInt_FMT " is replaced by the
766:   timestep counter

768: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectorySetDirname()`, `TSTrajectorySetUp()`
769: @*/
770: PetscErrorCode TSTrajectorySetFiletemplate(TSTrajectory tj, const char filetemplate[])
771: {
772:   const char *ptr = NULL, *ptr2 = NULL;

774:   PetscFunctionBegin;
776:   PetscAssertPointer(filetemplate, 2);
777:   PetscCheck(!tj->dirfiletemplate, PetscObjectComm((PetscObject)tj), PETSC_ERR_ARG_WRONGSTATE, "Cannot set filetemplate after TSTrajectory has been setup");

779:   PetscCheck(filetemplate[0], PetscObjectComm((PetscObject)tj), PETSC_ERR_USER, "-ts_trajectory_file_template requires a file name template, e.g. filename-%%06" PetscInt_FMT ".bin");
780:   /* Do some cursory validation of the input. */
781:   PetscCall(PetscStrstr(filetemplate, "%", (char **)&ptr));
782:   PetscCheck(ptr, PetscObjectComm((PetscObject)tj), PETSC_ERR_USER, "-ts_trajectory_file_template requires a file name template, e.g. filename-%%06" PetscInt_FMT ".bin");
783:   for (ptr++; ptr && *ptr; ptr++) {
784:     PetscCall(PetscStrchr(PetscInt_FMT "DiouxX", *ptr, (char **)&ptr2));
785:     PetscCheck(ptr2 || (*ptr >= '0' && *ptr <= '9'), PetscObjectComm((PetscObject)tj), PETSC_ERR_USER, "Invalid file template argument to -ts_trajectory_file_template, should look like filename-%%06" PetscInt_FMT ".bin");
786:     if (ptr2) break;
787:   }
788:   PetscCall(PetscFree(tj->filetemplate));
789:   PetscCall(PetscStrallocpy(filetemplate, &tj->filetemplate));
790:   PetscFunctionReturn(PETSC_SUCCESS);
791: }

793: /*@
794:   TSTrajectorySetFromOptions - Sets various `TSTrajectory` parameters from user options.

796:   Collective

798:   Input Parameters:
799: + tj - the `TSTrajectory` context obtained from `TSGetTrajectory()`
800: - ts - the `TS` context

802:   Options Database Keys:
803: + -ts_trajectory_type <type>             - basic, memory, singlefile, visualization
804: . -ts_trajectory_keep_files <true,false> - keep the files generated by the code after the program ends. This is true by default for singlefile and visualization
805: - -ts_trajectory_monitor                 - print `TSTrajectory` information

807:   Level: developer

809:   Note:
810:   This is not normally called directly by users

812: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectorySetUp()`
813: @*/
814: PetscErrorCode TSTrajectorySetFromOptions(TSTrajectory tj, TS ts)
815: {
816:   PetscBool set, flg;
817:   char      dirname[PETSC_MAX_PATH_LEN], filetemplate[PETSC_MAX_PATH_LEN];

819:   PetscFunctionBegin;
822:   PetscObjectOptionsBegin((PetscObject)tj);
823:   PetscCall(TSTrajectorySetTypeFromOptions_Private(PetscOptionsObject, tj, ts));
824:   PetscCall(PetscOptionsBool("-ts_trajectory_use_history", "Turn on/off usage of TSHistory", NULL, tj->usehistory, &tj->usehistory, NULL));
825:   PetscCall(PetscOptionsBool("-ts_trajectory_monitor", "Print checkpointing schedules", "TSTrajectorySetMonitor", tj->monitor ? PETSC_TRUE : PETSC_FALSE, &flg, &set));
826:   if (set) PetscCall(TSTrajectorySetMonitor(tj, flg));
827:   PetscCall(PetscOptionsInt("-ts_trajectory_reconstruction_order", "Interpolation order for reconstruction", NULL, tj->lag.order, &tj->lag.order, NULL));
828:   PetscCall(PetscOptionsBool("-ts_trajectory_reconstruction_caching", "Turn on/off caching of TSTrajectoryGetVecs input", NULL, tj->lag.caching, &tj->lag.caching, NULL));
829:   PetscCall(PetscOptionsBool("-ts_trajectory_adjointmode", "Instruct the trajectory that will be used in a TSAdjointSolve()", NULL, tj->adjoint_solve_mode, &tj->adjoint_solve_mode, NULL));
830:   PetscCall(PetscOptionsBool("-ts_trajectory_solution_only", "Checkpoint solution only", "TSTrajectorySetSolutionOnly", tj->solution_only, &tj->solution_only, NULL));
831:   PetscCall(PetscOptionsBool("-ts_trajectory_keep_files", "Keep any trajectory files generated during the run", "TSTrajectorySetKeepFiles", tj->keepfiles, &flg, &set));
832:   if (set) PetscCall(TSTrajectorySetKeepFiles(tj, flg));

834:   PetscCall(PetscOptionsString("-ts_trajectory_dirname", "Directory name for TSTrajectory file", "TSTrajectorySetDirname", NULL, dirname, sizeof(dirname) - 14, &set));
835:   if (set) PetscCall(TSTrajectorySetDirname(tj, dirname));

837:   PetscCall(PetscOptionsString("-ts_trajectory_file_template", "Template for TSTrajectory file name, use filename-%06" PetscInt_FMT ".bin", "TSTrajectorySetFiletemplate", NULL, filetemplate, sizeof(filetemplate), &set));
838:   if (set) PetscCall(TSTrajectorySetFiletemplate(tj, filetemplate));

840:   /* Handle specific TSTrajectory options */
841:   PetscTryTypeMethod(tj, setfromoptions, PetscOptionsObject);
842:   PetscOptionsEnd();
843:   PetscFunctionReturn(PETSC_SUCCESS);
844: }

846: /*@
847:   TSTrajectorySetUp - Sets up the internal data structures, e.g. stacks, for the later use
848:   of a `TS` `TSTrajectory`.

850:   Collective

852:   Input Parameters:
853: + tj - the `TSTrajectory` context
854: - ts - the TS context obtained from `TSCreate()`

856:   Level: developer

858: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`
859: @*/
860: PetscErrorCode TSTrajectorySetUp(TSTrajectory tj, TS ts)
861: {
862:   size_t s1, s2;

864:   PetscFunctionBegin;
865:   if (!tj) PetscFunctionReturn(PETSC_SUCCESS);
868:   if (tj->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);

870:   PetscCall(PetscLogEventBegin(TSTrajectory_SetUp, tj, ts, 0, 0));
871:   if (!((PetscObject)tj)->type_name) PetscCall(TSTrajectorySetType(tj, ts, TSTRAJECTORYBASIC));
872:   PetscTryTypeMethod(tj, setup, ts);

874:   tj->setupcalled = PETSC_TRUE;

876:   /* Set the counters to zero */
877:   tj->recomps    = 0;
878:   tj->diskreads  = 0;
879:   tj->diskwrites = 0;
880:   PetscCall(PetscStrlen(tj->dirname, &s1));
881:   PetscCall(PetscStrlen(tj->filetemplate, &s2));
882:   PetscCall(PetscFree(tj->dirfiletemplate));
883:   PetscCall(PetscMalloc((s1 + s2 + 10) * sizeof(char), &tj->dirfiletemplate));
884:   PetscCall(PetscSNPrintf(tj->dirfiletemplate, s1 + s2 + 10, "%s/%s", tj->dirname, tj->filetemplate));
885:   PetscCall(PetscLogEventEnd(TSTrajectory_SetUp, tj, ts, 0, 0));
886:   PetscFunctionReturn(PETSC_SUCCESS);
887: }

889: /*@
890:   TSTrajectorySetSolutionOnly - Tells the trajectory to store just the solution, and not any intermediate stage information

892:   Collective

894:   Input Parameters:
895: + tj            - the `TSTrajectory` context obtained with `TSGetTrajectory()`
896: - solution_only - the boolean flag

898:   Level: developer

900: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectoryGetSolutionOnly()`
901: @*/
902: PetscErrorCode TSTrajectorySetSolutionOnly(TSTrajectory tj, PetscBool solution_only)
903: {
904:   PetscFunctionBegin;
907:   tj->solution_only = solution_only;
908:   PetscFunctionReturn(PETSC_SUCCESS);
909: }

911: /*@
912:   TSTrajectoryGetSolutionOnly - Gets the value set with `TSTrajectorySetSolutionOnly()`.

914:   Logically Collective

916:   Input Parameter:
917: . tj - the `TSTrajectory` context

919:   Output Parameter:
920: . solution_only - the boolean flag

922:   Level: developer

924: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectorySetSolutionOnly()`
925: @*/
926: PetscErrorCode TSTrajectoryGetSolutionOnly(TSTrajectory tj, PetscBool *solution_only)
927: {
928:   PetscFunctionBegin;
930:   PetscAssertPointer(solution_only, 2);
931:   *solution_only = tj->solution_only;
932:   PetscFunctionReturn(PETSC_SUCCESS);
933: }

935: /*@
936:   TSTrajectoryGetUpdatedHistoryVecs - Get updated state and time-derivative history vectors.

938:   Collective

940:   Input Parameters:
941: + tj   - the `TSTrajectory` context
942: . ts   - the `TS` solver context
943: - time - the requested time

945:   Output Parameters:
946: + U    - state vector at given time (can be interpolated)
947: - Udot - time-derivative vector at given time (can be interpolated)

949:   Level: developer

951:   Notes:
952:   The vectors are interpolated if time does not match any time step stored in the `TSTrajectory()`. Pass NULL to not request a vector.

954:   This function differs from `TSTrajectoryGetVecs()` since the vectors obtained cannot be modified, and they need to be returned by
955:   calling `TSTrajectoryRestoreUpdatedHistoryVecs()`.

957: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectoryRestoreUpdatedHistoryVecs()`, `TSTrajectoryGetVecs()`
958: @*/
959: PetscErrorCode TSTrajectoryGetUpdatedHistoryVecs(TSTrajectory tj, TS ts, PetscReal time, Vec *U, Vec *Udot)
960: {
961:   PetscFunctionBegin;
965:   if (U) PetscAssertPointer(U, 4);
966:   if (Udot) PetscAssertPointer(Udot, 5);
967:   if (U && !tj->U) {
968:     DM dm;

970:     PetscCall(TSGetDM(ts, &dm));
971:     PetscCall(DMCreateGlobalVector(dm, &tj->U));
972:   }
973:   if (Udot && !tj->Udot) {
974:     DM dm;

976:     PetscCall(TSGetDM(ts, &dm));
977:     PetscCall(DMCreateGlobalVector(dm, &tj->Udot));
978:   }
979:   PetscCall(TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &time, U ? tj->U : NULL, Udot ? tj->Udot : NULL));
980:   if (U) {
981:     PetscCall(VecLockReadPush(tj->U));
982:     *U = tj->U;
983:   }
984:   if (Udot) {
985:     PetscCall(VecLockReadPush(tj->Udot));
986:     *Udot = tj->Udot;
987:   }
988:   PetscFunctionReturn(PETSC_SUCCESS);
989: }

991: /*@
992:   TSTrajectoryRestoreUpdatedHistoryVecs - Restores updated state and time-derivative history vectors obtained with `TSTrajectoryGetUpdatedHistoryVecs()`.

994:   Collective

996:   Input Parameters:
997: + tj   - the `TSTrajectory` context
998: . U    - state vector at given time (can be interpolated)
999: - Udot - time-derivative vector at given time (can be interpolated)

1001:   Level: developer

1003: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryGetUpdatedHistoryVecs()`
1004: @*/
1005: PetscErrorCode TSTrajectoryRestoreUpdatedHistoryVecs(TSTrajectory tj, Vec *U, Vec *Udot)
1006: {
1007:   PetscFunctionBegin;
1011:   PetscCheck(!U || *U == tj->U, PetscObjectComm((PetscObject)*U), PETSC_ERR_USER, "U was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1012:   PetscCheck(!Udot || *Udot == tj->Udot, PetscObjectComm((PetscObject)*Udot), PETSC_ERR_USER, "Udot was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1013:   if (U) {
1014:     PetscCall(VecLockReadPop(tj->U));
1015:     *U = NULL;
1016:   }
1017:   if (Udot) {
1018:     PetscCall(VecLockReadPop(tj->Udot));
1019:     *Udot = NULL;
1020:   }
1021:   PetscFunctionReturn(PETSC_SUCCESS);
1022: }