1 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
2
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 */
TSTrajectorySet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)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;
16
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;
38
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 }
48
TSTrajectorySetFromOptions_Basic(TSTrajectory tj,PetscOptionItems PetscOptionsObject)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 }
56
TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal * t)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;
63
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 }
99
TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)100 PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj, TS ts)
101 {
102 MPI_Comm comm;
103 PetscMPIInt rank;
104
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;
111
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 PetscCheck(!flg, comm, PETSC_ERR_SUP, "Directory %s not empty", tj->dirname);
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 }
123 }
124 PetscCall(PetscBarrier((PetscObject)tj));
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
TSTrajectoryDestroy_Basic(TSTrajectory tj)128 static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
129 {
130 TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic *)tj->data;
131
132 PetscFunctionBegin;
133 PetscCall(PetscViewerDestroy(&tjbasic->viewer));
134 PetscCall(PetscFree(tjbasic));
135 PetscFunctionReturn(PETSC_SUCCESS);
136 }
137
138 /*MC
139 TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
140
141 Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed.
142
143 This version saves the solutions at all the stages
144
145 $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
146
147 Level: intermediate
148
149 .seealso: [](ch_ts), `TSTrajectoryCreate()`, `TS`, `TSTrajectory`, `TSTrajectorySetType()`, `TSTrajectorySetDirname()`, `TSTrajectorySetFile()`,
150 `TSTrajectoryType`
151 M*/
TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)152 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj, TS ts)
153 {
154 TSTrajectory_Basic *tjbasic;
155
156 PetscFunctionBegin;
157 PetscCall(PetscNew(&tjbasic));
158
159 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)tj), &tjbasic->viewer));
160 PetscCall(PetscViewerSetType(tjbasic->viewer, PETSCVIEWERBINARY));
161 PetscCall(PetscViewerPushFormat(tjbasic->viewer, PETSC_VIEWER_NATIVE));
162 PetscCall(PetscViewerFileSetMode(tjbasic->viewer, FILE_MODE_WRITE));
163 tj->data = tjbasic;
164
165 tj->ops->set = TSTrajectorySet_Basic;
166 tj->ops->get = TSTrajectoryGet_Basic;
167 tj->ops->setup = TSTrajectorySetUp_Basic;
168 tj->ops->destroy = TSTrajectoryDestroy_Basic;
169 tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172