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, No Fortran Support
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: /*@
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: Options Database Key:
266: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
268: Level: intermediate
270: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryView`, `PetscObjectViewFromOptions()`, `TSTrajectoryCreate()`
271: @*/
272: PetscErrorCode TSTrajectoryViewFromOptions(TSTrajectory A, PetscObject obj, const char name[])
273: {
274: PetscFunctionBegin;
276: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*@
281: TSTrajectoryView - Prints information about the trajectory object
283: Collective
285: Input Parameters:
286: + tj - the `TSTrajectory` context obtained from `TSTrajectoryCreate()`
287: - viewer - visualization context
289: Options Database Key:
290: . -ts_trajectory_view - calls `TSTrajectoryView()` at end of `TSAdjointStep()`
292: Level: developer
294: Notes:
295: The available visualization contexts include
296: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
297: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
298: output where only the first processor opens
299: the file. All other processors send their
300: data to the first processor to print.
302: The user can open an alternative visualization context with
303: `PetscViewerASCIIOpen()` - output to a specified file.
305: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `PetscViewer`, `PetscViewerASCIIOpen()`
306: @*/
307: PetscErrorCode TSTrajectoryView(TSTrajectory tj, PetscViewer viewer)
308: {
309: PetscBool isascii;
311: PetscFunctionBegin;
313: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tj), &viewer));
315: PetscCheckSameComm(tj, 1, viewer, 2);
317: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
318: if (isascii) {
319: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tj, viewer));
320: PetscCall(PetscViewerASCIIPrintf(viewer, " total number of recomputations for adjoint calculation = %" PetscInt_FMT "\n", tj->recomps));
321: PetscCall(PetscViewerASCIIPrintf(viewer, " disk checkpoint reads = %" PetscInt_FMT "\n", tj->diskreads));
322: PetscCall(PetscViewerASCIIPrintf(viewer, " disk checkpoint writes = %" PetscInt_FMT "\n", tj->diskwrites));
323: PetscCall(PetscViewerASCIIPushTab(viewer));
324: PetscTryTypeMethod(tj, view, viewer);
325: PetscCall(PetscViewerASCIIPopTab(viewer));
326: }
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /*@C
331: TSTrajectorySetVariableNames - Sets the name of each component in the solution vector so that it may be saved with the trajectory
333: Collective
335: Input Parameters:
336: + ctx - the trajectory context
337: - names - the names of the components, final string must be `NULL`
339: Level: intermediate
341: Fortran Notes:
342: Fortran interface is not possible because of the string array argument
344: .seealso: [](ch_ts), `TSTrajectory`, `TSGetTrajectory()`
345: @*/
346: PetscErrorCode TSTrajectorySetVariableNames(TSTrajectory ctx, const char *const *names)
347: {
348: PetscFunctionBegin;
350: PetscAssertPointer(names, 2);
351: PetscCall(PetscStrArrayDestroy(&ctx->names));
352: PetscCall(PetscStrArrayallocpy(names, &ctx->names));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@C
357: TSTrajectorySetTransform - Solution vector will be transformed by provided function before being saved to disk
359: Collective
361: Input Parameters:
362: + tj - the `TSTrajectory` context
363: . transform - the transform function
364: . destroy - function to destroy the optional context
365: - tctx - optional context used by transform function
367: Level: intermediate
369: .seealso: [](ch_ts), `TSTrajectorySetVariableNames()`, `TSTrajectory`, `TSMonitorLGSetTransform()`
370: @*/
371: PetscErrorCode TSTrajectorySetTransform(TSTrajectory tj, PetscErrorCode (*transform)(void *, Vec, Vec *), PetscCtxDestroyFn *destroy, void *tctx)
372: {
373: PetscFunctionBegin;
375: tj->transform = transform;
376: tj->transformdestroy = destroy;
377: tj->transformctx = tctx;
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*@
382: TSTrajectoryCreate - This function creates an empty trajectory object used to store the time dependent solution of an ODE/DAE
384: Collective
386: Input Parameter:
387: . comm - the communicator
389: Output Parameter:
390: . tj - the trajectory object
392: Level: developer
394: Notes:
395: Usually one does not call this routine, it is called automatically when one calls `TSSetSaveTrajectory()`.
397: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSTrajectorySetUp()`, `TSTrajectoryDestroy()`, `TSTrajectorySetType()`, `TSTrajectorySetVariableNames()`, `TSGetTrajectory()`, `TSTrajectorySetKeepFiles()`
398: @*/
399: PetscErrorCode TSTrajectoryCreate(MPI_Comm comm, TSTrajectory *tj)
400: {
401: TSTrajectory t;
403: PetscFunctionBegin;
404: PetscAssertPointer(tj, 2);
405: PetscCall(TSInitializePackage());
407: PetscCall(PetscHeaderCreate(t, TSTRAJECTORY_CLASSID, "TSTrajectory", "Time stepping", "TS", comm, TSTrajectoryDestroy, TSTrajectoryView));
408: t->setupcalled = PETSC_FALSE;
409: PetscCall(TSHistoryCreate(comm, &t->tsh));
410: t->lag.order = 1;
411: t->lag.L = NULL;
412: t->lag.T = NULL;
413: t->lag.W = NULL;
414: t->lag.WW = NULL;
415: t->lag.TW = NULL;
416: t->lag.TT = NULL;
417: t->lag.caching = PETSC_TRUE;
418: t->lag.Ucached.id = 0;
419: t->lag.Ucached.state = -1;
420: t->lag.Ucached.time = PETSC_MIN_REAL;
421: t->lag.Ucached.step = PETSC_INT_MAX;
422: t->lag.Udotcached.id = 0;
423: t->lag.Udotcached.state = -1;
424: t->lag.Udotcached.time = PETSC_MIN_REAL;
425: t->lag.Udotcached.step = PETSC_INT_MAX;
426: t->adjoint_solve_mode = PETSC_TRUE;
427: t->solution_only = PETSC_FALSE;
428: t->keepfiles = PETSC_FALSE;
429: t->usehistory = PETSC_TRUE;
430: *tj = t;
431: PetscCall(TSTrajectorySetFiletemplate(t, "TS-%06" PetscInt_FMT ".bin"));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: /*@
436: TSTrajectorySetType - Sets the storage method to be used as in a trajectory
438: Collective
440: Input Parameters:
441: + tj - the `TSTrajectory` context
442: . ts - the `TS` context
443: - type - a known method
445: Options Database Key:
446: . -ts_trajectory_type type - Sets the method; use -help for a list of available methods (for instance, basic)
448: Level: developer
450: Developer Notes:
451: Why does this option require access to the `TS`
453: .seealso: [](ch_ts), `TSTrajectory`, `TS`, `TSTrajectoryCreate()`, `TSTrajectorySetFromOptions()`, `TSTrajectoryDestroy()`, `TSTrajectoryGetType()`
454: @*/
455: PetscErrorCode TSTrajectorySetType(TSTrajectory tj, TS ts, TSTrajectoryType type)
456: {
457: PetscErrorCode (*r)(TSTrajectory, TS);
458: PetscBool match;
460: PetscFunctionBegin;
462: PetscCall(PetscObjectTypeCompare((PetscObject)tj, type, &match));
463: if (match) PetscFunctionReturn(PETSC_SUCCESS);
465: PetscCall(PetscFunctionListFind(TSTrajectoryList, type, &r));
466: PetscCheck(r, PetscObjectComm((PetscObject)tj), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSTrajectory type: %s", type);
467: PetscTryTypeMethod(tj, destroy);
468: tj->ops->destroy = NULL;
469: PetscCall(PetscMemzero(tj->ops, sizeof(*tj->ops)));
471: PetscCall(PetscObjectChangeTypeName((PetscObject)tj, type));
472: PetscCall((*r)(tj, ts));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*@
477: TSTrajectoryGetType - Gets the trajectory type
479: Collective
481: Input Parameters:
482: + tj - the `TSTrajectory` context
483: - ts - the `TS` context
485: Output Parameter:
486: . type - a known method
488: Level: developer
490: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectorySetFromOptions()`, `TSTrajectoryDestroy()`, `TSTrajectorySetType()`
491: @*/
492: PetscErrorCode TSTrajectoryGetType(TSTrajectory tj, TS ts, TSTrajectoryType *type)
493: {
494: PetscFunctionBegin;
496: if (type) *type = ((PetscObject)tj)->type_name;
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory, TS);
501: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Singlefile(TSTrajectory, TS);
502: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Memory(TSTrajectory, TS);
503: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Visualization(TSTrajectory, TS);
505: /*@C
506: TSTrajectoryRegisterAll - Registers all of the `TSTrajectory` storage schecmes in the `TS` package.
508: Not Collective
510: Level: developer
512: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryRegister()`
513: @*/
514: PetscErrorCode TSTrajectoryRegisterAll(void)
515: {
516: PetscFunctionBegin;
517: if (TSTrajectoryRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
518: TSTrajectoryRegisterAllCalled = PETSC_TRUE;
520: PetscCall(TSTrajectoryRegister(TSTRAJECTORYBASIC, TSTrajectoryCreate_Basic));
521: PetscCall(TSTrajectoryRegister(TSTRAJECTORYSINGLEFILE, TSTrajectoryCreate_Singlefile));
522: PetscCall(TSTrajectoryRegister(TSTRAJECTORYMEMORY, TSTrajectoryCreate_Memory));
523: PetscCall(TSTrajectoryRegister(TSTRAJECTORYVISUALIZATION, TSTrajectoryCreate_Visualization));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@
528: TSTrajectoryReset - Resets a trajectory context
530: Collective
532: Input Parameter:
533: . tj - the `TSTrajectory` context obtained from `TSGetTrajectory()`
535: Level: developer
537: .seealso: [](ch_ts), `TS`, `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectorySetUp()`
538: @*/
539: PetscErrorCode TSTrajectoryReset(TSTrajectory tj)
540: {
541: PetscFunctionBegin;
542: if (!tj) PetscFunctionReturn(PETSC_SUCCESS);
544: PetscTryTypeMethod(tj, reset);
545: PetscCall(PetscFree(tj->dirfiletemplate));
546: PetscCall(TSHistoryDestroy(&tj->tsh));
547: PetscCall(TSHistoryCreate(PetscObjectComm((PetscObject)tj), &tj->tsh));
548: tj->setupcalled = PETSC_FALSE;
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: /*@
553: TSTrajectoryDestroy - Destroys a trajectory context
555: Collective
557: Input Parameter:
558: . tj - the `TSTrajectory` context obtained from `TSTrajectoryCreate()`
560: Level: developer
562: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectorySetUp()`
563: @*/
564: PetscErrorCode TSTrajectoryDestroy(TSTrajectory *tj)
565: {
566: PetscFunctionBegin;
567: if (!*tj) PetscFunctionReturn(PETSC_SUCCESS);
569: if (--((PetscObject)*tj)->refct > 0) {
570: *tj = NULL;
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: PetscCall(TSTrajectoryReset(*tj));
575: PetscCall(TSHistoryDestroy(&(*tj)->tsh));
576: PetscCall(VecDestroyVecs((*tj)->lag.order + 1, &(*tj)->lag.W));
577: PetscCall(PetscFree5((*tj)->lag.L, (*tj)->lag.T, (*tj)->lag.WW, (*tj)->lag.TT, (*tj)->lag.TW));
578: PetscCall(VecDestroy(&(*tj)->U));
579: PetscCall(VecDestroy(&(*tj)->Udot));
581: if ((*tj)->transformdestroy) PetscCall((*(*tj)->transformdestroy)(&(*tj)->transformctx));
582: PetscTryTypeMethod(*tj, destroy);
583: if (!((*tj)->keepfiles)) {
584: PetscMPIInt rank;
585: MPI_Comm comm;
587: PetscCall(PetscObjectGetComm((PetscObject)*tj, &comm));
588: PetscCallMPI(MPI_Comm_rank(comm, &rank));
589: if (rank == 0 && (*tj)->dirname) { /* we own the directory, so we run PetscRMTree on it */
590: PetscCall(PetscRMTree((*tj)->dirname));
591: }
592: }
593: PetscCall(PetscStrArrayDestroy(&(*tj)->names));
594: PetscCall(PetscFree((*tj)->dirname));
595: PetscCall(PetscFree((*tj)->filetemplate));
596: PetscCall(PetscHeaderDestroy(tj));
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /*
601: TSTrajectorySetTypeFromOptions_Private - Sets the type of `TSTrajectory` from user options.
603: Collective
605: Input Parameter:
606: + tj - the `TSTrajectory` context
607: - ts - the TS context
609: Options Database Key:
610: . -ts_trajectory_type (basic|memory|singlefile|visualization) - the method of storing the trajectory
612: Level: developer
614: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryType`, `TSTrajectorySetFromOptions()`, `TSTrajectorySetType()`
615: */
616: static PetscErrorCode TSTrajectorySetTypeFromOptions_Private(PetscOptionItems PetscOptionsObject, TSTrajectory tj, TS ts)
617: {
618: PetscBool opt;
619: const char *defaultType;
620: char typeName[256];
622: PetscFunctionBegin;
623: if (((PetscObject)tj)->type_name) defaultType = ((PetscObject)tj)->type_name;
624: else defaultType = TSTRAJECTORYBASIC;
626: PetscCall(TSTrajectoryRegisterAll());
627: PetscCall(PetscOptionsFList("-ts_trajectory_type", "TSTrajectory method", "TSTrajectorySetType", TSTrajectoryList, defaultType, typeName, 256, &opt));
628: if (opt) {
629: PetscCall(TSTrajectorySetType(tj, ts, typeName));
630: } else {
631: PetscCall(TSTrajectorySetType(tj, ts, defaultType));
632: }
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*@
637: TSTrajectorySetUseHistory - Use `TSHistory` in `TSTrajectory`
639: Collective
641: Input Parameters:
642: + tj - the `TSTrajectory` context
643: - flg - `PETSC_TRUE` to save, `PETSC_FALSE` to disable
645: Options Database Key:
646: . -ts_trajectory_use_history - have it use `TSHistory`
648: Level: advanced
650: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectorySetUp()`
651: @*/
652: PetscErrorCode TSTrajectorySetUseHistory(TSTrajectory tj, PetscBool flg)
653: {
654: PetscFunctionBegin;
657: tj->usehistory = flg;
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: /*@
662: TSTrajectorySetMonitor - Monitor the schedules generated by the `TSTrajectory` checkpointing controller
664: Collective
666: Input Parameters:
667: + tj - the `TSTrajectory` context
668: - flg - `PETSC_TRUE` to active a monitor, `PETSC_FALSE` to disable
670: Options Database Key:
671: . -ts_trajectory_monitor - print `TSTrajectory` information
673: Level: developer
675: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectorySetUp()`
676: @*/
677: PetscErrorCode TSTrajectorySetMonitor(TSTrajectory tj, PetscBool flg)
678: {
679: PetscFunctionBegin;
682: if (flg) tj->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)tj));
683: else tj->monitor = NULL;
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: /*@
688: TSTrajectorySetKeepFiles - Keep the files generated by the `TSTrajectory` once the program is done
690: Collective
692: Input Parameters:
693: + tj - the `TSTrajectory` context
694: - flg - `PETSC_TRUE` to save, `PETSC_FALSE` to disable
696: Options Database Key:
697: . -ts_trajectory_keep_files - have it keep the files
699: Level: advanced
701: Note:
702: 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.
704: .seealso: [](ch_ts), `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectorySetUp()`, `TSTrajectorySetMonitor()`
705: @*/
706: PetscErrorCode TSTrajectorySetKeepFiles(TSTrajectory tj, PetscBool flg)
707: {
708: PetscFunctionBegin;
711: tj->keepfiles = flg;
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: /*@
716: TSTrajectorySetDirname - Specify the name of the directory where `TSTrajectory` disk checkpoints are stored.
718: Collective
720: Input Parameters:
721: + tj - the `TSTrajectory` context
722: - dirname - the directory name
724: Options Database Key:
725: . -ts_trajectory_dirname - set the directory name
727: Level: developer
729: Notes:
730: The final location of the files is determined by dirname/filetemplate where filetemplate was provided by `TSTrajectorySetFiletemplate()`
732: If this is not called `TSTrajectory` selects a unique new name for the directory
734: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectorySetFiletemplate()`, `TSTrajectorySetUp()`
735: @*/
736: PetscErrorCode TSTrajectorySetDirname(TSTrajectory tj, const char dirname[])
737: {
738: PetscBool flg;
740: PetscFunctionBegin;
742: PetscCall(PetscStrcmp(tj->dirname, dirname, &flg));
743: PetscCheck(flg || !tj->dirfiletemplate, PetscObjectComm((PetscObject)tj), PETSC_ERR_ARG_WRONGSTATE, "Cannot set directoryname after TSTrajectory has been setup");
744: PetscCall(PetscFree(tj->dirname));
745: PetscCall(PetscStrallocpy(dirname, &tj->dirname));
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: /*@
750: TSTrajectorySetFiletemplate - Specify the name template for the files storing `TSTrajectory` checkpoints.
752: Collective
754: Input Parameters:
755: + tj - the `TSTrajectory` context
756: - filetemplate - the template
758: Options Database Key:
759: . -ts_trajectory_file_template - set the file name template
761: Level: developer
763: Notes:
764: The name template should be of the form, for example filename-%06" PetscInt_FMT ".bin It should not begin with a leading /
766: 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
767: timestep counter
769: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectorySetDirname()`, `TSTrajectorySetUp()`
770: @*/
771: PetscErrorCode TSTrajectorySetFiletemplate(TSTrajectory tj, const char filetemplate[])
772: {
773: const char *ptr = NULL, *ptr2 = NULL;
775: PetscFunctionBegin;
777: PetscAssertPointer(filetemplate, 2);
778: PetscCheck(!tj->dirfiletemplate, PetscObjectComm((PetscObject)tj), PETSC_ERR_ARG_WRONGSTATE, "Cannot set filetemplate after TSTrajectory has been setup");
780: 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");
781: /* Do some cursory validation of the input. */
782: PetscCall(PetscStrstr(filetemplate, "%", (char **)&ptr));
783: PetscCheck(ptr, PetscObjectComm((PetscObject)tj), PETSC_ERR_USER, "-ts_trajectory_file_template requires a file name template, e.g. filename-%%06" PetscInt_FMT ".bin");
784: for (ptr++; ptr && *ptr; ptr++) {
785: PetscCall(PetscStrchr(PetscInt_FMT "DiouxX", *ptr, (char **)&ptr2));
786: 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");
787: if (ptr2) break;
788: }
789: PetscCall(PetscFree(tj->filetemplate));
790: PetscCall(PetscStrallocpy(filetemplate, &tj->filetemplate));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: /*@
795: TSTrajectorySetFromOptions - Sets various `TSTrajectory` parameters from user options.
797: Collective
799: Input Parameters:
800: + tj - the `TSTrajectory` context obtained from `TSGetTrajectory()`
801: - ts - the `TS` context
803: Options Database Keys:
804: + -ts_trajectory_type type - basic, memory, singlefile, visualization
805: . -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
806: - -ts_trajectory_monitor - print `TSTrajectory` information
808: Level: developer
810: Note:
811: This is not normally called directly by users
813: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectorySetUp()`
814: @*/
815: PetscErrorCode TSTrajectorySetFromOptions(TSTrajectory tj, TS ts)
816: {
817: PetscBool set, flg;
818: char dirname[PETSC_MAX_PATH_LEN], filetemplate[PETSC_MAX_PATH_LEN];
820: PetscFunctionBegin;
823: PetscObjectOptionsBegin((PetscObject)tj);
824: PetscCall(TSTrajectorySetTypeFromOptions_Private(PetscOptionsObject, tj, ts));
825: PetscCall(PetscOptionsBool("-ts_trajectory_use_history", "Turn on/off usage of TSHistory", NULL, tj->usehistory, &tj->usehistory, NULL));
826: PetscCall(PetscOptionsBool("-ts_trajectory_monitor", "Print checkpointing schedules", "TSTrajectorySetMonitor", tj->monitor ? PETSC_TRUE : PETSC_FALSE, &flg, &set));
827: if (set) PetscCall(TSTrajectorySetMonitor(tj, flg));
828: PetscCall(PetscOptionsInt("-ts_trajectory_reconstruction_order", "Interpolation order for reconstruction", NULL, tj->lag.order, &tj->lag.order, NULL));
829: PetscCall(PetscOptionsBool("-ts_trajectory_reconstruction_caching", "Turn on/off caching of TSTrajectoryGetVecs input", NULL, tj->lag.caching, &tj->lag.caching, NULL));
830: 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));
831: PetscCall(PetscOptionsBool("-ts_trajectory_solution_only", "Checkpoint solution only", "TSTrajectorySetSolutionOnly", tj->solution_only, &tj->solution_only, NULL));
832: PetscCall(PetscOptionsBool("-ts_trajectory_keep_files", "Keep any trajectory files generated during the run", "TSTrajectorySetKeepFiles", tj->keepfiles, &flg, &set));
833: if (set) PetscCall(TSTrajectorySetKeepFiles(tj, flg));
835: PetscCall(PetscOptionsString("-ts_trajectory_dirname", "Directory name for TSTrajectory file", "TSTrajectorySetDirname", NULL, dirname, sizeof(dirname) - 14, &set));
836: if (set) PetscCall(TSTrajectorySetDirname(tj, dirname));
838: PetscCall(PetscOptionsString("-ts_trajectory_file_template", "Template for TSTrajectory file name, use filename-%06" PetscInt_FMT ".bin", "TSTrajectorySetFiletemplate", NULL, filetemplate, sizeof(filetemplate), &set));
839: if (set) PetscCall(TSTrajectorySetFiletemplate(tj, filetemplate));
841: /* Handle specific TSTrajectory options */
842: PetscTryTypeMethod(tj, setfromoptions, PetscOptionsObject);
843: PetscOptionsEnd();
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: /*@
848: TSTrajectorySetUp - Sets up the internal data structures, e.g. stacks, for the later use
849: of a `TS` `TSTrajectory`.
851: Collective
853: Input Parameters:
854: + tj - the `TSTrajectory` context
855: - ts - the TS context obtained from `TSCreate()`
857: Level: developer
859: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`
860: @*/
861: PetscErrorCode TSTrajectorySetUp(TSTrajectory tj, TS ts)
862: {
863: size_t s1, s2;
865: PetscFunctionBegin;
866: if (!tj) PetscFunctionReturn(PETSC_SUCCESS);
869: if (tj->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
871: PetscCall(PetscLogEventBegin(TSTrajectory_SetUp, tj, ts, 0, 0));
872: if (!((PetscObject)tj)->type_name) PetscCall(TSTrajectorySetType(tj, ts, TSTRAJECTORYBASIC));
873: PetscTryTypeMethod(tj, setup, ts);
875: tj->setupcalled = PETSC_TRUE;
877: /* Set the counters to zero */
878: tj->recomps = 0;
879: tj->diskreads = 0;
880: tj->diskwrites = 0;
881: PetscCall(PetscStrlen(tj->dirname, &s1));
882: PetscCall(PetscStrlen(tj->filetemplate, &s2));
883: PetscCall(PetscFree(tj->dirfiletemplate));
884: PetscCall(PetscMalloc((s1 + s2 + 10) * sizeof(char), &tj->dirfiletemplate));
885: PetscCall(PetscSNPrintf(tj->dirfiletemplate, s1 + s2 + 10, "%s/%s", tj->dirname, tj->filetemplate));
886: PetscCall(PetscLogEventEnd(TSTrajectory_SetUp, tj, ts, 0, 0));
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: /*@
891: TSTrajectorySetSolutionOnly - Tells the trajectory to store just the solution, and not any intermediate stage information
893: Collective
895: Input Parameters:
896: + tj - the `TSTrajectory` context obtained with `TSGetTrajectory()`
897: - solution_only - the boolean flag
899: Level: developer
901: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectoryGetSolutionOnly()`
902: @*/
903: PetscErrorCode TSTrajectorySetSolutionOnly(TSTrajectory tj, PetscBool solution_only)
904: {
905: PetscFunctionBegin;
908: tj->solution_only = solution_only;
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: /*@
913: TSTrajectoryGetSolutionOnly - Gets the value set with `TSTrajectorySetSolutionOnly()`.
915: Logically Collective
917: Input Parameter:
918: . tj - the `TSTrajectory` context
920: Output Parameter:
921: . solution_only - the boolean flag
923: Level: developer
925: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectorySetSolutionOnly()`
926: @*/
927: PetscErrorCode TSTrajectoryGetSolutionOnly(TSTrajectory tj, PetscBool *solution_only)
928: {
929: PetscFunctionBegin;
931: PetscAssertPointer(solution_only, 2);
932: *solution_only = tj->solution_only;
933: PetscFunctionReturn(PETSC_SUCCESS);
934: }
936: /*@
937: TSTrajectoryGetUpdatedHistoryVecs - Get updated state and time-derivative history vectors.
939: Collective
941: Input Parameters:
942: + tj - the `TSTrajectory` context
943: . ts - the `TS` solver context
944: - time - the requested time
946: Output Parameters:
947: + U - state vector at given time (can be interpolated)
948: - Udot - time-derivative vector at given time (can be interpolated)
950: Level: developer
952: Notes:
953: The vectors are interpolated if time does not match any time step stored in the `TSTrajectory()`. Pass NULL to not request a vector.
955: This function differs from `TSTrajectoryGetVecs()` since the vectors obtained cannot be modified, and they need to be returned by
956: calling `TSTrajectoryRestoreUpdatedHistoryVecs()`.
958: .seealso: [](ch_ts), `TSTrajectory`, `TSSetSaveTrajectory()`, `TSTrajectoryCreate()`, `TSTrajectoryDestroy()`, `TSTrajectoryRestoreUpdatedHistoryVecs()`, `TSTrajectoryGetVecs()`
959: @*/
960: PetscErrorCode TSTrajectoryGetUpdatedHistoryVecs(TSTrajectory tj, TS ts, PetscReal time, Vec *U, Vec *Udot)
961: {
962: PetscFunctionBegin;
966: if (U) PetscAssertPointer(U, 4);
967: if (Udot) PetscAssertPointer(Udot, 5);
968: if (U && !tj->U) {
969: DM dm;
971: PetscCall(TSGetDM(ts, &dm));
972: PetscCall(DMCreateGlobalVector(dm, &tj->U));
973: }
974: if (Udot && !tj->Udot) {
975: DM dm;
977: PetscCall(TSGetDM(ts, &dm));
978: PetscCall(DMCreateGlobalVector(dm, &tj->Udot));
979: }
980: PetscCall(TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &time, U ? tj->U : NULL, Udot ? tj->Udot : NULL));
981: if (U) {
982: PetscCall(VecLockReadPush(tj->U));
983: *U = tj->U;
984: }
985: if (Udot) {
986: PetscCall(VecLockReadPush(tj->Udot));
987: *Udot = tj->Udot;
988: }
989: PetscFunctionReturn(PETSC_SUCCESS);
990: }
992: /*@
993: TSTrajectoryRestoreUpdatedHistoryVecs - Restores updated state and time-derivative history vectors obtained with `TSTrajectoryGetUpdatedHistoryVecs()`.
995: Collective
997: Input Parameters:
998: + tj - the `TSTrajectory` context
999: . U - state vector at given time (can be interpolated)
1000: - Udot - time-derivative vector at given time (can be interpolated)
1002: Level: developer
1004: .seealso: [](ch_ts), `TSTrajectory`, `TSTrajectoryGetUpdatedHistoryVecs()`
1005: @*/
1006: PetscErrorCode TSTrajectoryRestoreUpdatedHistoryVecs(TSTrajectory tj, Vec *U, Vec *Udot)
1007: {
1008: PetscFunctionBegin;
1012: PetscCheck(!U || *U == tj->U, PetscObjectComm((PetscObject)*U), PETSC_ERR_USER, "U was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1013: PetscCheck(!Udot || *Udot == tj->Udot, PetscObjectComm((PetscObject)*Udot), PETSC_ERR_USER, "Udot was not obtained from TSTrajectoryGetUpdatedHistoryVecs()");
1014: if (U) {
1015: PetscCall(VecLockReadPop(tj->U));
1016: *U = NULL;
1017: }
1018: if (Udot) {
1019: PetscCall(VecLockReadPop(tj->Udot));
1020: *Udot = NULL;
1021: }
1022: PetscFunctionReturn(PETSC_SUCCESS);
1023: }