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

  4: typedef struct {
  5:   PetscViewer viewer;
  6: } TSTrajectory_Basic;

  8: static PetscErrorCode TSTrajectorySet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
  9: {
 10:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
 11:   char               filename[PETSC_MAX_PATH_LEN];
 12:   PetscErrorCode     ierr;

 15:   PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
 16:   PetscViewerFileSetName(tjbasic->viewer,filename);
 17:   PetscViewerSetUp(tjbasic->viewer);
 18:   VecView(X,tjbasic->viewer);
 19:   PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL,PETSC_FALSE);
 20:   if (stepnum && !tj->solution_only) {
 21:     Vec       *Y;
 22:     PetscReal tprev;
 23:     PetscInt  ns,i;

 25:     TSGetStages(ts,&ns,&Y);
 26:     for (i=0;i<ns;i++) {
 27:       VecView(Y[i],tjbasic->viewer);
 28:     }
 29:     TSGetPrevTime(ts,&tprev);
 30:     PetscViewerBinaryWrite(tjbasic->viewer,&tprev,1,PETSC_REAL,PETSC_FALSE);
 31:   }
 32:   return(0);
 33: }

 35: static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj)
 36: {

 40:   PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type");
 41:   PetscOptionsTail();
 42:   return(0);
 43: }

 45: static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
 46: {
 47:   PetscViewer    viewer;
 48:   char           filename[PETSC_MAX_PATH_LEN];
 50:   Vec            Sol;

 53:   PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
 54:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 55:   TSGetSolution(ts,&Sol);
 56:   PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
 57:   VecLoad(Sol,viewer);
 58:   PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);
 59:   if (stepnum != 0 && !tj->solution_only) {
 60:     Vec         *Y;
 61:     PetscInt    Nr,i;
 62:     PetscReal   timepre;

 64:     TSGetStages(ts,&Nr,&Y);
 65:     for (i=0;i<Nr ;i++) {
 66:       VecLoad(Y[i],viewer);
 67:     }
 68:     PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);
 69:     if (tj->adjoint_solve_mode) {
 70:       TSSetTimeStep(ts,-(*t)+timepre);
 71:     }
 72:   }
 73:   PetscViewerDestroy(&viewer);
 74:   return(0);
 75: }

 77: static PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)
 78: {
 79:   MPI_Comm       comm;
 80:   PetscMPIInt    rank;

 84:   PetscObjectGetComm((PetscObject)tj,&comm);
 85:   MPI_Comm_rank(comm,&rank);
 86:   if (!rank) {
 87:     const char* dir = tj->dirname;
 88:     PetscBool   flg;

 90:     /* I don't like running PetscRMTree on a directory */
 91:     PetscTestDirectory(dir,'w',&flg);
 92:     if (!flg) {
 93:       PetscTestFile(dir,'r',&flg);
 94:       if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir);
 95:       PetscMkdir(dir);
 96:     } else SETERRQ1(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname);
 97:   }
 98:   PetscBarrier((PetscObject)tj);
 99:   return(0);
100: }

102: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
103: {
104:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
105:   PetscErrorCode     ierr;

108:   PetscViewerDestroy(&tjbasic->viewer);
109:   PetscFree(tjbasic);
110:   return(0);
111: }

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

116:       Saves each timestep into a seperate file named SA-data/SA-%06d.bin. The file name can be changed.

118:       This version saves the solutions at all the stages

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

122:   Level: intermediate

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

126: M*/
127: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
128: {
129:   TSTrajectory_Basic *tjbasic;
130:   PetscErrorCode     ierr;

133:   PetscNew(&tjbasic);

135:   PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);
136:   PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);
137:   PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);
138:   PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);
139:   tj->data = tjbasic;

141:   tj->ops->set            = TSTrajectorySet_Basic;
142:   tj->ops->get            = TSTrajectoryGet_Basic;
143:   tj->ops->setup          = TSTrajectorySetUp_Basic;
144:   tj->ops->destroy        = TSTrajectoryDestroy_Basic;
145:   tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
146:   return(0);
147: }