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;

 18:   PetscSNPrintf(filename, sizeof(filename), tj->dirfiletemplate, stepnum);
 19:   PetscViewerFileSetName(tjbasic->viewer, filename); /* this triggers PetscViewer to be set up again */
 20:   PetscViewerSetUp(tjbasic->viewer);
 21:   VecView(X, tjbasic->viewer);
 22:   PetscViewerBinaryWrite(tjbasic->viewer, &time, 1, PETSC_REAL);
 23:   if (stepnum && !tj->solution_only) {
 24:     Vec      *Y;
 25:     PetscReal tprev;
 26:     TSGetStages(ts, &ns, &Y);
 27:     for (i = 0; i < ns; i++) {
 28:       /* 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. */
 29:       if (ts->stifflyaccurate && i == ns - 1) continue;
 30:       VecView(Y[i], tjbasic->viewer);
 31:     }
 32:     TSGetPrevTime(ts, &tprev);
 33:     PetscViewerBinaryWrite(tjbasic->viewer, &tprev, 1, PETSC_REAL);
 34:   }
 35:   /* Tangent linear sensitivities needed by second-order adjoint */
 36:   if (ts->forward_solve) {
 37:     Mat A, *S;

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

 49: static PetscErrorCode TSTrajectorySetFromOptions_Basic(TSTrajectory tj, PetscOptionItems *PetscOptionsObject)
 50: {
 51:   PetscOptionsHeadBegin(PetscOptionsObject, "TS trajectory options for Basic type");
 52:   PetscOptionsHeadEnd();
 53:   return 0;
 54: }

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

 63:   PetscSNPrintf(filename, sizeof(filename), tj->dirfiletemplate, stepnum);
 64:   PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj), filename, FILE_MODE_READ, &viewer);
 65:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
 66:   TSGetSolution(ts, &Sol);
 67:   VecLoad(Sol, viewer);
 68:   PetscViewerBinaryRead(viewer, t, 1, NULL, PETSC_REAL);
 69:   if (stepnum && !tj->solution_only) {
 70:     Vec      *Y;
 71:     PetscReal timepre;
 72:     TSGetStages(ts, &ns, &Y);
 73:     for (i = 0; i < ns; i++) {
 74:       /* 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. */
 75:       if (ts->stifflyaccurate && i == ns - 1) continue;
 76:       VecLoad(Y[i], viewer);
 77:     }
 78:     PetscViewerBinaryRead(viewer, &timepre, 1, NULL, PETSC_REAL);
 79:     if (tj->adjoint_solve_mode) TSSetTimeStep(ts, -(*t) + timepre);
 80:   }
 81:   /* Tangent linear sensitivities needed by second-order adjoint */
 82:   if (ts->forward_solve) {
 83:     if (!ts->stifflyaccurate) {
 84:       Mat A;
 85:       TSForwardGetSensitivities(ts, NULL, &A);
 86:       MatLoad(A, viewer);
 87:     }
 88:     if (stepnum) {
 89:       Mat *S;
 90:       TSForwardGetStages(ts, &ns, &S);
 91:       for (i = 0; i < ns; i++) MatLoad(S[i], viewer);
 92:     }
 93:   }
 94:   PetscViewerDestroy(&viewer);
 95:   return 0;
 96: }

 98: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj, TS ts)
 99: {
100:   MPI_Comm    comm;
101:   PetscMPIInt rank;

103:   PetscObjectGetComm((PetscObject)tj, &comm);
104:   MPI_Comm_rank(comm, &rank);
105:   if (rank == 0) {
106:     char     *dir = tj->dirname;
107:     PetscBool flg;

109:     if (!dir) {
110:       char dtempname[16] = "TS-data-XXXXXX";
111:       PetscMkdtemp(dtempname);
112:       PetscStrallocpy(dtempname, &tj->dirname);
113:     } else {
114:       PetscTestDirectory(dir, 'w', &flg);
115:       if (!flg) {
116:         PetscTestFile(dir, 'r', &flg);
118:         PetscMkdir(dir);
119:       } else SETERRQ(comm, PETSC_ERR_SUP, "Directory %s not empty", tj->dirname);
120:     }
121:   }
122:   PetscBarrier((PetscObject)tj);
123:   return 0;
124: }

126: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
127: {
128:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic *)tj->data;

130:   PetscViewerDestroy(&tjbasic->viewer);
131:   PetscFree(tjbasic);
132:   return 0;
133: }

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

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

140:       This version saves the solutions at all the stages

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

144:   Level: intermediate

146: .seealso: [](chapter_ts), `TSTrajectoryCreate()`, `TS`, `TSTrajectory`, `TSTrajectorySetType()`, `TSTrajectorySetDirname()`, `TSTrajectorySetFile()`,
147:           `TSTrajectoryType`
148: M*/
149: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj, TS ts)
150: {
151:   TSTrajectory_Basic *tjbasic;

153:   PetscNew(&tjbasic);

155:   PetscViewerCreate(PetscObjectComm((PetscObject)tj), &tjbasic->viewer);
156:   PetscViewerSetType(tjbasic->viewer, PETSCVIEWERBINARY);
157:   PetscViewerPushFormat(tjbasic->viewer, PETSC_VIEWER_NATIVE);
158:   PetscViewerFileSetMode(tjbasic->viewer, FILE_MODE_WRITE);
159:   tj->data = tjbasic;

161:   tj->ops->set            = TSTrajectorySet_Basic;
162:   tj->ops->get            = TSTrajectoryGet_Basic;
163:   tj->ops->setup          = TSTrajectorySetUp_Basic;
164:   tj->ops->destroy        = TSTrajectoryDestroy_Basic;
165:   tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
166:   return 0;
167: }