xref: /petsc/src/ts/trajectory/impls/basic/trajbasic.c (revision 07c0461f92f2c1aea27ba8cc5bb53c8785e502d5)
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];
22   PetscReal      tprev;
23   PetscErrorCode ierr;
24   MPI_Comm       comm;
25 
26   PetscFunctionBegin;
27   ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr);
28   ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr);
29   if (stepnum == 0) {
30     PetscMPIInt rank;
31     ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
32     if (!rank) {
33       ierr = PetscRMTree("SA-data");CHKERRQ(ierr);
34       ierr = PetscMkdir("SA-data");CHKERRQ(ierr);
35     }
36     ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-%06d.bin",stepnum);CHKERRQ(ierr);
37     ierr = OutputBIN(comm,filename,&viewer);CHKERRQ(ierr);
38     ierr = VecView(X,viewer);CHKERRQ(ierr);
39     ierr = PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
40     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
41     PetscFunctionReturn(0);
42   }
43   ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-%06d.bin",stepnum);CHKERRQ(ierr);
44   ierr = OutputBIN(comm,filename,&viewer);CHKERRQ(ierr);
45   ierr = VecView(X,viewer);CHKERRQ(ierr);
46   ierr = PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
47 
48   ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr);
49   for (i=0;i<ns;i++) {
50     ierr = VecView(Y[i],viewer);CHKERRQ(ierr);
51   }
52 
53   ierr = TSGetPrevTime(ts,&tprev);CHKERRQ(ierr);
54   ierr = PetscViewerBinaryWrite(viewer,&tprev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr);
55 
56   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
57   PetscFunctionReturn(0);
58 }
59 
60 static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t)
61 {
62   Vec            Sol,*Y;
63   PetscInt       Nr,i;
64   PetscViewer    viewer;
65   PetscReal      timepre;
66   char           filename[PETSC_MAX_PATH_LEN];
67   PetscErrorCode ierr;
68 
69   PetscFunctionBegin;
70   ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr);
71   ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-%06d.bin",stepnum);CHKERRQ(ierr);
72   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
73 
74   ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr);
75   ierr = VecLoad(Sol,viewer);CHKERRQ(ierr);
76 
77   ierr = PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);CHKERRQ(ierr);
78 
79   if (stepnum != 0) {
80     ierr = TSGetStages(ts,&Nr,&Y);CHKERRQ(ierr);
81     for (i=0;i<Nr ;i++) {
82       ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr);
83     }
84     ierr = PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);CHKERRQ(ierr);
85     ierr = TSSetTimeStep(ts,-(*t)+timepre);CHKERRQ(ierr);
86   }
87 
88   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
89   PetscFunctionReturn(0);
90 }
91 
92 /*MC
93       TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file
94 
95       Saves each timestep into a seperate file in SA-data/SA-%06d.bin
96 
97       This version saves the solutions at all the stages
98 
99       $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format
100 
101   Level: intermediate
102 
103 .seealso:  TSTrajectoryCreate(), TS, TSTrajectorySetType()
104 
105 M*/
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