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