Actual source code: trajbasic.c


  2: #include <petsc/private/tsimpl.h>

  4: typedef struct {
  5:   PetscViewer viewer;
  6: } TSTrajectory_Basic;
  7: /*
  8:   For n-th time step, TSTrajectorySet_Basic always saves the solution X(t_n) and the current time t_n,
  9:   and optionally saves the stage values Y[] between t_{n-1} and t_n, the previous time t_{n-1}, and
 10:   forward stage sensitivities S[] = dY[]/dp.
 11: */
 12: static PetscErrorCode TSTrajectorySet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
 13: {
 14:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
 15:   char               filename[PETSC_MAX_PATH_LEN];
 16:   PetscInt           ns,i;
 17:   PetscErrorCode     ierr;

 20:   PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
 21:   PetscViewerFileSetName(tjbasic->viewer,filename); /* this triggers PetscViewer to be set up again */
 22:   PetscViewerSetUp(tjbasic->viewer);
 23:   VecView(X,tjbasic->viewer);
 24:   PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL);
 25:   if (stepnum && !tj->solution_only) {
 26:     Vec       *Y;
 27:     PetscReal tprev;
 28:     TSGetStages(ts,&ns,&Y);
 29:     for (i=0; i<ns; i++) {
 30:       /* For stiffly accurate TS methods, the last stage Y[ns-1] is the same as the solution X, thus does not need to be saved again. */
 31:       if (ts->stifflyaccurate && i == ns-1) continue;
 32:       VecView(Y[i],tjbasic->viewer);
 33:     }
 34:     TSGetPrevTime(ts,&tprev);
 35:     PetscViewerBinaryWrite(tjbasic->viewer,&tprev,1,PETSC_REAL);
 36:   }
 37:   /* Tangent linear sensitivities needed by second-order adjoint */
 38:   if (ts->forward_solve) {
 39:     Mat A,*S;

 41:     TSForwardGetSensitivities(ts,NULL,&A);
 42:     MatView(A,tjbasic->viewer);
 43:     if (stepnum) {
 44:       TSForwardGetStages(ts,&ns,&S);
 45:       for (i=0; i<ns; i++) {
 46:         MatView(S[i],tjbasic->viewer);
 47:       }
 48:     }
 49:   }
 50:   return(0);
 51: }

 53: static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj)
 54: {

 58:   PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type");
 59:   PetscOptionsTail();
 60:   return(0);
 61: }

 63: static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
 64: {
 65:   PetscViewer    viewer;
 66:   char           filename[PETSC_MAX_PATH_LEN];
 68:   Vec            Sol;
 69:   PetscInt       ns,i;

 72:   PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
 73:   PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj),filename,FILE_MODE_READ,&viewer);
 74:   PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
 75:   TSGetSolution(ts,&Sol);
 76:   VecLoad(Sol,viewer);
 77:   PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);
 78:   if (stepnum && !tj->solution_only) {
 79:     Vec       *Y;
 80:     PetscReal timepre;
 81:     TSGetStages(ts,&ns,&Y);
 82:     for (i=0; i<ns; i++) {
 83:       /* For stiffly accurate TS methods, the last stage Y[ns-1] is the same as the solution X, thus does not need to be loaded again. */
 84:       if (ts->stifflyaccurate && i == ns-1) continue;
 85:       VecLoad(Y[i],viewer);
 86:     }
 87:     PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);
 88:     if (tj->adjoint_solve_mode) {
 89:       TSSetTimeStep(ts,-(*t)+timepre);
 90:     }
 91:   }
 92:   /* Tangent linear sensitivities needed by second-order adjoint */
 93:   if (ts->forward_solve) {

 95:     if (!ts->stifflyaccurate) {
 96:       Mat A;
 97:       TSForwardGetSensitivities(ts,NULL,&A);
 98:       MatLoad(A,viewer);
 99:     }
100:     if (stepnum) {
101:       Mat *S;
102:       TSForwardGetStages(ts,&ns,&S);
103:       for (i=0; i<ns; i++) {
104:         MatLoad(S[i],viewer);
105:       }
106:     }
107:   }
108:   PetscViewerDestroy(&viewer);
109:   return(0);
110: }

112: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)
113: {
114:   MPI_Comm       comm;
115:   PetscMPIInt    rank;

119:   PetscObjectGetComm((PetscObject)tj,&comm);
120:   MPI_Comm_rank(comm,&rank);
121:   if (rank == 0) {
122:     char      *dir = tj->dirname;
123:     PetscBool flg;

125:     if (!dir) {
126:       char dtempname[16] = "TS-data-XXXXXX";
127:       PetscMkdtemp(dtempname);
128:       PetscStrallocpy(dtempname,&tj->dirname);
129:     } else {
130:       PetscTestDirectory(dir,'w',&flg);
131:       if (!flg) {
132:         PetscTestFile(dir,'r',&flg);
133:         if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir);
134:         PetscMkdir(dir);
135:       } else SETERRQ1(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname);
136:     }
137:   }
138:   PetscBarrier((PetscObject)tj);
139:   return(0);
140: }

142: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
143: {
144:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
145:   PetscErrorCode     ierr;

148:   PetscViewerDestroy(&tjbasic->viewer);
149:   PetscFree(tjbasic);
150:   return(0);
151: }

153: /*MC
154:       TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file

156:       Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed.

158:       This version saves the solutions at all the stages

160:       $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format

162:   Level: intermediate

164: .seealso:  TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile()

166: M*/
167: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
168: {
169:   TSTrajectory_Basic *tjbasic;
170:   PetscErrorCode     ierr;

173:   PetscNew(&tjbasic);

175:   PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);
176:   PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);
177:   PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);
178:   PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);
179:   tj->data = tjbasic;

181:   tj->ops->set            = TSTrajectorySet_Basic;
182:   tj->ops->get            = TSTrajectoryGet_Basic;
183:   tj->ops->setup          = TSTrajectorySetUp_Basic;
184:   tj->ops->destroy        = TSTrajectoryDestroy_Basic;
185:   tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
186:   return(0);
187: }