Actual source code: trajbasic.c
petsc-3.13.6 2020-09-29
2: #include <petsc/private/tsimpl.h>
4: typedef struct {
5: PetscViewer viewer;
6: } TSTrajectory_Basic;
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;
16: PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
17: PetscViewerFileSetName(tjbasic->viewer,filename); /* this triggers PetscViewer to be set up again */
18: PetscViewerSetUp(tjbasic->viewer);
19: VecView(X,tjbasic->viewer);
20: PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL);
21: if (stepnum && !tj->solution_only) {
22: Vec *Y;
23: PetscReal tprev;
25: TSGetStages(ts,&ns,&Y);
26: for (i=0; i<ns; i++) {
27: VecView(Y[i],tjbasic->viewer);
28: }
29: TSGetPrevTime(ts,&tprev);
30: PetscViewerBinaryWrite(tjbasic->viewer,&tprev,1,PETSC_REAL);
31: }
32: /* Tangent linear sensitivities needed by second-order adjoint */
33: if (ts->forward_solve) {
34: Mat A,*S;
36: TSForwardGetSensitivities(ts,NULL,&A);
37: MatView(A,tjbasic->viewer);
38: if (stepnum) {
39: TSForwardGetStages(ts,&ns,&S);
40: for (i=0; i<ns; i++) {
41: MatView(S[i],tjbasic->viewer);
42: }
43: }
44: }
45: return(0);
46: }
48: static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj)
49: {
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];
63: Vec Sol;
64: PetscInt ns,i;
67: PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);
68: PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj),filename,FILE_MODE_READ,&viewer);
69: TSGetSolution(ts,&Sol);
70: PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
71: VecLoad(Sol,viewer);
72: PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);
73: if (stepnum && !tj->solution_only) {
74: Vec *Y;
75: PetscReal timepre;
77: TSGetStages(ts,&ns,&Y);
78: for (i=0; i<ns; i++) {
79: VecLoad(Y[i],viewer);
80: }
81: PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);
82: if (tj->adjoint_solve_mode) {
83: TSSetTimeStep(ts,-(*t)+timepre);
84: }
85: }
86: /* Tangent linear sensitivities needed by second-order adjoint */
87: if (ts->forward_solve) {
88: Mat A,*S;
90: TSForwardGetSensitivities(ts,NULL,&A);
91: MatLoad(A,viewer);
92: if (stepnum) {
93: TSForwardGetStages(ts,&ns,&S);
94: for (i=0; i<ns; i++) {
95: MatLoad(S[i],viewer);
96: }
97: }
98: }
99: PetscViewerDestroy(&viewer);
100: return(0);
101: }
103: PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)
104: {
105: MPI_Comm comm;
106: PetscMPIInt rank;
110: PetscObjectGetComm((PetscObject)tj,&comm);
111: MPI_Comm_rank(comm,&rank);
112: if (!rank) {
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);
124: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir);
125: PetscMkdir(dir);
126: } else SETERRQ1(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;
136: PetscErrorCode ierr;
139: PetscViewerDestroy(&tjbasic->viewer);
140: PetscFree(tjbasic);
141: return(0);
142: }
144: /*MC
145: TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
147: Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed.
149: This version saves the solutions at all the stages
151: $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
153: Level: intermediate
155: .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile()
157: M*/
158: PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
159: {
160: TSTrajectory_Basic *tjbasic;
161: PetscErrorCode ierr;
164: PetscNew(&tjbasic);
166: PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);
167: PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);
168: PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);
169: PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);
170: tj->data = tjbasic;
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: return(0);
178: }