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