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