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