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