1 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 4 typedef struct { 5 PetscViewer viewer; 6 } TSTrajectory_Basic; 7 8 static PetscErrorCode TSTrajectorySet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal time,Vec X) 9 { 10 TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data; 11 char filename[PETSC_MAX_PATH_LEN]; 12 PetscInt ns,i; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 ierr = PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);CHKERRQ(ierr); 17 ierr = PetscViewerFileSetName(tjbasic->viewer,filename);CHKERRQ(ierr); /* this triggers PetscViewer to be set up again */ 18 ierr = PetscViewerSetUp(tjbasic->viewer);CHKERRQ(ierr); 19 ierr = VecView(X,tjbasic->viewer);CHKERRQ(ierr); 20 ierr = PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 21 if (stepnum && !tj->solution_only) { 22 Vec *Y; 23 PetscReal tprev; 24 25 ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr); 26 for (i=0; i<ns; i++) { 27 ierr = VecView(Y[i],tjbasic->viewer);CHKERRQ(ierr); 28 } 29 ierr = TSGetPrevTime(ts,&tprev);CHKERRQ(ierr); 30 ierr = PetscViewerBinaryWrite(tjbasic->viewer,&tprev,1,PETSC_REAL,PETSC_FALSE);CHKERRQ(ierr); 31 } 32 /* Tangent linear sensitivities needed by second-order adjoint */ 33 if (ts->forward_solve) { 34 Mat A,*S; 35 36 ierr = TSForwardGetSensitivities(ts,NULL,&A);CHKERRQ(ierr); 37 ierr = MatView(A,tjbasic->viewer);CHKERRQ(ierr); 38 if (stepnum) { 39 ierr = TSForwardGetStages(ts,&ns,&S);CHKERRQ(ierr); 40 for (i=0; i<ns; i++) { 41 ierr = MatView(S[i],tjbasic->viewer);CHKERRQ(ierr); 42 } 43 } 44 } 45 PetscFunctionReturn(0); 46 } 47 48 static PetscErrorCode TSTrajectorySetFromOptions_Basic(PetscOptionItems *PetscOptionsObject,TSTrajectory tj) 49 { 50 PetscErrorCode ierr; 51 52 PetscFunctionBegin; 53 ierr = PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type");CHKERRQ(ierr); 54 ierr = PetscOptionsTail();CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode TSTrajectoryGet_Basic(TSTrajectory tj,TS ts,PetscInt stepnum,PetscReal *t) 59 { 60 PetscViewer viewer; 61 char filename[PETSC_MAX_PATH_LEN]; 62 PetscErrorCode ierr; 63 Vec Sol; 64 PetscInt ns,i; 65 66 PetscFunctionBegin; 67 ierr = PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum);CHKERRQ(ierr); 68 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 69 ierr = TSGetSolution(ts,&Sol);CHKERRQ(ierr); 70 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);CHKERRQ(ierr); 71 ierr = VecLoad(Sol,viewer);CHKERRQ(ierr); 72 ierr = PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL);CHKERRQ(ierr); 73 if (stepnum && !tj->solution_only) { 74 Vec *Y; 75 PetscReal timepre; 76 77 ierr = TSGetStages(ts,&ns,&Y);CHKERRQ(ierr); 78 for (i=0; i<ns; i++) { 79 ierr = VecLoad(Y[i],viewer);CHKERRQ(ierr); 80 } 81 ierr = PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL);CHKERRQ(ierr); 82 if (tj->adjoint_solve_mode) { 83 ierr = TSSetTimeStep(ts,-(*t)+timepre);CHKERRQ(ierr); 84 } 85 } 86 /* Tangent linear sensitivities needed by second-order adjoint */ 87 if (ts->forward_solve) { 88 Mat A,*S; 89 90 ierr = TSForwardGetSensitivities(ts,NULL,&A);CHKERRQ(ierr); 91 ierr = MatLoad(A,viewer);CHKERRQ(ierr); 92 if (stepnum) { 93 ierr = TSForwardGetStages(ts,&ns,&S);CHKERRQ(ierr); 94 for (i=0; i<ns; i++) { 95 ierr = MatLoad(S[i],viewer);CHKERRQ(ierr); 96 } 97 } 98 } 99 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 103 static PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts) 104 { 105 MPI_Comm comm; 106 PetscMPIInt rank; 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 ierr = PetscObjectGetComm((PetscObject)tj,&comm);CHKERRQ(ierr); 111 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 112 if (!rank) { 113 const char* dir = tj->dirname; 114 PetscBool flg; 115 116 /* I don't like running PetscRMTree on a directory */ 117 ierr = PetscTestDirectory(dir,'w',&flg);CHKERRQ(ierr); 118 if (!flg) { 119 ierr = PetscTestFile(dir,'r',&flg);CHKERRQ(ierr); 120 if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir); 121 ierr = PetscMkdir(dir);CHKERRQ(ierr); 122 } else SETERRQ1(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname); 123 } 124 ierr = PetscBarrier((PetscObject)tj);CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127 128 static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj) 129 { 130 TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data; 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 ierr = PetscViewerDestroy(&tjbasic->viewer);CHKERRQ(ierr); 135 ierr = PetscFree(tjbasic);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 } 138 139 /*MC 140 TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file 141 142 Saves each timestep into a seperate file named SA-data/SA-%06d.bin. The file name can be changed. 143 144 This version saves the solutions at all the stages 145 146 $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format 147 148 Level: intermediate 149 150 .seealso: TSTrajectoryCreate(), TS, TSTrajectorySetType(), TSTrajectorySetDirname(), TSTrajectorySetFile() 151 152 M*/ 153 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts) 154 { 155 TSTrajectory_Basic *tjbasic; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ierr = PetscNew(&tjbasic);CHKERRQ(ierr); 160 161 ierr = PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer);CHKERRQ(ierr); 162 ierr = PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY);CHKERRQ(ierr); 163 ierr = PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE);CHKERRQ(ierr); 164 ierr = PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE);CHKERRQ(ierr); 165 tj->data = tjbasic; 166 167 tj->ops->set = TSTrajectorySet_Basic; 168 tj->ops->get = TSTrajectoryGet_Basic; 169 tj->ops->setup = TSTrajectorySetUp_Basic; 170 tj->ops->destroy = TSTrajectoryDestroy_Basic; 171 tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic; 172 PetscFunctionReturn(0); 173 } 174