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