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(TSSetMaxSteps(ts, 10)); 39 PetscCall(TSSetSaveTrajectory(ts)); 40 PetscCall(TSGetTrajectory(ts, &tj)); 41 PetscCall(TSTrajectorySetType(tj, ts, TSTRAJECTORYBASIC)); 42 PetscCall(TSTrajectorySetFromOptions(tj, ts)); 43 PetscCall(TSTrajectorySetSolutionOnly(tj, PETSC_TRUE)); 44 PetscCall(TSTrajectorySetUp(tj, ts)); 45 PetscCall(PetscOptionsGetBool(NULL, NULL, "-check", &check, NULL)); 46 PetscCall(PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL)); 47 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-interptimes", times, &Nt, NULL)); 48 sort = PETSC_FALSE; 49 PetscCall(PetscOptionsGetBool(NULL, NULL, "-sorttimes", &sort, NULL)); 50 if (sort) PetscCall(PetscSortReal(10, TT)); 51 sort = PETSC_FALSE; 52 PetscCall(PetscOptionsGetBool(NULL, NULL, "-sortkeys", &sort, NULL)); 53 if (sort) PetscCall(PetscSortInt(10, II)); 54 p = PetscMax(p, -p); 55 56 /* populate trajectory */ 57 for (i = 0; i < 10; i++) { 58 PetscCall(VecSet(W, func(p, TT[i]))); 59 PetscCall(TSSetStepNumber(ts, II[i])); 60 PetscCall(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; 65 66 PetscCall(TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, W, Wdot)); 67 PetscCall(VecGetArrayRead(W, &aW)); 68 PetscCall(VecGetArrayRead(Wdot, &aWdot)); 69 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0]))); 70 PetscCall(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 PetscCall(VecRestoreArrayRead(W, &aW)); 74 PetscCall(VecRestoreArrayRead(Wdot, &aWdot)); 75 PetscCheck(!check || serr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state %g > %g", (double)serr, (double)tol); 76 PetscCheck(!check || derr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state derivative %g > %g", (double)derr, (double)tol); 77 } 78 for (i = Nt - 1; i >= 0; i--) { 79 PetscReal testtime = times[i], serr; 80 const PetscScalar *aW; 81 82 PetscCall(TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, W, NULL)); 83 PetscCall(VecGetArrayRead(W, &aW)); 84 PetscCall(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 PetscCall(VecRestoreArrayRead(W, &aW)); 87 PetscCheck(!check || serr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state %g > %g", (double)serr, (double)tol); 88 } 89 for (i = Nt - 1; i >= 0; i--) { 90 PetscReal testtime = times[i], derr; 91 const PetscScalar *aWdot; 92 93 PetscCall(TSTrajectoryGetVecs(tj, ts, PETSC_DECIDE, &testtime, NULL, Wdot)); 94 PetscCall(VecGetArrayRead(Wdot, &aWdot)); 95 PetscCall(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 PetscCall(VecRestoreArrayRead(Wdot, &aWdot)); 98 PetscCheck(!check || derr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state derivative %g > %g", (double)derr, (double)tol); 99 } 100 for (i = 0; i < Nt; i++) { 101 PetscReal testtime = times[i], serr, derr; 102 const PetscScalar *aW, *aWdot; 103 Vec hW, hWdot; 104 105 PetscCall(TSTrajectoryGetUpdatedHistoryVecs(tj, ts, testtime, &hW, &hWdot)); 106 PetscCall(VecGetArrayRead(hW, &aW)); 107 PetscCall(VecGetArrayRead(hWdot, &aWdot)); 108 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " f(%g) = %g (reconstructed %g)\n", (double)testtime, (double)PetscRealPart(func(p, testtime)), (double)PetscRealPart(aW[0]))); 109 PetscCall(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 PetscCall(VecRestoreArrayRead(hW, &aW)); 113 PetscCall(VecRestoreArrayRead(hWdot, &aWdot)); 114 PetscCall(TSTrajectoryRestoreUpdatedHistoryVecs(tj, &hW, &hWdot)); 115 PetscCheck(!check || serr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state %g > %g", (double)serr, (double)tol); 116 PetscCheck(!check || derr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state derivative %g > %g", (double)derr, (double)tol); 117 } 118 119 /* Test on-the-fly reconstruction */ 120 PetscCall(TSDestroy(&ts)); 121 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 122 PetscCall(TSSetSolution(ts, W2)); 123 PetscCall(TSSetMaxSteps(ts, 10)); 124 PetscCall(TSSetSaveTrajectory(ts)); 125 PetscCall(TSGetTrajectory(ts, &tj)); 126 PetscCall(TSTrajectorySetType(tj, ts, TSTRAJECTORYBASIC)); 127 PetscCall(TSTrajectorySetFromOptions(tj, ts)); 128 PetscCall(TSTrajectorySetSolutionOnly(tj, PETSC_TRUE)); 129 PetscCall(TSTrajectorySetUp(tj, ts)); 130 use1 = PETSC_FALSE; 131 use2 = PETSC_TRUE; 132 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_state", &use1, NULL)); 133 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_der", &use2, NULL)); 134 PetscCall(PetscSortReal(10, TT)); 135 for (i = 0; i < 10; i++) { 136 PetscCall(TSSetStepNumber(ts, i)); 137 PetscCall(VecSet(W, func(p, TT[i]))); 138 PetscCall(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; 143 144 PetscCall(TSTrajectoryGetUpdatedHistoryVecs(tj, ts, testtime, use1 ? &hW : NULL, use2 ? &hWdot : NULL)); 145 if (use1) { 146 PetscCall(VecGetArrayRead(hW, &aW)); 147 PetscCall(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 PetscCall(VecRestoreArrayRead(hW, &aW)); 150 PetscCheck(!check || serr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state %g > %g", (double)serr, (double)tol); 151 } 152 if (use2) { 153 PetscCall(VecGetArrayRead(hWdot, &aWdot)); 154 PetscCall(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 PetscCall(VecRestoreArrayRead(hWdot, &aWdot)); 157 PetscCheck(!check || i < p || derr <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in state derivative %g > %g", (double)derr, (double)tol); 158 } 159 PetscCall(TSTrajectoryRestoreUpdatedHistoryVecs(tj, use1 ? &hW : NULL, use2 ? &hWdot : NULL)); 160 } 161 } 162 PetscCall(TSRemoveTrajectory(ts)); 163 PetscCall(TSDestroy(&ts)); 164 PetscCall(VecDestroy(&W)); 165 PetscCall(VecDestroy(&W2)); 166 PetscCall(VecDestroy(&Wdot)); 167 PetscCall(PetscFinalize()); 168 return 0; 169 } 170 171 /*TEST 172 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 177 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 182 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 187 188 TEST*/ 189