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 PetscCall(PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum)); 20 PetscCall(PetscViewerFileSetName(tjbasic->viewer,filename)); /* this triggers PetscViewer to be set up again */ 21 PetscCall(PetscViewerSetUp(tjbasic->viewer)); 22 PetscCall(VecView(X,tjbasic->viewer)); 23 PetscCall(PetscViewerBinaryWrite(tjbasic->viewer,&time,1,PETSC_REAL)); 24 if (stepnum && !tj->solution_only) { 25 Vec *Y; 26 PetscReal tprev; 27 PetscCall(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 PetscCall(VecView(Y[i],tjbasic->viewer)); 32 } 33 PetscCall(TSGetPrevTime(ts,&tprev)); 34 PetscCall(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 PetscCall(TSForwardGetSensitivities(ts,NULL,&A)); 41 PetscCall(MatView(A,tjbasic->viewer)); 42 if (stepnum) { 43 PetscCall(TSForwardGetStages(ts,&ns,&S)); 44 for (i=0; i<ns; i++) { 45 PetscCall(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 PetscCall(PetscOptionsHead(PetscOptionsObject,"TS trajectory options for Basic type")); 56 PetscCall(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 PetscCall(PetscSNPrintf(filename,sizeof(filename),tj->dirfiletemplate,stepnum)); 69 PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj),filename,FILE_MODE_READ,&viewer)); 70 PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE)); 71 PetscCall(TSGetSolution(ts,&Sol)); 72 PetscCall(VecLoad(Sol,viewer)); 73 PetscCall(PetscViewerBinaryRead(viewer,t,1,NULL,PETSC_REAL)); 74 if (stepnum && !tj->solution_only) { 75 Vec *Y; 76 PetscReal timepre; 77 PetscCall(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 PetscCall(VecLoad(Y[i],viewer)); 82 } 83 PetscCall(PetscViewerBinaryRead(viewer,&timepre,1,NULL,PETSC_REAL)); 84 if (tj->adjoint_solve_mode) { 85 PetscCall(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 PetscCall(TSForwardGetSensitivities(ts,NULL,&A)); 94 PetscCall(MatLoad(A,viewer)); 95 } 96 if (stepnum) { 97 Mat *S; 98 PetscCall(TSForwardGetStages(ts,&ns,&S)); 99 for (i=0; i<ns; i++) { 100 PetscCall(MatLoad(S[i],viewer)); 101 } 102 } 103 } 104 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject)tj,&comm)); 115 PetscCallMPI(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 PetscCall(PetscMkdtemp(dtempname)); 123 PetscCall(PetscStrallocpy(dtempname,&tj->dirname)); 124 } else { 125 PetscCall(PetscTestDirectory(dir,'w',&flg)); 126 if (!flg) { 127 PetscCall(PetscTestFile(dir,'r',&flg)); 128 PetscCheck(!flg,PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir); 129 PetscCall(PetscMkdir(dir)); 130 } else SETERRQ(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname); 131 } 132 } 133 PetscCall(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 PetscCall(PetscViewerDestroy(&tjbasic->viewer)); 143 PetscCall(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 PetscCall(PetscNew(&tjbasic)); 167 168 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer)); 169 PetscCall(PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY)); 170 PetscCall(PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE)); 171 PetscCall(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