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 PetscOptionsHeadBegin(PetscOptionsObject,"TS trajectory options for Basic type"); 56 PetscOptionsHeadEnd(); 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) PetscCall(TSSetTimeStep(ts,-(*t)+timepre)); 85 } 86 /* Tangent linear sensitivities needed by second-order adjoint */ 87 if (ts->forward_solve) { 88 89 if (!ts->stifflyaccurate) { 90 Mat A; 91 PetscCall(TSForwardGetSensitivities(ts,NULL,&A)); 92 PetscCall(MatLoad(A,viewer)); 93 } 94 if (stepnum) { 95 Mat *S; 96 PetscCall(TSForwardGetStages(ts,&ns,&S)); 97 for (i=0; i<ns; i++) { 98 PetscCall(MatLoad(S[i],viewer)); 99 } 100 } 101 } 102 PetscCall(PetscViewerDestroy(&viewer)); 103 PetscFunctionReturn(0); 104 } 105 106 PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj,TS ts) 107 { 108 MPI_Comm comm; 109 PetscMPIInt rank; 110 111 PetscFunctionBegin; 112 PetscCall(PetscObjectGetComm((PetscObject)tj,&comm)); 113 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 114 if (rank == 0) { 115 char *dir = tj->dirname; 116 PetscBool flg; 117 118 if (!dir) { 119 char dtempname[16] = "TS-data-XXXXXX"; 120 PetscCall(PetscMkdtemp(dtempname)); 121 PetscCall(PetscStrallocpy(dtempname,&tj->dirname)); 122 } else { 123 PetscCall(PetscTestDirectory(dir,'w',&flg)); 124 if (!flg) { 125 PetscCall(PetscTestFile(dir,'r',&flg)); 126 PetscCheck(!flg,PETSC_COMM_SELF,PETSC_ERR_USER,"Specified path is a file - not a dir: %s",dir); 127 PetscCall(PetscMkdir(dir)); 128 } else SETERRQ(comm,PETSC_ERR_SUP,"Directory %s not empty",tj->dirname); 129 } 130 } 131 PetscCall(PetscBarrier((PetscObject)tj)); 132 PetscFunctionReturn(0); 133 } 134 135 static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj) 136 { 137 TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic*)tj->data; 138 139 PetscFunctionBegin; 140 PetscCall(PetscViewerDestroy(&tjbasic->viewer)); 141 PetscCall(PetscFree(tjbasic)); 142 PetscFunctionReturn(0); 143 } 144 145 /*MC 146 TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file 147 148 Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed. 149 150 This version saves the solutions at all the stages 151 152 $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format 153 154 Level: intermediate 155 156 .seealso: `TSTrajectoryCreate()`, `TS`, `TSTrajectorySetType()`, `TSTrajectorySetDirname()`, `TSTrajectorySetFile()` 157 158 M*/ 159 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj,TS ts) 160 { 161 TSTrajectory_Basic *tjbasic; 162 163 PetscFunctionBegin; 164 PetscCall(PetscNew(&tjbasic)); 165 166 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)tj),&tjbasic->viewer)); 167 PetscCall(PetscViewerSetType(tjbasic->viewer,PETSCVIEWERBINARY)); 168 PetscCall(PetscViewerPushFormat(tjbasic->viewer,PETSC_VIEWER_NATIVE)); 169 PetscCall(PetscViewerFileSetMode(tjbasic->viewer,FILE_MODE_WRITE)); 170 tj->data = tjbasic; 171 172 tj->ops->set = TSTrajectorySet_Basic; 173 tj->ops->get = TSTrajectoryGet_Basic; 174 tj->ops->setup = TSTrajectorySetUp_Basic; 175 tj->ops->destroy = TSTrajectoryDestroy_Basic; 176 tj->ops->setfromoptions = TSTrajectorySetFromOptions_Basic; 177 PetscFunctionReturn(0); 178 } 179