Actual source code: ex13.c

  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>

  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; }

 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;

 23:   PetscInitialize(&argc,&argv,(char*)0,help);
 24:   VecCreate(PETSC_COMM_WORLD,&W);
 25:   VecSetSizes(W,1,PETSC_DECIDE);
 26:   VecSetUp(W);
 27:   VecDuplicate(W,&Wdot);
 28:   VecDuplicate(W,&W2);
 29:   TSCreate(PETSC_COMM_WORLD,&ts);
 30:   TSSetSolution(ts,W2);
 31:   TSSetMaxSteps(ts,10);
 32:   TSSetSaveTrajectory(ts);
 33:   TSGetTrajectory(ts,&tj);
 34:   TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
 35:   TSTrajectorySetFromOptions(tj,ts);
 36:   TSTrajectorySetSolutionOnly(tj,PETSC_TRUE);
 37:   TSTrajectorySetUp(tj,ts);
 38:   PetscOptionsGetBool(NULL,NULL,"-check",&check,NULL);
 39:   PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
 40:   PetscOptionsGetRealArray(NULL,NULL,"-interptimes",times,&Nt,NULL);
 41:   sort = PETSC_FALSE;
 42:   PetscOptionsGetBool(NULL,NULL,"-sorttimes",&sort,NULL);
 43:   if (sort) {
 44:     PetscSortReal(10,TT);
 45:   }
 46:   sort = PETSC_FALSE;
 47:   PetscOptionsGetBool(NULL,NULL,"-sortkeys",&sort,NULL);
 48:   if (sort) {
 49:     PetscSortInt(10,II);
 50:   }
 51:   p = PetscMax(p,-p);

 53:   /* populate trajectory */
 54:   for (i = 0; i < 10; i++) {
 55:     VecSet(W,func(p,TT[i]));
 56:     TSSetStepNumber(ts,II[i]);
 57:     TSTrajectorySet(tj,ts,II[i],TT[i],W);
 58:   }
 59:   for (i = 0; i < Nt; i++) {
 60:     PetscReal testtime = times[i], serr, derr;
 61:     const PetscScalar *aW,*aWdot;

 63:     TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&testtime,W,Wdot);
 64:     VecGetArrayRead(W,&aW);
 65:     VecGetArrayRead(Wdot,&aWdot);
 66:     PetscPrintf(PETSC_COMM_WORLD," f(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(func(p,testtime)),(double)PetscRealPart(aW[0]));
 67:     PetscPrintf(PETSC_COMM_WORLD,"df(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(dfunc(p,testtime)),(double)PetscRealPart(aWdot[0]));
 68:     serr = PetscAbsScalar(func(p,testtime)-aW[0]);
 69:     derr = PetscAbsScalar(dfunc(p,testtime)-aWdot[0]);
 70:     VecRestoreArrayRead(W,&aW);
 71:     VecRestoreArrayRead(Wdot,&aWdot);
 74:   }
 75:   for (i = Nt-1; i >= 0; i--) {
 76:     PetscReal         testtime = times[i], serr;
 77:     const PetscScalar *aW;

 79:     TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&testtime,W,NULL);
 80:     VecGetArrayRead(W,&aW);
 81:     PetscPrintf(PETSC_COMM_WORLD," f(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(func(p,testtime)),(double)PetscRealPart(aW[0]));
 82:     serr = PetscAbsScalar(func(p,testtime)-aW[0]);
 83:     VecRestoreArrayRead(W,&aW);
 85:   }
 86:   for (i = Nt-1; i >= 0; i--) {
 87:     PetscReal         testtime = times[i], derr;
 88:     const PetscScalar *aWdot;

 90:     TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&testtime,NULL,Wdot);
 91:     VecGetArrayRead(Wdot,&aWdot);
 92:     PetscPrintf(PETSC_COMM_WORLD,"df(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(dfunc(p,testtime)),(double)PetscRealPart(aWdot[0]));
 93:     derr = PetscAbsScalar(dfunc(p,testtime)-aWdot[0]);
 94:     VecRestoreArrayRead(Wdot,&aWdot);
 96:   }
 97:   for (i = 0; i < Nt; i++) {
 98:     PetscReal         testtime = times[i], serr, derr;
 99:     const PetscScalar *aW,*aWdot;
100:     Vec               hW,hWdot;

102:     TSTrajectoryGetUpdatedHistoryVecs(tj,ts,testtime,&hW,&hWdot);
103:     VecGetArrayRead(hW,&aW);
104:     VecGetArrayRead(hWdot,&aWdot);
105:     PetscPrintf(PETSC_COMM_WORLD," f(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(func(p,testtime)),(double)PetscRealPart(aW[0]));
106:     PetscPrintf(PETSC_COMM_WORLD,"df(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(dfunc(p,testtime)),(double)PetscRealPart(aWdot[0]));
107:     serr = PetscAbsScalar(func(p,testtime)-aW[0]);
108:     derr = PetscAbsScalar(dfunc(p,testtime)-aWdot[0]);
109:     VecRestoreArrayRead(hW,&aW);
110:     VecRestoreArrayRead(hWdot,&aWdot);
111:     TSTrajectoryRestoreUpdatedHistoryVecs(tj,&hW,&hWdot);
114:   }

116:   /* Test on-the-fly reconstruction */
117:   TSDestroy(&ts);
118:   TSCreate(PETSC_COMM_WORLD,&ts);
119:   TSSetSolution(ts,W2);
120:   TSSetMaxSteps(ts,10);
121:   TSSetSaveTrajectory(ts);
122:   TSGetTrajectory(ts,&tj);
123:   TSTrajectorySetType(tj,ts,TSTRAJECTORYBASIC);
124:   TSTrajectorySetFromOptions(tj,ts);
125:   TSTrajectorySetSolutionOnly(tj,PETSC_TRUE);
126:   TSTrajectorySetUp(tj,ts);
127:   use1 = PETSC_FALSE;
128:   use2 = PETSC_TRUE;
129:   PetscOptionsGetBool(NULL,NULL,"-use_state",&use1,NULL);
130:   PetscOptionsGetBool(NULL,NULL,"-use_der",&use2,NULL);
131:   PetscSortReal(10,TT);
132:   for (i = 0; i < 10; i++) {

134:     TSSetStepNumber(ts,i);
135:     VecSet(W,func(p,TT[i]));
136:     TSTrajectorySet(tj,ts,i,TT[i],W);
137:     if (i) {
138:       const PetscScalar *aW,*aWdot;
139:       Vec               hW,hWdot;
140:       PetscReal         testtime = TT[i], serr, derr;

142:       TSTrajectoryGetUpdatedHistoryVecs(tj,ts,testtime,use1 ? &hW : NULL,use2 ? &hWdot : NULL);
143:       if (use1) {
144:         VecGetArrayRead(hW,&aW);
145:         PetscPrintf(PETSC_COMM_WORLD," f(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(func(p,testtime)),(double)PetscRealPart(aW[0]));
146:         serr = PetscAbsScalar(func(p,testtime)-aW[0]);
147:         VecRestoreArrayRead(hW,&aW);
149:       }
150:       if (use2) {
151:         VecGetArrayRead(hWdot,&aWdot);
152:         PetscPrintf(PETSC_COMM_WORLD,"df(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(dfunc(p,testtime)),(double)PetscRealPart(aWdot[0]));
153:         derr = PetscAbsScalar(dfunc(p,testtime)-aWdot[0]);
154:         VecRestoreArrayRead(hWdot,&aWdot);
156:       }
157:       TSTrajectoryRestoreUpdatedHistoryVecs(tj,use1 ? &hW : NULL,use2 ? &hWdot : NULL);
158:     }
159:   }
160:   TSRemoveTrajectory(ts);
161:   TSDestroy(&ts);
162:   VecDestroy(&W);
163:   VecDestroy(&W2);
164:   VecDestroy(&Wdot);
165:   PetscFinalize();
166:   return 0;
167: }

169: /*TEST

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

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

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

186: TEST*/