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 = TSGetStepNumber(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(tj->dirname);CHKERRQ(ierr); 34 ierr = PetscMkdir(tj->dirname);CHKERRQ(ierr); 35 } 36 ierr = PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,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),tj->dirfiletemplate,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 = PetscSNPrintf(filename,sizeof filename,tj->dirfiletemplate,stepnum);CHKERRQ(ierr); 71 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 72 73 ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr); 74 ierr = VecLoad(Sol,viewer);CHKERRQ(ierr); 75 76 ierr = PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);CHKERRQ(ierr); 77 78 if (stepnum != 0) { 79 ierr = TSGetStages(ts,&Nr,&Y);CHKERRQ(ierr); 80 for (i=0;i<Nr ;i++) { 81 ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr); 82 } 83 ierr = PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);CHKERRQ(ierr); 84 ierr = TSSetTimeStep(ts,-(*t)+timepre);CHKERRQ(ierr); 85 } 86 87 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 88 PetscFunctionReturn(0); 89 } 90 91 static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj) 92 { 93 PetscErrorCode ierr; 94 PetscMPIInt rank; 95 MPI_Comm comm; 96 97 PetscFunctionBegin; 98 if (!tj->keepfiles) { 99 ierr = PetscObjectGetComm((PetscObject)tj,&comm);CHKERRQ(ierr); 100 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 101 if (!rank) { 102 ierr = PetscRMTree(tj->dirname);CHKERRQ(ierr); 103 } 104 } 105 PetscFunctionReturn(0); 106 } 107 108 /*MC 109 TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file 110 111 Saves each timestep into a seperate file named SA-data/SA-%06d.bin. The file name can be changed. 112 113 This version saves the solutions at all the stages 114 115 $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format 116 117 Level: intermediate 118 119 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile() 120 121 M*/ 122 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts) 123 { 124 PetscFunctionBegin; 125 tj->keepfiles = PETSC_FALSE; 126 tj->ops->set = TSTrajectorySet_Basic; 127 tj->ops->get = TSTrajectoryGet_Basic; 128 tj->ops->destroy = TSTrajectoryDestroy_Basic; 129 PetscFunctionReturn(0); 130 } 131