Actual source code: trajbasic.c
1: #include <petsc/private/tsimpl.h>
3: typedef struct {
4: PetscViewer viewer;
5: } TSTrajectory_Basic;
6: /*
7: For n-th time step, TSTrajectorySet_Basic always saves the solution X(t_n) and the current time t_n,
8: and optionally saves the stage values Y[] between t_{n-1} and t_n, the previous time t_{n-1}, and
9: forward stage sensitivities S[] = dY[]/dp.
10: */
11: static PetscErrorCode TSTrajectorySet_Basic(TSTrajectory tj, TS ts, PetscInt stepnum, PetscReal time, Vec X)
12: {
13: TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic *)tj->data;
14: char filename[PETSC_MAX_PATH_LEN];
15: PetscInt ns, i;
17: PetscFunctionBegin;
18: PetscCall(PetscSNPrintf(filename, sizeof(filename), tj->dirfiletemplate, stepnum));
19: PetscCall(PetscViewerFileSetName(tjbasic->viewer, filename)); /* this triggers PetscViewer to be set up again */
20: PetscCall(PetscViewerSetUp(tjbasic->viewer));
21: PetscCall(VecView(X, tjbasic->viewer));
22: PetscCall(PetscViewerBinaryWrite(tjbasic->viewer, &time, 1, PETSC_REAL));
23: if (stepnum && !tj->solution_only) {
24: Vec *Y;
25: PetscReal tprev;
26: PetscCall(TSGetStages(ts, &ns, &Y));
27: for (i = 0; i < ns; i++) {
28: /* 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. */
29: if (ts->stifflyaccurate && i == ns - 1) continue;
30: PetscCall(VecView(Y[i], tjbasic->viewer));
31: }
32: PetscCall(TSGetPrevTime(ts, &tprev));
33: PetscCall(PetscViewerBinaryWrite(tjbasic->viewer, &tprev, 1, PETSC_REAL));
34: }
35: /* Tangent linear sensitivities needed by second-order adjoint */
36: if (ts->forward_solve) {
37: Mat A, *S;
39: PetscCall(TSForwardGetSensitivities(ts, NULL, &A));
40: PetscCall(MatView(A, tjbasic->viewer));
41: if (stepnum) {
42: PetscCall(TSForwardGetStages(ts, &ns, &S));
43: for (i = 0; i < ns; i++) PetscCall(MatView(S[i], tjbasic->viewer));
44: }
45: }
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode TSTrajectorySetFromOptions_Basic(TSTrajectory tj, PetscOptionItems *PetscOptionsObject)
50: {
51: PetscFunctionBegin;
52: PetscOptionsHeadBegin(PetscOptionsObject, "TS trajectory options for Basic type");
53: PetscOptionsHeadEnd();
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj, TS ts, PetscInt stepnum, PetscReal *t)
58: {
59: PetscViewer viewer;
60: char filename[PETSC_MAX_PATH_LEN];
61: Vec Sol;
62: PetscInt ns, i;
64: PetscFunctionBegin;
65: PetscCall(PetscSNPrintf(filename, sizeof(filename), tj->dirfiletemplate, stepnum));
66: PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj), filename, FILE_MODE_READ, &viewer));
67: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
68: PetscCall(TSGetSolution(ts, &Sol));
69: PetscCall(VecLoad(Sol, viewer));
70: PetscCall(PetscViewerBinaryRead(viewer, t, 1, NULL, PETSC_REAL));
71: if (stepnum && !tj->solution_only) {
72: Vec *Y;
73: PetscReal timepre;
74: PetscCall(TSGetStages(ts, &ns, &Y));
75: for (i = 0; i < ns; i++) {
76: /* 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. */
77: if (ts->stifflyaccurate && i == ns - 1) continue;
78: PetscCall(VecLoad(Y[i], viewer));
79: }
80: PetscCall(PetscViewerBinaryRead(viewer, &timepre, 1, NULL, PETSC_REAL));
81: if (tj->adjoint_solve_mode) PetscCall(TSSetTimeStep(ts, -(*t) + timepre));
82: }
83: /* Tangent linear sensitivities needed by second-order adjoint */
84: if (ts->forward_solve) {
85: if (!ts->stifflyaccurate) {
86: Mat A;
87: PetscCall(TSForwardGetSensitivities(ts, NULL, &A));
88: PetscCall(MatLoad(A, viewer));
89: }
90: if (stepnum) {
91: Mat *S;
92: PetscCall(TSForwardGetStages(ts, &ns, &S));
93: for (i = 0; i < ns; i++) PetscCall(MatLoad(S[i], viewer));
94: }
95: }
96: PetscCall(PetscViewerDestroy(&viewer));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj, TS ts)
101: {
102: MPI_Comm comm;
103: PetscMPIInt rank;
105: PetscFunctionBegin;
106: PetscCall(PetscObjectGetComm((PetscObject)tj, &comm));
107: PetscCallMPI(MPI_Comm_rank(comm, &rank));
108: if (rank == 0) {
109: char *dir = tj->dirname;
110: PetscBool flg;
112: if (!dir) {
113: char dtempname[16] = "TS-data-XXXXXX";
114: PetscCall(PetscMkdtemp(dtempname));
115: PetscCall(PetscStrallocpy(dtempname, &tj->dirname));
116: } else {
117: PetscCall(PetscTestDirectory(dir, 'w', &flg));
118: if (!flg) {
119: PetscCall(PetscTestFile(dir, 'r', &flg));
120: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_USER, "Specified path is a file - not a dir: %s", dir);
121: PetscCall(PetscMkdir(dir));
122: } else SETERRQ(comm, PETSC_ERR_SUP, "Directory %s not empty", tj->dirname);
123: }
124: }
125: PetscCall(PetscBarrier((PetscObject)tj));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
130: {
131: TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic *)tj->data;
133: PetscFunctionBegin;
134: PetscCall(PetscViewerDestroy(&tjbasic->viewer));
135: PetscCall(PetscFree(tjbasic));
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: /*MC
140: TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
142: Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed.
144: This version saves the solutions at all the stages
146: $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
148: Level: intermediate
150: .seealso: [](ch_ts), `TSTrajectoryCreate()`, `TS`, `TSTrajectory`, `TSTrajectorySetType()`, `TSTrajectorySetDirname()`, `TSTrajectorySetFile()`,
151: `TSTrajectoryType`
152: M*/
153: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj, TS ts)
154: {
155: TSTrajectory_Basic *tjbasic;
157: PetscFunctionBegin;
158: PetscCall(PetscNew(&tjbasic));
160: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)tj), &tjbasic->viewer));
161: PetscCall(PetscViewerSetType(tjbasic->viewer, PETSCVIEWERBINARY));
162: PetscCall(PetscViewerPushFormat(tjbasic->viewer, PETSC_VIEWER_NATIVE));
163: PetscCall(PetscViewerFileSetMode(tjbasic->viewer, FILE_MODE_WRITE));
164: tj->data = tjbasic;
166: tj->ops->set = TSTrajectorySet_Basic;
167: tj->ops->get = TSTrajectoryGet_Basic;
168: tj->ops->setup = TSTrajectorySetUp_Basic;
169: tj->ops->destroy = TSTrajectoryDestroy_Basic;
170: tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }