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++) {
44: MatView(S[i],tjbasic->viewer);
45: }
46: }
47: }
48: return 0;
49: }
51: static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj)
52: {
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];
62: Vec Sol;
63: PetscInt ns,i;
65: PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
66: PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj),filename,FILE_MODE_READ,&viewer);
67: PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
68: TSGetSolution(ts,&Sol);
69: VecLoad(Sol,viewer);
70: PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);
71: if (stepnum && !tj->solution_only) {
72: Vec *Y;
73: PetscReal timepre;
74: TSGetStages(ts,&ns,&Y);
75: for (i=0; i<ns; i++) {
76: /* 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. */
77: if (ts->stifflyaccurate && i == ns-1) continue;
78: VecLoad(Y[i],viewer);
79: }
80: PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);
81: if (tj->adjoint_solve_mode) {
82: TSSetTimeStep(ts,-(*t)+timepre);
83: }
84: }
85: /* Tangent linear sensitivities needed by second-order adjoint */
86: if (ts->forward_solve) {
88: if (!ts->stifflyaccurate) {
89: Mat A;
90: TSForwardGetSensitivities(ts,NULL,&A);
91: MatLoad(A,viewer);
92: }
93: if (stepnum) {
94: Mat *S;
95: TSForwardGetStages(ts,&ns,&S);
96: for (i=0; i<ns; i++) {
97: MatLoad(S[i],viewer);
98: }
99: }
100: }
101: PetscViewerDestroy(&viewer);
102: return 0;
103: }
105: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)
106: {
107: MPI_Comm comm;
108: PetscMPIInt rank;
110: PetscObjectGetComm((PetscObject)tj,&comm);
111: MPI_Comm_rank(comm,&rank);
112: if (rank == 0) {
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);
125: PetscMkdir(dir);
126: } else SETERRQ(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;
137: PetscViewerDestroy(&tjbasic->viewer);
138: PetscFree(tjbasic);
139: return 0;
140: }
142: /*MC
143: TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
145: Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed.
147: This version saves the solutions at all the stages
149: $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
151: Level: intermediate
153: .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile()
155: M*/
156: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
157: {
158: TSTrajectory_Basic *tjbasic;
160: PetscNew(&tjbasic);
162: PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);
163: PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);
164: PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);
165: PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);
166: tj->data = tjbasic;
168: tj->ops->set = TSTrajectorySet_Basic;
169: tj->ops->get = TSTrajectoryGet_Basic;
170: tj->ops->setup = TSTrajectorySetUp_Basic;
171: tj->ops->destroy = TSTrajectoryDestroy_Basic;
172: tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
173: return 0;
174: }