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++) PetscCall(MatView(S[i], tjbasic->viewer)); 45 } 46 } 47 PetscFunctionReturn(PETSC_SUCCESS); 48 } 49 50 static PetscErrorCode TSTrajectorySetFromOptions_Basic(TSTrajectory tj, PetscOptionItems *PetscOptionsObject) 51 { 52 PetscFunctionBegin; 53 PetscOptionsHeadBegin(PetscOptionsObject, "TS trajectory options for Basic type"); 54 PetscOptionsHeadEnd(); 55 PetscFunctionReturn(PETSC_SUCCESS); 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 Vec Sol; 63 PetscInt ns, i; 64 65 PetscFunctionBegin; 66 PetscCall(PetscSNPrintf(filename, sizeof(filename), tj->dirfiletemplate, stepnum)); 67 PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)tj), filename, FILE_MODE_READ, &viewer)); 68 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE)); 69 PetscCall(TSGetSolution(ts, &Sol)); 70 PetscCall(VecLoad(Sol, viewer)); 71 PetscCall(PetscViewerBinaryRead(viewer, t, 1, NULL, PETSC_REAL)); 72 if (stepnum && !tj->solution_only) { 73 Vec *Y; 74 PetscReal timepre; 75 PetscCall(TSGetStages(ts, &ns, &Y)); 76 for (i = 0; i < ns; i++) { 77 /* 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. */ 78 if (ts->stifflyaccurate && i == ns - 1) continue; 79 PetscCall(VecLoad(Y[i], viewer)); 80 } 81 PetscCall(PetscViewerBinaryRead(viewer, &timepre, 1, NULL, PETSC_REAL)); 82 if (tj->adjoint_solve_mode) PetscCall(TSSetTimeStep(ts, -(*t) + timepre)); 83 } 84 /* Tangent linear sensitivities needed by second-order adjoint */ 85 if (ts->forward_solve) { 86 if (!ts->stifflyaccurate) { 87 Mat A; 88 PetscCall(TSForwardGetSensitivities(ts, NULL, &A)); 89 PetscCall(MatLoad(A, viewer)); 90 } 91 if (stepnum) { 92 Mat *S; 93 PetscCall(TSForwardGetStages(ts, &ns, &S)); 94 for (i = 0; i < ns; i++) PetscCall(MatLoad(S[i], viewer)); 95 } 96 } 97 PetscCall(PetscViewerDestroy(&viewer)); 98 PetscFunctionReturn(PETSC_SUCCESS); 99 } 100 101 PetscErrorCode TSTrajectorySetUp_Basic(TSTrajectory tj, TS ts) 102 { 103 MPI_Comm comm; 104 PetscMPIInt rank; 105 106 PetscFunctionBegin; 107 PetscCall(PetscObjectGetComm((PetscObject)tj, &comm)); 108 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 109 if (rank == 0) { 110 char *dir = tj->dirname; 111 PetscBool flg; 112 113 if (!dir) { 114 char dtempname[16] = "TS-data-XXXXXX"; 115 PetscCall(PetscMkdtemp(dtempname)); 116 PetscCall(PetscStrallocpy(dtempname, &tj->dirname)); 117 } else { 118 PetscCall(PetscTestDirectory(dir, 'w', &flg)); 119 if (!flg) { 120 PetscCall(PetscTestFile(dir, 'r', &flg)); 121 PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_USER, "Specified path is a file - not a dir: %s", dir); 122 PetscCall(PetscMkdir(dir)); 123 } else SETERRQ(comm, PETSC_ERR_SUP, "Directory %s not empty", tj->dirname); 124 } 125 } 126 PetscCall(PetscBarrier((PetscObject)tj)); 127 PetscFunctionReturn(PETSC_SUCCESS); 128 } 129 130 static PetscErrorCode TSTrajectoryDestroy_Basic(TSTrajectory tj) 131 { 132 TSTrajectory_Basic *tjbasic = (TSTrajectory_Basic *)tj->data; 133 134 PetscFunctionBegin; 135 PetscCall(PetscViewerDestroy(&tjbasic->viewer)); 136 PetscCall(PetscFree(tjbasic)); 137 PetscFunctionReturn(PETSC_SUCCESS); 138 } 139 140 /*MC 141 TSTRAJECTORYBASIC - Stores each solution of the ODE/DAE in a file 142 143 Saves each timestep into a separate file named TS-data-XXXXXX/TS-%06d.bin. The file name can be changed. 144 145 This version saves the solutions at all the stages 146 147 $PETSC_DIR/share/petsc/matlab/PetscReadBinaryTrajectory.m can read in files created with this format 148 149 Level: intermediate 150 151 .seealso: [](ch_ts), `TSTrajectoryCreate()`, `TS`, `TSTrajectory`, `TSTrajectorySetType()`, `TSTrajectorySetDirname()`, `TSTrajectorySetFile()`, 152 `TSTrajectoryType` 153 M*/ 154 PETSC_EXTERN PetscErrorCode TSTrajectoryCreate_Basic(TSTrajectory tj, TS ts) 155 { 156 TSTrajectory_Basic *tjbasic; 157 158 PetscFunctionBegin; 159 PetscCall(PetscNew(&tjbasic)); 160 161 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)tj), &tjbasic->viewer)); 162 PetscCall(PetscViewerSetType(tjbasic->viewer, PETSCVIEWERBINARY)); 163 PetscCall(PetscViewerPushFormat(tjbasic->viewer, PETSC_VIEWER_NATIVE)); 164 PetscCall(PetscViewerFileSetMode(tjbasic->viewer, FILE_MODE_WRITE)); 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(PETSC_SUCCESS); 173 } 174