Actual source code: ex13.c
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>
9: PetscScalar func(PetscInt p, PetscReal t)
10: {
11: return p ? t * func(p - 1, t) : 1.0;
12: }
13: PetscScalar dfunc(PetscInt p, PetscReal t)
14: {
15: return p > 0 ? (PetscReal)p * func(p - 1, t) : 0.0;
16: }
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;
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);
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;
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;
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;
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;
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: }
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;
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: }
173: /*TEST
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
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
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
190: TEST*/