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