1 static char help[] = "Tests TSTrajectoryGetVecs. \n\n";
2 /*
3 This example tests TSTrajectory and the ability of TSTrajectoryGetVecs
4 to reconstructs states and derivatives via interpolation (if necessary).
5 It also tests TSTrajectory{Get|Restore}UpdatedHistoryVecs
6 */
7 #include <petscts.h>
8
func(PetscInt p,PetscReal t)9 PetscScalar func(PetscInt p, PetscReal t)
10 {
11 return p ? t * func(p - 1, t) : 1.0;
12 }
dfunc(PetscInt p,PetscReal t)13 PetscScalar dfunc(PetscInt p, PetscReal t)
14 {
15 return p > 0 ? (PetscReal)p * func(p - 1, t) : 0.0;
16 }
17
main(int argc,char ** argv)18 int main(int argc, char **argv)
19 {
20 TS ts;
21 Vec W, W2, Wdot;
22 TSTrajectory tj;
23 PetscReal times[10], tol = PETSC_SMALL;
24 PetscReal TT[10] = {2, 9, 1, 3, 6, 7, 5, 10, 4, 8};
25 PetscInt i, p = 1, Nt = 10;
26 PetscInt II[10] = {1, 4, 9, 2, 3, 6, 5, 8, 0, 7};
27 PetscBool sort, use1, use2, check = PETSC_FALSE;
28
29 PetscFunctionBeginUser;
30 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
31 PetscCall(VecCreate(PETSC_COMM_WORLD, &W));
32 PetscCall(VecSetSizes(W, 1, PETSC_DECIDE));
33 PetscCall(VecSetUp(W));
34 PetscCall(VecDuplicate(W, &Wdot));
35 PetscCall(VecDuplicate(W, &W2));
36 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
37 PetscCall(TSSetSolution(ts, W2));
38 PetscCall(TSSetMaxTime(ts, 10.));
39 PetscCall(TSSetMaxSteps(ts, 10));
40 PetscCall(TSSetSaveTrajectory(ts));
41 PetscCall(TSGetTrajectory(ts, &tj));
42 PetscCall(TSTrajectorySetType(tj, ts, TSTRAJECTORYBASIC));
43 PetscCall(TSTrajectorySetFromOptions(tj, ts));
44 PetscCall(TSTrajectorySetSolutionOnly(tj, PETSC_TRUE));
45 PetscCall(TSTrajectorySetUp(tj, ts));
46 PetscCall(PetscOptionsGetBool(NULL, NULL, "-check", &check, NULL));
47 PetscCall(PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL));
48 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-interptimes", times, &Nt, NULL));
49 sort = PETSC_FALSE;
50 PetscCall(PetscOptionsGetBool(NULL, NULL, "-sorttimes", &sort, NULL));
51 if (sort) PetscCall(PetscSortReal(10, TT));
52 sort = PETSC_FALSE;
53 PetscCall(PetscOptionsGetBool(NULL, NULL, "-sortkeys", &sort, NULL));
54 if (sort) PetscCall(PetscSortInt(10, II));
55 p = PetscMax(p, -p);
56
57 /* populate trajectory */
58 for (i = 0; i < 10; i++) {
59 PetscCall(VecSet(W, func(p, TT[i])));
60 PetscCall(TSSetStepNumber(ts, II[i]));
61 PetscCall(TSTrajectorySet(tj, ts, II[i], TT[i], W));
62 }
63 for (i = 0; i < Nt; i++) {
64 PetscReal testtime = times[i], serr, derr;
65 const PetscScalar *aW, *aWdot;
66
67 PetscCall(TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, W, Wdot));
68 PetscCall(VecGetArrayRead(W, &aW));
69 PetscCall(VecGetArrayRead(Wdot, &aWdot));
70 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0])));
71 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "df(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(dfunc(p, testtime)), (double)PetscRealPart(aWdot[0])));
72 serr = PetscAbsScalar(func(p, testtime) - aW[0]);
73 derr = PetscAbsScalar(dfunc(p, testtime) - aWdot[0]);
74 PetscCall(VecRestoreArrayRead(W, &aW));
75 PetscCall(VecRestoreArrayRead(Wdot, &aWdot));
76 PetscCheck(!check || serr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state %g > %g", (double)serr, (double)tol);
77 PetscCheck(!check || derr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state derivative %g > %g", (double)derr, (double)tol);
78 }
79 for (i = Nt - 1; i >= 0; i--) {
80 PetscReal testtime = times[i], serr;
81 const PetscScalar *aW;
82
83 PetscCall(TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, W, NULL));
84 PetscCall(VecGetArrayRead(W, &aW));
85 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0])));
86 serr = PetscAbsScalar(func(p, testtime) - aW[0]);
87 PetscCall(VecRestoreArrayRead(W, &aW));
88 PetscCheck(!check || serr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state %g > %g", (double)serr, (double)tol);
89 }
90 for (i = Nt - 1; i >= 0; i--) {
91 PetscReal testtime = times[i], derr;
92 const PetscScalar *aWdot;
93
94 PetscCall(TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, NULL, Wdot));
95 PetscCall(VecGetArrayRead(Wdot, &aWdot));
96 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "df(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(dfunc(p, testtime)), (double)PetscRealPart(aWdot[0])));
97 derr = PetscAbsScalar(dfunc(p, testtime) - aWdot[0]);
98 PetscCall(VecRestoreArrayRead(Wdot, &aWdot));
99 PetscCheck(!check || derr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state derivative %g > %g", (double)derr, (double)tol);
100 }
101 for (i = 0; i < Nt; i++) {
102 PetscReal testtime = times[i], serr, derr;
103 const PetscScalar *aW, *aWdot;
104 Vec hW, hWdot;
105
106 PetscCall(TSTrajectoryGetUpdatedHistoryVecs(tj, ts, testtime, &hW, &hWdot));
107 PetscCall(VecGetArrayRead(hW, &aW));
108 PetscCall(VecGetArrayRead(hWdot, &aWdot));
109 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0])));
110 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "df(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(dfunc(p, testtime)), (double)PetscRealPart(aWdot[0])));
111 serr = PetscAbsScalar(func(p, testtime) - aW[0]);
112 derr = PetscAbsScalar(dfunc(p, testtime) - aWdot[0]);
113 PetscCall(VecRestoreArrayRead(hW, &aW));
114 PetscCall(VecRestoreArrayRead(hWdot, &aWdot));
115 PetscCall(TSTrajectoryRestoreUpdatedHistoryVecs(tj, &hW, &hWdot));
116 PetscCheck(!check || serr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state %g > %g", (double)serr, (double)tol);
117 PetscCheck(!check || derr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state derivative %g > %g", (double)derr, (double)tol);
118 }
119
120 /* Test on-the-fly reconstruction */
121 PetscCall(TSDestroy(&ts));
122 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
123 PetscCall(TSSetSolution(ts, W2));
124 PetscCall(TSSetMaxTime(ts, 10.));
125 PetscCall(TSSetMaxSteps(ts, 10));
126 PetscCall(TSSetSaveTrajectory(ts));
127 PetscCall(TSGetTrajectory(ts, &tj));
128 PetscCall(TSTrajectorySetType(tj, ts, TSTRAJECTORYBASIC));
129 PetscCall(TSTrajectorySetFromOptions(tj, ts));
130 PetscCall(TSTrajectorySetSolutionOnly(tj, PETSC_TRUE));
131 PetscCall(TSTrajectorySetUp(tj, ts));
132 use1 = PETSC_FALSE;
133 use2 = PETSC_TRUE;
134 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_state", &use1, NULL));
135 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_der", &use2, NULL));
136 PetscCall(PetscSortReal(10, TT));
137 for (i = 0; i < 10; i++) {
138 PetscCall(TSSetStepNumber(ts, i));
139 PetscCall(VecSet(W, func(p, TT[i])));
140 PetscCall(TSTrajectorySet(tj, ts, i, TT[i], W));
141 if (i) {
142 const PetscScalar *aW, *aWdot;
143 Vec hW, hWdot;
144 PetscReal testtime = TT[i], serr, derr;
145
146 PetscCall(TSTrajectoryGetUpdatedHistoryVecs(tj, ts, testtime, use1 ? &hW : NULL, use2 ? &hWdot : NULL));
147 if (use1) {
148 PetscCall(VecGetArrayRead(hW, &aW));
149 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0])));
150 serr = PetscAbsScalar(func(p, testtime) - aW[0]);
151 PetscCall(VecRestoreArrayRead(hW, &aW));
152 PetscCheck(!check || serr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state %g > %g", (double)serr, (double)tol);
153 }
154 if (use2) {
155 PetscCall(VecGetArrayRead(hWdot, &aWdot));
156 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "df(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(dfunc(p, testtime)), (double)PetscRealPart(aWdot[0])));
157 derr = PetscAbsScalar(dfunc(p, testtime) - aWdot[0]);
158 PetscCall(VecRestoreArrayRead(hWdot, &aWdot));
159 PetscCheck(!check || i < p || derr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state derivative %g > %g", (double)derr, (double)tol);
160 }
161 PetscCall(TSTrajectoryRestoreUpdatedHistoryVecs(tj, use1 ? &hW : NULL, use2 ? &hWdot : NULL));
162 }
163 }
164 PetscCall(TSRemoveTrajectory(ts));
165 PetscCall(TSDestroy(&ts));
166 PetscCall(VecDestroy(&W));
167 PetscCall(VecDestroy(&W2));
168 PetscCall(VecDestroy(&Wdot));
169 PetscCall(PetscFinalize());
170 return 0;
171 }
172
173 /*TEST
174
175 test:
176 suffix: 1
177 requires: !single
178 args: -ts_trajectory_monitor -p 3 -ts_trajectory_reconstruction_order 3 -interptimes 1,9.9,3,1.1,1.1,5.6 -check
179
180 test:
181 suffix: 2
182 requires: !single
183 args: -sortkeys -ts_trajectory_monitor -ts_trajectory_type memory -p 3 -ts_trajectory_reconstruction_order 3 -ts_trajectory_adjointmode 0 -interptimes 1,9.9,3,1.1,1.1,5.6 -check
184
185 test:
186 suffix: 3
187 requires: !single
188 args: -ts_trajectory_monitor -p 3 -ts_trajectory_reconstruction_order 5 -interptimes 1,9.9,3,1.1,1.1,5.6 -check
189
190 TEST*/
191