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