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 PetscFunctionBeginUser; 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 19 #undef __FUNCT__ 20 #define __FUNCT__ "TSTrajectorySet_Basic" 21 PetscErrorCode TSTrajectorySet_Basic(TSTrajectory jac,TS ts,PetscInt stepnum,PetscReal time,Vec X) 22 { 23 PetscViewer viewer; 24 PetscInt ns,i; 25 Vec *Y; 26 char filename[PETSC_MAX_PATH_LEN]; 27 PetscReal tprev; 28 PetscErrorCode ierr; 29 30 PetscFunctionBeginUser; 31 if (stepnum == 0) { 32 #if defined(PETSC_HAVE_POPEN) 33 ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr); 34 if (stepnum == 0) { 35 PetscMPIInt rank; 36 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);CHKERRQ(ierr); 37 if (!rank) { 38 char command[PETSC_MAX_PATH_LEN]; 39 FILE *fd; 40 int err; 41 42 ierr = PetscMemzero(command,sizeof(command));CHKERRQ(ierr); 43 ierr = PetscSNPrintf(command,PETSC_MAX_PATH_LEN,"rm -fr %s","SA-data");CHKERRQ(ierr); 44 ierr = PetscPOpen(PETSC_COMM_SELF,NULL,command,"r",&fd);CHKERRQ(ierr); 45 ierr = PetscPClose(PETSC_COMM_SELF,fd,&err);CHKERRQ(ierr); 46 ierr = PetscSNPrintf(command,PETSC_MAX_PATH_LEN,"mkdir %s","SA-data");CHKERRQ(ierr); 47 ierr = PetscPOpen(PETSC_COMM_SELF,NULL,command,"r",&fd);CHKERRQ(ierr); 48 ierr = PetscPClose(PETSC_COMM_SELF,fd,&err);CHKERRQ(ierr); 49 } 50 } 51 #endif 52 PetscFunctionReturn(0); 53 } 54 ierr = TSGetPrevTime(ts,&tprev);CHKERRQ(ierr); 55 ierr = TSGetTotalSteps(ts,&stepnum);CHKERRQ(ierr); 56 ierr = PetscSNPrintf(filename,sizeof(filename),"SA-data/SA-%06d.bin",stepnum);CHKERRQ(ierr); 57 ierr = OutputBIN(filename,&viewer);CHKERRQ(ierr); 58 ierr = VecView(X,viewer);CHKERRQ(ierr); 59 /* ierr = PetscRealView(1,&time,viewer);CHKERRQ(ierr); */ 60 ierr = PetscViewerBinaryWrite(viewer,&tprev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 61 /* ierr = PetscViewerBinaryWrite(viewer,&h ,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); */ 62 ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr); 63 64 for (i=0;i<ns;i++) { 65 ierr = VecView(Y[i],viewer);CHKERRQ(ierr); 66 } 67 68 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "TSTrajectoryGet_Basic" 74 PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory jac,TS ts,PetscInt step,PetscReal t) 75 { 76 PetscReal ptime; 77 Vec Sol,*Y; 78 PetscInt Nr,i; 79 PetscViewer viewer; 80 PetscReal timepre; 81 char filename[PETSC_MAX_PATH_LEN]; 82 PetscErrorCode ierr; 83 84 PetscFunctionBeginUser; 85 ierr = TSGetTotalSteps(ts,&step);CHKERRQ(ierr); 86 ierr = PetscSNPrintf(filename,sizeof filename,"SA-data/SA-%06d.bin",step);CHKERRQ(ierr); 87 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 88 89 ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr); 90 ierr = VecLoad(Sol,viewer);CHKERRQ(ierr); 91 92 Nr = 1; 93 /* ierr = PetscRealLoad(Nr,&Nr,&timepre,viewer);CHKERRQ(ierr); */ 94 ierr = PetscViewerBinaryRead(viewer,&timepre,1,PETSC_REAL);CHKERRQ(ierr); 95 96 ierr = TSGetStages(ts,&Nr,&Y);CHKERRQ(ierr); 97 for (i=0;i<Nr ;i++) { 98 ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr); 99 } 100 101 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 102 103 ierr = TSGetTime(ts,&ptime);CHKERRQ(ierr); 104 ierr = TSSetTimeStep(ts,-ptime+timepre);CHKERRQ(ierr); 105 106 PetscFunctionReturn(0); 107 } 108 109 /*MC 110 TSTRAJECTPRYBASIC - Stores each solution of the ODE/ADE in a file 111 112 Level: intermediate 113 114 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType() 115 116 M*/ 117 #undef __FUNCT__ 118 #define __FUNCT__ "TSTrajectoryCreate_Basic" 119 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory ts) 120 { 121 PetscFunctionBegin; 122 ts->ops->set = TSTrajectorySet_Basic; 123 ts->ops->get = TSTrajectoryGet_Basic; 124 PetscFunctionReturn(0); 125 } 126