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 { 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 } 17 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; 28 29 PetscFunctionBeginUser; 30 PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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); 56 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; 66 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; 82 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; 93 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; 105 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 } 119 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; 145 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 } 172 173 /*TEST 174 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 179 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 184 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 189 190 TEST*/ 191