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