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