Actual source code: ex13.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      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;

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

     54:   /* populate trajectory */
     55:   for (i = 0; i < 10; i++) {
     56:     VecSet(W,func(p,TT[i]));
     57:     TSSetStepNumber(ts,II[i]);
     58:     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;

     64:     TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&testtime,W,Wdot);
     65:     VecGetArrayRead(W,&aW);
     66:     VecGetArrayRead(Wdot,&aWdot);
     67:     PetscPrintf(PETSC_COMM_WORLD," f(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(func(p,testtime)),(double)PetscRealPart(aW[0]));
     68:     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:     VecRestoreArrayRead(W,&aW);
     72:     VecRestoreArrayRead(Wdot,&aWdot);
     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;

     80:     TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&testtime,W,NULL);
     81:     VecGetArrayRead(W,&aW);
     82:     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:     VecRestoreArrayRead(W,&aW);
     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;

     91:     TSTrajectoryGetVecs(tj,ts,PETSC_DECIDE,&testtime,NULL,Wdot);
     92:     VecGetArrayRead(Wdot,&aWdot);
     93:     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:     VecRestoreArrayRead(Wdot,&aWdot);
     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;

    103:     TSTrajectoryGetUpdatedHistoryVecs(tj,ts,testtime,&hW,&hWdot);
    104:     VecGetArrayRead(hW,&aW);
    105:     VecGetArrayRead(hWdot,&aWdot);
    106:     PetscPrintf(PETSC_COMM_WORLD," f(%g) = %g (reconstructed %g)\n",(double)testtime,(double)PetscRealPart(func(p,testtime)),(double)PetscRealPart(aW[0]));
    107:     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:     VecRestoreArrayRead(hW,&aW);
    111:     VecRestoreArrayRead(hWdot,&aWdot);
    112:     TSTrajectoryRestoreUpdatedHistoryVecs(tj,&hW,&hWdot);
    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:   }

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

    135:     TSSetStepNumber(ts,i);
    136:     VecSet(W,func(p,TT[i]));
    137:     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;

    143:       TSTrajectoryGetUpdatedHistoryVecs(tj,ts,testtime,use1 ? &hW : NULL,use2 ? &hWdot : NULL);
    144:       if (use1) {
    145:         VecGetArrayRead(hW,&aW);
    146:         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:         VecRestoreArrayRead(hW,&aW);
    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:         VecGetArrayRead(hWdot,&aWdot);
    153:         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:         VecRestoreArrayRead(hWdot,&aWdot);
    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:       TSTrajectoryRestoreUpdatedHistoryVecs(tj,use1 ? &hW : NULL,use2 ? &hWdot : NULL);
    159:     }
    160:   }
    161:   TSDestroy(&ts);
    162:   VecDestroy(&W);
    163:   VecDestroy(&W2);
    164:   VecDestroy(&Wdot);
    165:   PetscFinalize();
    166:   return ierr;
    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*/