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;
30: PetscInitialize(&argc, &argv, (char *)0, help);
31: VecCreate(PETSC_COMM_WORLD, &W);
32: VecSetSizes(W, 1, PETSC_DECIDE);
33: VecSetUp(W);
34: VecDuplicate(W, &Wdot);
35: VecDuplicate(W, &W2);
36: TSCreate(PETSC_COMM_WORLD, &ts);
37: TSSetSolution(ts, W2);
38: TSSetMaxSteps(ts, 10);
39: TSSetSaveTrajectory(ts);
40: TSGetTrajectory(ts, &tj);
41: TSTrajectorySetType(tj, ts, TSTRAJECTORYBASIC);
42: TSTrajectorySetFromOptions(tj, ts);
43: TSTrajectorySetSolutionOnly(tj, PETSC_TRUE);
44: TSTrajectorySetUp(tj, ts);
45: PetscOptionsGetBool(NULL, NULL, "-check", &check, NULL);
46: PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL);
47: PetscOptionsGetRealArray(NULL, NULL, "-interptimes", times, &Nt, NULL);
48: sort = PETSC_FALSE;
49: PetscOptionsGetBool(NULL, NULL, "-sorttimes", &sort, NULL);
50: if (sort) PetscSortReal(10, TT);
51: sort = PETSC_FALSE;
52: PetscOptionsGetBool(NULL, NULL, "-sortkeys", &sort, NULL);
53: if (sort) PetscSortInt(10, II);
54: p = PetscMax(p, -p);
56: /* populate trajectory */
57: for (i = 0; i < 10; i++) {
58: VecSet(W, func(p, TT[i]));
59: TSSetStepNumber(ts, II[i]);
60: TSTrajectorySet(tj, ts, II[i], TT[i], W);
61: }
62: for (i = 0; i < Nt; i++) {
63: PetscReal testtime = times[i], serr, derr;
64: const PetscScalar *aW, *aWdot;
66: TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, W, Wdot);
67: VecGetArrayRead(W, &aW);
68: VecGetArrayRead(Wdot, &aWdot);
69: PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0]));
70: PetscPrintf(PETSC_COMM_WORLD, "df(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(dfunc(p, testtime)), (double)PetscRealPart(aWdot[0]));
71: serr = PetscAbsScalar(func(p, testtime) - aW[0]);
72: derr = PetscAbsScalar(dfunc(p, testtime) - aWdot[0]);
73: VecRestoreArrayRead(W, &aW);
74: VecRestoreArrayRead(Wdot, &aWdot);
77: }
78: for (i = Nt - 1; i >= 0; i--) {
79: PetscReal testtime = times[i], serr;
80: const PetscScalar *aW;
82: TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, W, NULL);
83: VecGetArrayRead(W, &aW);
84: PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0]));
85: serr = PetscAbsScalar(func(p, testtime) - aW[0]);
86: VecRestoreArrayRead(W, &aW);
88: }
89: for (i = Nt - 1; i >= 0; i--) {
90: PetscReal testtime = times[i], derr;
91: const PetscScalar *aWdot;
93: TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, NULL, Wdot);
94: VecGetArrayRead(Wdot, &aWdot);
95: PetscPrintf(PETSC_COMM_WORLD, "df(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(dfunc(p, testtime)), (double)PetscRealPart(aWdot[0]));
96: derr = PetscAbsScalar(dfunc(p, testtime) - aWdot[0]);
97: VecRestoreArrayRead(Wdot, &aWdot);
99: }
100: for (i = 0; i < Nt; i++) {
101: PetscReal testtime = times[i], serr, derr;
102: const PetscScalar *aW, *aWdot;
103: Vec hW, hWdot;
105: TSTrajectoryGetUpdatedHistoryVecs(tj, ts, testtime, &hW, &hWdot);
106: VecGetArrayRead(hW, &aW);
107: VecGetArrayRead(hWdot, &aWdot);
108: PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0]));
109: PetscPrintf(PETSC_COMM_WORLD, "df(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(dfunc(p, testtime)), (double)PetscRealPart(aWdot[0]));
110: serr = PetscAbsScalar(func(p, testtime) - aW[0]);
111: derr = PetscAbsScalar(dfunc(p, testtime) - aWdot[0]);
112: VecRestoreArrayRead(hW, &aW);
113: VecRestoreArrayRead(hWdot, &aWdot);
114: TSTrajectoryRestoreUpdatedHistoryVecs(tj, &hW, &hWdot);
117: }
119: /* Test on-the-fly reconstruction */
120: TSDestroy(&ts);
121: TSCreate(PETSC_COMM_WORLD, &ts);
122: TSSetSolution(ts, W2);
123: TSSetMaxSteps(ts, 10);
124: TSSetSaveTrajectory(ts);
125: TSGetTrajectory(ts, &tj);
126: TSTrajectorySetType(tj, ts, TSTRAJECTORYBASIC);
127: TSTrajectorySetFromOptions(tj, ts);
128: TSTrajectorySetSolutionOnly(tj, PETSC_TRUE);
129: TSTrajectorySetUp(tj, ts);
130: use1 = PETSC_FALSE;
131: use2 = PETSC_TRUE;
132: PetscOptionsGetBool(NULL, NULL, "-use_state", &use1, NULL);
133: PetscOptionsGetBool(NULL, NULL, "-use_der", &use2, NULL);
134: PetscSortReal(10, TT);
135: for (i = 0; i < 10; i++) {
136: TSSetStepNumber(ts, i);
137: VecSet(W, func(p, TT[i]));
138: TSTrajectorySet(tj, ts, i, TT[i], W);
139: if (i) {
140: const PetscScalar *aW, *aWdot;
141: Vec hW, hWdot;
142: PetscReal testtime = TT[i], serr, derr;
144: TSTrajectoryGetUpdatedHistoryVecs(tj, ts, testtime, use1 ? &hW : NULL, use2 ? &hWdot : NULL);
145: if (use1) {
146: VecGetArrayRead(hW, &aW);
147: PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0]));
148: serr = PetscAbsScalar(func(p, testtime) - aW[0]);
149: VecRestoreArrayRead(hW, &aW);
151: }
152: if (use2) {
153: VecGetArrayRead(hWdot, &aWdot);
154: PetscPrintf(PETSC_COMM_WORLD, "df(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(dfunc(p, testtime)), (double)PetscRealPart(aWdot[0]));
155: derr = PetscAbsScalar(dfunc(p, testtime) - aWdot[0]);
156: VecRestoreArrayRead(hWdot, &aWdot);
158: }
159: TSTrajectoryRestoreUpdatedHistoryVecs(tj, use1 ? &hW : NULL, use2 ? &hWdot : NULL);
160: }
161: }
162: TSRemoveTrajectory(ts);
163: TSDestroy(&ts);
164: VecDestroy(&W);
165: VecDestroy(&W2);
166: VecDestroy(&Wdot);
167: PetscFinalize();
168: return 0;
169: }
171: /*TEST
173: test:
174: suffix: 1
175: requires: !single
176: args: -ts_trajectory_monitor -p 3 -ts_trajectory_reconstruction_order 3 -interptimes 1,9.9,3,1.1,1.1,5.6 -check
178: test:
179: suffix: 2
180: requires: !single
181: 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
183: test:
184: suffix: 3
185: requires: !single
186: args: -ts_trajectory_monitor -p 3 -ts_trajectory_reconstruction_order 5 -interptimes 1,9.9,3,1.1,1.1,5.6 -check
188: TEST*/