xref: /petsc/src/ts/trajectory/impls/basic/trajbasic.c (revision 36a9e3b9f6565ce1252c167e0dc4a4cf71b0f2ec)
1 
2 #include <petsc/private/tsimpl.h>        /*I "petscts.h"  I*/
3 
4 typedef struct {
5   PetscViewer viewer;
6 } TSTrajectory_Basic;
7 
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   PetscErrorCode     ierr;
13 
14   PetscFunctionBegin;
15   ierr = PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);CHKERRQ(ierr);
16   ierr = PetscViewerFileSetName(tjbasic->viewer,filename);CHKERRQ(ierr);
17   ierr = PetscViewerSetUp(tjbasic->viewer);CHKERRQ(ierr);
18   ierr = VecView(X,tjbasic->viewer);CHKERRQ(ierr);
19   ierr = PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
20   if (stepnum && !tj->solution_only) {
21     Vec       *Y;
22     PetscReal tprev;
23     PetscInt  ns,i;
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,PETSC_FALSE);CHKERRQ(ierr);
31   }
32   PetscFunctionReturn(0);
33 }
34 
35 static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj)
36 {
37   PetscErrorCode ierr;
38 
39   PetscFunctionBegin;
40   ierr = PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type");CHKERRQ(ierr);
41   ierr = PetscOptionsTail();CHKERRQ(ierr);
42   PetscFunctionReturn(0);
43 }
44 
45 static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
46 {
47   PetscViewer    viewer;
48   char           filename[PETSC_MAX_PATH_LEN];
49   PetscErrorCode ierr;
50   Vec            Sol;
51 
52   PetscFunctionBegin;
53   ierr = PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);CHKERRQ(ierr);
54   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
55   ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr);
56   ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);CHKERRQ(ierr);
57   ierr = VecLoad(Sol,viewer);CHKERRQ(ierr);
58   ierr = PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);CHKERRQ(ierr);
59   if (stepnum != 0 && !tj->solution_only) {
60     Vec         *Y;
61     PetscInt    Nr,i;
62     PetscReal   timepre;
63 
64     ierr = TSGetStages(ts,&Nr,&Y);CHKERRQ(ierr);
65     for (i=0;i<Nr ;i++) {
66       ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr);
67     }
68     ierr = PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);CHKERRQ(ierr);
69     if (tj->adjoint_solve_mode) {
70       ierr = TSSetTimeStep(ts,-(*t)+timepre);CHKERRQ(ierr);
71     }
72   }
73   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
74   PetscFunctionReturn(0);
75 }
76 
77 static PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts)
78 {
79   MPI_Comm       comm;
80   PetscMPIInt    rank;
81   PetscErrorCode ierr;
82 
83   PetscFunctionBegin;
84   ierr = PetscObjectGetComm((PetscObject)tj,&comm);CHKERRQ(ierr);
85   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
86   if (!rank) {
87     const char* dir = tj->dirname;
88     PetscBool   flg;
89 
90     /* I don't like running PetscRMTree on a directory */
91     ierr = PetscTestDirectory(dir,'w',&flg);CHKERRQ(ierr);
92     if (!flg) {
93       ierr = PetscTestFile(dir,'r',&flg);CHKERRQ(ierr);
94       if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir);
95       ierr = PetscMkdir(dir);CHKERRQ(ierr);
96     } else SETERRQ1(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname);
97   }
98   ierr = PetscBarrier((PetscObject)tj);CHKERRQ(ierr);
99   PetscFunctionReturn(0);
100 }
101 
102 static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
103 {
104   TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data;
105   PetscErrorCode     ierr;
106 
107   PetscFunctionBegin;
108   ierr = PetscViewerDestroy(&tjbasic->viewer);CHKERRQ(ierr);
109   ierr = PetscFree(tjbasic);CHKERRQ(ierr);
110   PetscFunctionReturn(0);
111 }
112 
113 /*MC
114       TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
115 
116       Saves each timestep into a seperate file named SA-data/SA-%06d.bin. The file name can be changed.
117 
118       This version saves the solutions at all the stages
119 
120       $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
121 
122   Level: intermediate
123 
124 .seealso:  TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile()
125 
126 M*/
127 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
128 {
129   TSTrajectory_Basic *tjbasic;
130   PetscErrorCode     ierr;
131 
132   PetscFunctionBegin;
133   ierr = PetscNew(&tjbasic);CHKERRQ(ierr);
134 
135   ierr = PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);CHKERRQ(ierr);
136   ierr = PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
137   ierr = PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);CHKERRQ(ierr);
138   ierr = PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
139   tj->data = tjbasic;
140 
141   tj->ops->set            = TSTrajectorySet_Basic;
142   tj->ops->get            = TSTrajectoryGet_Basic;
143   tj->ops->setup          = TSTrajectorySetUp_Basic;
144   tj->ops->destroy        = TSTrajectoryDestroy_Basic;
145   tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic;
146   PetscFunctionReturn(0);
147 }
148