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: PetscSNPrintf(filename, sizeof(filename), tj->dirfiletemplate, stepnum);
19: PetscViewerFileSetName(tjbasic->viewer, filename); /* this triggers PetscViewer to be set up again */
20: PetscViewerSetUp(tjbasic->viewer);
21: VecView(X, tjbasic->viewer);
22: PetscViewerBinaryWrite(tjbasic->viewer, &time, 1, PETSC_REAL);
23: if (stepnum && !tj->solution_only) {
24: Vec *Y;
25: PetscReal tprev;
26: 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: VecView(Y[i], tjbasic->viewer);
31: }
32: TSGetPrevTime(ts, &tprev);
33: 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: TSForwardGetSensitivities(ts, NULL, &A);
40: MatView(A, tjbasic->viewer);
41: if (stepnum) {
42: TSForwardGetStages(ts, &ns, &S);
43: for (i = 0; i < ns; i++) MatView(S[i], tjbasic->viewer);
44: }
45: }
46: return 0;
47: }
49: static PetscErrorCode TSTrajectorySetFromOptions_Basic(TSTrajectory tj, PetscOptionItems *PetscOptionsObject)
50: {
51: PetscOptionsHeadBegin(PetscOptionsObject, "TS trajectory options for Basic type");
52: PetscOptionsHeadEnd();
53: return 0;
54: }
56: static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj, TS ts, PetscInt stepnum, PetscReal *t)
57: {
58: PetscViewer viewer;
59: char filename[PETSC_MAX_PATH_LEN];
60: Vec Sol;
61: PetscInt ns, i;
63: PetscSNPrintf(filename, sizeof(filename), tj->dirfiletemplate, stepnum);
64: PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj), filename, FILE_MODE_READ, &viewer);
65: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
66: TSGetSolution(ts, &Sol);
67: VecLoad(Sol, viewer);
68: PetscViewerBinaryRead(viewer, t, 1, NULL, PETSC_REAL);
69: if (stepnum && !tj->solution_only) {
70: Vec *Y;
71: PetscReal timepre;
72: TSGetStages(ts, &ns, &Y);
73: for (i = 0; i < ns; i++) {
74: /* 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. */
75: if (ts->stifflyaccurate && i == ns - 1) continue;
76: VecLoad(Y[i], viewer);
77: }
78: PetscViewerBinaryRead(viewer, &timepre, 1, NULL, PETSC_REAL);
79: if (tj->adjoint_solve_mode) TSSetTimeStep(ts, -(*t) + timepre);
80: }
81: /* Tangent linear sensitivities needed by second-order adjoint */
82: if (ts->forward_solve) {
83: if (!ts->stifflyaccurate) {
84: Mat A;
85: TSForwardGetSensitivities(ts, NULL, &A);
86: MatLoad(A, viewer);
87: }
88: if (stepnum) {
89: Mat *S;
90: TSForwardGetStages(ts, &ns, &S);
91: for (i = 0; i < ns; i++) MatLoad(S[i], viewer);
92: }
93: }
94: PetscViewerDestroy(&viewer);
95: return 0;
96: }
98: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj, TS ts)
99: {
100: MPI_Comm comm;
101: PetscMPIInt rank;
103: PetscObjectGetComm((PetscObject)tj, &comm);
104: MPI_Comm_rank(comm, &rank);
105: if (rank == 0) {
106: char *dir = tj->dirname;
107: PetscBool flg;
109: if (!dir) {
110: char dtempname[16] = "TS-data-XXXXXX";
111: PetscMkdtemp(dtempname);
112: PetscStrallocpy(dtempname, &tj->dirname);
113: } else {
114: PetscTestDirectory(dir, 'w', &flg);
115: if (!flg) {
116: PetscTestFile(dir, 'r', &flg);
118: PetscMkdir(dir);
119: } else SETERRQ(comm, PETSC_ERR_SUP, "Directory %s not empty", tj->dirname);
120: }
121: }
122: PetscBarrier((PetscObject)tj);
123: return 0;
124: }
126: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
127: {
128: TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic *)tj->data;
130: PetscViewerDestroy(&tjbasic->viewer);
131: PetscFree(tjbasic);
132: return 0;
133: }
135: /*MC
136: TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
138: Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed.
140: This version saves the solutions at all the stages
142: $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
144: Level: intermediate
146: .seealso: [](chapter_ts), `TSTrajectoryCreate()`, `TS`, `TSTrajectory`, `TSTrajectorySetType()`, `TSTrajectorySetDirname()`, `TSTrajectorySetFile()`,
147: `TSTrajectoryType`
148: M*/
149: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj, TS ts)
150: {
151: TSTrajectory_Basic *tjbasic;
153: PetscNew(&tjbasic);
155: PetscViewerCreate(PetscObjectComm((PetscObject)tj), &tjbasic->viewer);
156: PetscViewerSetType(tjbasic->viewer, PETSCVIEWERBINARY);
157: PetscViewerPushFormat(tjbasic->viewer, PETSC_VIEWER_NATIVE);
158: PetscViewerFileSetMode(tjbasic->viewer, FILE_MODE_WRITE);
159: tj->data = tjbasic;
161: tj->ops->set = TSTrajectorySet_Basic;
162: tj->ops->get = TSTrajectoryGet_Basic;
163: tj->ops->setup = TSTrajectorySetUp_Basic;
164: tj->ops->destroy = TSTrajectoryDestroy_Basic;
165: tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
166: return 0;
167: }