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:   PetscInt           ns,i;
 13:   PetscErrorCode     ierr;

 16:   PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
 17:   PetscViewerFileSetName(tjbasic->viewer,filename); /* this triggers PetscViewer to be set up again */
 18:   PetscViewerSetUp(tjbasic->viewer);
 19:   VecView(X,tjbasic->viewer);
 20:   PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL);
 21:   if (stepnum && !tj->solution_only) {
 22:     Vec       *Y;
 23:     PetscReal tprev;

 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);
 31:   }
 32:   /* Tangent linear sensitivities needed by second-order adjoint */
 33:   if (ts->forward_solve) {
 34:     Mat A,*S;

 36:     TSForwardGetSensitivities(ts,NULL,&A);
 37:     MatView(A,tjbasic->viewer);
 38:     if (stepnum) {
 39:       TSForwardGetStages(ts,&ns,&S);
 40:       for (i=0; i<ns; i++) {
 41:         MatView(S[i],tjbasic->viewer);
 42:       }
 43:     }
 44:   }
 45:   return(0);
 46: }

 48: static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj)
 49: {

 53:   PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type");
 54:   PetscOptionsTail();
 55:   return(0);
 56: }

 58: static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
 59: {
 60:   PetscViewer    viewer;
 61:   char           filename[PETSC_MAX_PATH_LEN];
 63:   Vec            Sol;
 64:   PetscInt       ns,i;

 67:   PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
 68:   PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj),filename,FILE_MODE_READ,&viewer);
 69:   TSGetSolution(ts,&Sol);
 70:   PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
 71:   VecLoad(Sol,viewer);
 72:   PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);
 73:   if (stepnum && !tj->solution_only) {
 74:     Vec       *Y;
 75:     PetscReal timepre;

 77:     TSGetStages(ts,&ns,&Y);
 78:     for (i=0; i<ns; i++) {
 79:       VecLoad(Y[i],viewer);
 80:     }
 81:     PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);
 82:     if (tj->adjoint_solve_mode) {
 83:       TSSetTimeStep(ts,-(*t)+timepre);
 84:     }
 85:   }
 86:   /* Tangent linear sensitivities needed by second-order adjoint */
 87:   if (ts->forward_solve) {
 88:     Mat A,*S;

 90:     TSForwardGetSensitivities(ts,NULL,&A);
 91:     MatLoad(A,viewer);
 92:     if (stepnum) {
 93:       TSForwardGetStages(ts,&ns,&S);
 94:       for (i=0; i<ns; i++) {
 95:         MatLoad(S[i],viewer);
 96:       }
 97:     }
 98:   }
 99:   PetscViewerDestroy(&viewer);
100:   return(0);
101: }

103: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)
104: {
105:   MPI_Comm       comm;
106:   PetscMPIInt    rank;

110:   PetscObjectGetComm((PetscObject)tj,&comm);
111:   MPI_Comm_rank(comm,&rank);
112:   if (!rank) {
113:     char      *dir = tj->dirname;
114:     PetscBool flg;

116:     if (!dir) {
117:       char dtempname[16] = "TS-data-XXXXXX";
118:       PetscMkdtemp(dtempname);
119:       PetscStrallocpy(dtempname,&tj->dirname);
120:     } else {
121:       PetscTestDirectory(dir,'w',&flg);
122:       if (!flg) {
123:         PetscTestFile(dir,'r',&flg);
124:         if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir);
125:         PetscMkdir(dir);
126:       } else SETERRQ1(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname);
127:     }
128:   }
129:   PetscBarrier((PetscObject)tj);
130:   return(0);
131: }

133: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
134: {
135:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
136:   PetscErrorCode     ierr;

139:   PetscViewerDestroy(&tjbasic->viewer);
140:   PetscFree(tjbasic);
141:   return(0);
142: }

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

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

149:       This version saves the solutions at all the stages

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

153:   Level: intermediate

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

157: M*/
158: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
159: {
160:   TSTrajectory_Basic *tjbasic;
161:   PetscErrorCode     ierr;

164:   PetscNew(&tjbasic);

166:   PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);
167:   PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);
168:   PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);
169:   PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);
170:   tj->data = tjbasic;

172:   tj->ops->set            = TSTrajectorySet_Basic;
173:   tj->ops->get            = TSTrajectoryGet_Basic;
174:   tj->ops->setup          = TSTrajectorySetUp_Basic;
175:   tj->ops->destroy        = TSTrajectoryDestroy_Basic;
176:   tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
177:   return(0);
178: }