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*/