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