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