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

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

 50: static PetscErrorCode TSTrajectorySetFromOptions_Basic(TSTrajectory tj, PetscOptionItems *PetscOptionsObject)
 51: {
 52:   PetscFunctionBegin;
 53:   PetscOptionsHeadBegin(PetscOptionsObject, "TS trajectory options for Basic type");
 54:   PetscOptionsHeadEnd();
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 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];
 62:   Vec         Sol;
 63:   PetscInt    ns, i;

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

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

106:   PetscFunctionBegin;
107:   PetscCall(PetscObjectGetComm((PetscObject)tj, &comm));
108:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
109:   if (rank == 0) {
110:     char     *dir = tj->dirname;
111:     PetscBool flg;

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

130: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
131: {
132:   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic *)tj->data;

134:   PetscFunctionBegin;
135:   PetscCall(PetscViewerDestroy(&tjbasic->viewer));
136:   PetscCall(PetscFree(tjbasic));
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

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

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

145:       This version saves the solutions at all the stages

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

149:   Level: intermediate

151: .seealso: [](ch_ts), `TSTrajectoryCreate()`, `TS`, `TSTrajectory`, `TSTrajectorySetType()`, `TSTrajectorySetDirname()`, `TSTrajectorySetFile()`,
152:           `TSTrajectoryType`
153: M*/
154: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj, TS ts)
155: {
156:   TSTrajectory_Basic *tjbasic;

158:   PetscFunctionBegin;
159:   PetscCall(PetscNew(&tjbasic));

161:   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)tj), &tjbasic->viewer));
162:   PetscCall(PetscViewerSetType(tjbasic->viewer, PETSCVIEWERBINARY));
163:   PetscCall(PetscViewerPushFormat(tjbasic->viewer, PETSC_VIEWER_NATIVE));
164:   PetscCall(PetscViewerFileSetMode(tjbasic->viewer, FILE_MODE_WRITE));
165:   tj->data = tjbasic;

167:   tj->ops->set            = TSTrajectorySet_Basic;
168:   tj->ops->get            = TSTrajectoryGet_Basic;
169:   tj->ops->setup          = TSTrajectorySetUp_Basic;
170:   tj->ops->destroy        = TSTrajectoryDestroy_Basic;
171:   tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }