xref: /petsc/src/ts/trajectory/impls/basic/trajbasic.c (revision e3aab17eafbaa3898314c4c7870cc364cef2d1b4)
1 
2 #include <petsc/private/tsimpl.h>        /*I "petscts.h"  I*/
3 
4 static PetscErrorCode OutputBIN(MPI_Comm comm,const char *filename,PetscViewer *viewer)
5 {
6   PetscErrorCode ierr;
7 
8   PetscFunctionBegin;
9   ierr = PetscViewerCreate(comm,viewer);CHKERRQ(ierr);
10   ierr = PetscViewerSetType(*viewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
11   ierr = PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
12   ierr = PetscViewerFileSetName(*viewer,filename);CHKERRQ(ierr);
13   PetscFunctionReturn(0);
14 }
15 
16 static PetscErrorCode TSTrajectorySet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
17 {
18   PetscViewer    viewer;
19   PetscInt       ns,i;
20   Vec            *Y;
21   char           filename[PETSC_MAX_PATH_LEN],ftemplate[PETSC_MAX_PATH_LEN];
22   PetscReal      tprev;
23   PetscErrorCode ierr;
24   MPI_Comm       comm;
25 
26   PetscFunctionBegin;
27   ierr = PetscSNPrintf(ftemplate,sizeof(ftemplate),"%s%s",tj->dirname,tj->filetemplate);CHKERRQ(ierr);
28   ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr);
29   ierr = TSGetStepNumber(ts,&stepnum);CHKERRQ(ierr);
30   if (stepnum == 0) {
31     PetscMPIInt rank;
32     ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
33     if (!rank) {
34       ierr = PetscRMTree(tj->dirname);CHKERRQ(ierr);
35       ierr = PetscMkdir(tj->dirname);CHKERRQ(ierr);
36     }
37     ierr = PetscSNPrintf(filename,sizeof(filename),ftemplate,stepnum);CHKERRQ(ierr);
38     ierr = OutputBIN(comm,filename,&viewer);CHKERRQ(ierr);
39     ierr = VecView(X,viewer);CHKERRQ(ierr);
40     ierr = PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
41     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
42     PetscFunctionReturn(0);
43   }
44   ierr = PetscSNPrintf(filename,sizeof(filename),ftemplate,stepnum);CHKERRQ(ierr);
45   ierr = OutputBIN(comm,filename,&viewer);CHKERRQ(ierr);
46   ierr = VecView(X,viewer);CHKERRQ(ierr);
47   ierr = PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
48 
49   ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr);
50   for (i=0;i<ns;i++) {
51     ierr = VecView(Y[i],viewer);CHKERRQ(ierr);
52   }
53 
54   ierr = TSGetPrevTime(ts,&tprev);CHKERRQ(ierr);
55   ierr = PetscViewerBinaryWrite(viewer,&tprev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
56 
57   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
58   PetscFunctionReturn(0);
59 }
60 
61 static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
62 {
63   Vec            Sol,*Y;
64   PetscInt       Nr,i;
65   PetscViewer    viewer;
66   PetscReal      timepre;
67   char           filename[PETSC_MAX_PATH_LEN],ftemplate[PETSC_MAX_PATH_LEN];
68   PetscErrorCode ierr;
69 
70   PetscFunctionBegin;
71   ierr = PetscSNPrintf(ftemplate,sizeof(ftemplate),"%s%s",tj->dirname,tj->filetemplate);CHKERRQ(ierr);
72   ierr = PetscSNPrintf(filename,sizeof filename,ftemplate,stepnum);CHKERRQ(ierr);
73   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
74 
75   ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr);
76   ierr = VecLoad(Sol,viewer);CHKERRQ(ierr);
77 
78   ierr = PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);CHKERRQ(ierr);
79 
80   if (stepnum != 0) {
81     ierr = TSGetStages(ts,&Nr,&Y);CHKERRQ(ierr);
82     for (i=0;i<Nr ;i++) {
83       ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr);
84     }
85     ierr = PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);CHKERRQ(ierr);
86     ierr = TSSetTimeStep(ts,-(*t)+timepre);CHKERRQ(ierr);
87   }
88 
89   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
90   PetscFunctionReturn(0);
91 }
92 
93 static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj)
94 {
95   PetscErrorCode ierr;
96   PetscMPIInt    rank;
97   MPI_Comm       comm;
98 
99   PetscFunctionBegin;
100   if (!tj->keepfiles) {
101     ierr = PetscObjectGetComm((PetscObject)tj,&comm);CHKERRQ(ierr);
102     ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
103     if (!rank) {
104       ierr = PetscRMTree(tj->dirname);CHKERRQ(ierr);
105     }
106   }
107   PetscFunctionReturn(0);
108 }
109 
110 /*MC
111       TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
112 
113       Saves each timestep into a seperate file named SA-data/SA-%06d.bin. The file name can be changed.
114 
115       This version saves the solutions at all the stages
116 
117       $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
118 
119   Level: intermediate
120 
121 .seealso:  TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile()
122 
123 M*/
124 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
125 {
126   PetscFunctionBegin;
127   tj->keepfiles    = PETSC_FALSE;
128   tj->ops->set     = TSTrajectorySet_Basic;
129   tj->ops->get     = TSTrajectoryGet_Basic;
130   tj->ops->destroy = TSTrajectoryDestroy_Basic;
131   PetscFunctionReturn(0);
132 }
133