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;
17: PetscErrorCode ierr;
20: PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
21: PetscViewerFileSetName(tjbasic->viewer,filename); /* this triggers PetscViewer to be set up again */
22: PetscViewerSetUp(tjbasic->viewer);
23: VecView(X,tjbasic->viewer);
24: PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL);
25: if (stepnum && !tj->solution_only) {
26: Vec *Y;
27: PetscReal tprev;
28: TSGetStages(ts,&ns,&Y);
29: for (i=0; i<ns; i++) {
30: /* 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. */
31: if (ts->stifflyaccurate && i == ns-1) continue;
32: VecView(Y[i],tjbasic->viewer);
33: }
34: TSGetPrevTime(ts,&tprev);
35: PetscViewerBinaryWrite(tjbasic->viewer,&tprev,1,PETSC_REAL);
36: }
37: /* Tangent linear sensitivities needed by second-order adjoint */
38: if (ts->forward_solve) {
39: Mat A,*S;
41: TSForwardGetSensitivities(ts,NULL,&A);
42: MatView(A,tjbasic->viewer);
43: if (stepnum) {
44: TSForwardGetStages(ts,&ns,&S);
45: for (i=0; i<ns; i++) {
46: MatView(S[i],tjbasic->viewer);
47: }
48: }
49: }
50: return(0);
51: }
53: static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj)
54: {
58: PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type");
59: PetscOptionsTail();
60: return(0);
61: }
63: static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
64: {
65: PetscViewer viewer;
66: char filename[PETSC_MAX_PATH_LEN];
68: Vec Sol;
69: PetscInt ns,i;
72: PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
73: PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj),filename,FILE_MODE_READ,&viewer);
74: PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
75: TSGetSolution(ts,&Sol);
76: VecLoad(Sol,viewer);
77: PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);
78: if (stepnum && !tj->solution_only) {
79: Vec *Y;
80: PetscReal timepre;
81: TSGetStages(ts,&ns,&Y);
82: for (i=0; i<ns; i++) {
83: /* 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. */
84: if (ts->stifflyaccurate && i == ns-1) continue;
85: VecLoad(Y[i],viewer);
86: }
87: PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);
88: if (tj->adjoint_solve_mode) {
89: TSSetTimeStep(ts,-(*t)+timepre);
90: }
91: }
92: /* Tangent linear sensitivities needed by second-order adjoint */
93: if (ts->forward_solve) {
95: if (!ts->stifflyaccurate) {
96: Mat A;
97: TSForwardGetSensitivities(ts,NULL,&A);
98: MatLoad(A,viewer);
99: }
100: if (stepnum) {
101: Mat *S;
102: TSForwardGetStages(ts,&ns,&S);
103: for (i=0; i<ns; i++) {
104: MatLoad(S[i],viewer);
105: }
106: }
107: }
108: PetscViewerDestroy(&viewer);
109: return(0);
110: }
112: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)
113: {
114: MPI_Comm comm;
115: PetscMPIInt rank;
119: PetscObjectGetComm((PetscObject)tj,&comm);
120: MPI_Comm_rank(comm,&rank);
121: if (rank == 0) {
122: char *dir = tj->dirname;
123: PetscBool flg;
125: if (!dir) {
126: char dtempname[16] = "TS-data-XXXXXX";
127: PetscMkdtemp(dtempname);
128: PetscStrallocpy(dtempname,&tj->dirname);
129: } else {
130: PetscTestDirectory(dir,'w',&flg);
131: if (!flg) {
132: PetscTestFile(dir,'r',&flg);
133: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir);
134: PetscMkdir(dir);
135: } else SETERRQ1(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname);
136: }
137: }
138: PetscBarrier((PetscObject)tj);
139: return(0);
140: }
142: static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
143: {
144: TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
145: PetscErrorCode ierr;
148: PetscViewerDestroy(&tjbasic->viewer);
149: PetscFree(tjbasic);
150: return(0);
151: }
153: /*MC
154: TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
156: Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed.
158: This version saves the solutions at all the stages
160: $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
162: Level: intermediate
164: .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile()
166: M*/
167: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
168: {
169: TSTrajectory_Basic *tjbasic;
170: PetscErrorCode ierr;
173: PetscNew(&tjbasic);
175: PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);
176: PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);
177: PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);
178: PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);
179: tj->data = tjbasic;
181: tj->ops->set = TSTrajectorySet_Basic;
182: tj->ops->get = TSTrajectoryGet_Basic;
183: tj->ops->setup = TSTrajectorySetUp_Basic;
184: tj->ops->destroy = TSTrajectoryDestroy_Basic;
185: tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
186: return(0);
187: }