xref: /petsc/src/ts/trajectory/impls/basic/trajbasic.c (revision e611a964e9853b74d61a56642fe9d06a6e51780f)
1 
2 #include <petsc/private/tsimpl.h>        /*I "petscts.h"  I*/
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "OutputBIN"
6 static PetscErrorCode OutputBIN(const char *filename,PetscViewer *viewer)
7 {
8   PetscErrorCode ierr;
9 
10   PetscFunctionBegin;
11   ierr = PetscViewerCreate(PETSC_COMM_WORLD,viewer);CHKERRQ(ierr);
12   ierr = PetscViewerSetType(*viewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
13   ierr = PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
14   ierr = PetscViewerFileSetName(*viewer,filename);CHKERRQ(ierr);
15   PetscFunctionReturn(0);
16 }
17 
18 #undef __FUNCT__
19 #define __FUNCT__ "TSTrajectorySet_Basic"
20 static PetscErrorCode TSTrajectorySet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X)
21 {
22   PetscViewer    viewer;
23   PetscInt       ns,i;
24   Vec            *Y;
25   char           filename[PETSC_MAX_PATH_LEN];
26   PetscReal      tprev;
27   PetscErrorCode ierr;
28 
29   PetscFunctionBegin;
30   ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr);
31   if (stepnum == 0) {
32     PetscMPIInt rank;
33     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);CHKERRQ(ierr);
34     if (!rank) {
35       ierr = PetscRMTree("SA-data");CHKERRQ(ierr);
36       ierr = PetscMkdir("SA-data");CHKERRQ(ierr);
37     }
38     ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-%06d.bin",stepnum);CHKERRQ(ierr);
39     ierr = OutputBIN(filename,&viewer);CHKERRQ(ierr);
40     ierr = VecView(X,viewer);CHKERRQ(ierr);
41     ierr = PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
42     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
43     PetscFunctionReturn(0);
44   }
45   ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-%06d.bin",stepnum);CHKERRQ(ierr);
46   ierr = OutputBIN(filename,&viewer);CHKERRQ(ierr);
47   ierr = VecView(X,viewer);CHKERRQ(ierr);
48   ierr = PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
49 
50   ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr);
51   for (i=0;i<ns;i++) {
52     ierr = VecView(Y[i],viewer);CHKERRQ(ierr);
53   }
54 
55   ierr = TSGetPrevTime(ts,&tprev);CHKERRQ(ierr);
56   ierr = PetscViewerBinaryWrite(viewer,&tprev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
57 
58   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
59   PetscFunctionReturn(0);
60 }
61 
62 #undef __FUNCT__
63 #define __FUNCT__ "TSTrajectoryGet_Basic"
64 static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
65 {
66   Vec            Sol,*Y;
67   PetscInt       Nr,i;
68   PetscViewer    viewer;
69   PetscReal      timepre;
70   char           filename[PETSC_MAX_PATH_LEN];
71   PetscErrorCode ierr;
72 
73   PetscFunctionBegin;
74   ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr);
75   ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-%06d.bin",stepnum);CHKERRQ(ierr);
76   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
77 
78   ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr);
79   ierr = VecLoad(Sol,viewer);CHKERRQ(ierr);
80 
81   ierr = PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);CHKERRQ(ierr);
82 
83   if (stepnum != 0) {
84     ierr = TSGetStages(ts,&Nr,&Y);CHKERRQ(ierr);
85     for (i=0;i<Nr ;i++) {
86       ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr);
87     }
88     ierr = PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);CHKERRQ(ierr);
89     ierr = TSSetTimeStep(ts,-(*t)+timepre);CHKERRQ(ierr);
90   }
91 
92   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
93   PetscFunctionReturn(0);
94 }
95 
96 /*MC
97       TSTRAJECTORYBASIC - Stores each solution of the ODE/ADE in a file
98 
99   Level: intermediate
100 
101 .seealso:  TSTrajectoryCreate(), TS, TSTrajectorySetType()
102 
103 M*/
104 #undef __FUNCT__
105 #define __FUNCT__ "TSTrajectoryCreate_Basic"
106 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts)
107 {
108   PetscFunctionBegin;
109   tj->ops->set  = TSTrajectorySet_Basic;
110   tj->ops->get  = TSTrajectoryGet_Basic;
111   PetscFunctionReturn(0);
112 }
113