1 #include <petsc/private/tshistoryimpl.h> 2 3 /* These macros can be moved to petscimpl.h eventually */ 4 #if defined(PETSC_USE_DEBUG) 5 6 #define PetscValidLogicalCollectiveIntComm(a, b, c) \ 7 do { \ 8 PetscInt b1[2], b2[2]; \ 9 b1[0] = -b; \ 10 b1[1] = b; \ 11 PetscCall(MPIU_Allreduce(b1, b2, 2, MPIU_INT, MPI_MAX, a)); \ 12 PetscCheck(-b2[0] == b2[1], a, PETSC_ERR_ARG_WRONG, "Int value must be same on all processes, argument # %d", c); \ 13 } while (0) 14 15 #define PetscValidLogicalCollectiveBoolComm(a, b, c) \ 16 do { \ 17 PetscMPIInt b1[2], b2[2]; \ 18 b1[0] = -(PetscMPIInt)b; \ 19 b1[1] = (PetscMPIInt)b; \ 20 PetscCall(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, a)); \ 21 PetscCheck(-b2[0] == b2[1], a, PETSC_ERR_ARG_WRONG, "Bool value must be same on all processes, argument # %d", c); \ 22 } while (0) 23 24 #define PetscValidLogicalCollectiveRealComm(a, b, c) \ 25 do { \ 26 PetscReal b1[3], b2[3]; \ 27 if (PetscIsNanReal(b)) { \ 28 b1[2] = 1; \ 29 } else { \ 30 b1[2] = 0; \ 31 }; \ 32 b1[0] = -b; \ 33 b1[1] = b; \ 34 PetscCallMPI(MPI_Allreduce(b1, b2, 3, MPIU_REAL, MPIU_MAX, a)); \ 35 PetscCheck((b2[2] == 1) || PetscEqualReal(-b2[0], b2[1]), a, PETSC_ERR_ARG_WRONG, "Real value must be same on all processes, argument # %d", c); \ 36 } while (0) 37 38 #else 39 40 #define PetscValidLogicalCollectiveRealComm(a, b, c) \ 41 do { \ 42 } while (0) 43 #define PetscValidLogicalCollectiveIntComm(a, b, c) \ 44 do { \ 45 } while (0) 46 #define PetscValidLogicalCollectiveBoolComm(a, b, c) \ 47 do { \ 48 } while (0) 49 50 #endif 51 52 struct _n_TSHistory { 53 MPI_Comm comm; /* used for runtime collective checks */ 54 PetscReal *hist; /* time history */ 55 PetscInt *hist_id; /* stores the stepid in time history */ 56 size_t n; /* current number of steps registered */ 57 PetscBool sorted; /* if the history is sorted in ascending order */ 58 size_t c; /* current capacity of history */ 59 size_t s; /* reallocation size */ 60 }; 61 62 PetscErrorCode TSHistoryGetNumSteps(TSHistory tsh, PetscInt *n) 63 { 64 PetscFunctionBegin; 65 PetscValidIntPointer(n, 2); 66 *n = tsh->n; 67 PetscFunctionReturn(PETSC_SUCCESS); 68 } 69 70 PetscErrorCode TSHistoryUpdate(TSHistory tsh, PetscInt id, PetscReal time) 71 { 72 PetscFunctionBegin; 73 if (tsh->n == tsh->c) { /* reallocation */ 74 tsh->c += tsh->s; 75 PetscCall(PetscRealloc(tsh->c * sizeof(*tsh->hist), &tsh->hist)); 76 PetscCall(PetscRealloc(tsh->c * sizeof(*tsh->hist_id), &tsh->hist_id)); 77 } 78 tsh->sorted = (PetscBool)(tsh->sorted && (tsh->n ? time >= tsh->hist[tsh->n - 1] : PETSC_TRUE)); 79 #if defined(PETSC_USE_DEBUG) 80 if (tsh->n) { /* id should be unique */ 81 PetscInt loc, *ids; 82 83 PetscCall(PetscMalloc1(tsh->n, &ids)); 84 PetscCall(PetscArraycpy(ids, tsh->hist_id, tsh->n)); 85 PetscCall(PetscSortInt(tsh->n, ids)); 86 PetscCall(PetscFindInt(id, tsh->n, ids, &loc)); 87 PetscCall(PetscFree(ids)); 88 PetscCheck(loc < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "History id should be unique"); 89 } 90 #endif 91 tsh->hist[tsh->n] = time; 92 tsh->hist_id[tsh->n] = id; 93 tsh->n += 1; 94 PetscFunctionReturn(PETSC_SUCCESS); 95 } 96 97 PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t) 98 { 99 PetscFunctionBegin; 100 if (!t) PetscFunctionReturn(PETSC_SUCCESS); 101 PetscValidRealPointer(t, 4); 102 if (!tsh->sorted) { 103 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id)); 104 tsh->sorted = PETSC_TRUE; 105 } 106 PetscCheck(step >= 0 && step < (PetscInt)tsh->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Given time step %" PetscInt_FMT " does not match any in history [0,%" PetscInt_FMT "]", step, (PetscInt)tsh->n); 107 if (!backward) *t = tsh->hist[step]; 108 else *t = tsh->hist[tsh->n - step - 1]; 109 PetscFunctionReturn(PETSC_SUCCESS); 110 } 111 112 PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt) 113 { 114 PetscFunctionBegin; 115 if (!dt) PetscFunctionReturn(PETSC_SUCCESS); 116 PetscValidRealPointer(dt, 4); 117 if (!tsh->sorted) { 118 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id)); 119 tsh->sorted = PETSC_TRUE; 120 } 121 PetscCheck(step >= 0 && step <= (PetscInt)tsh->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Given time step %" PetscInt_FMT " does not match any in history [0,%" PetscInt_FMT "]", step, (PetscInt)tsh->n); 122 if (!backward) *dt = tsh->hist[PetscMin(step + 1, (PetscInt)tsh->n - 1)] - tsh->hist[PetscMin(step, (PetscInt)tsh->n - 1)]; 123 else *dt = tsh->hist[PetscMax((PetscInt)tsh->n - step - 1, 0)] - tsh->hist[PetscMax((PetscInt)tsh->n - step - 2, 0)]; 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc) 128 { 129 PetscFunctionBegin; 130 PetscValidIntPointer(loc, 3); 131 if (!tsh->sorted) { 132 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id)); 133 tsh->sorted = PETSC_TRUE; 134 } 135 PetscCall(PetscFindReal(time, tsh->n, tsh->hist, PETSC_SMALL, loc)); 136 PetscFunctionReturn(PETSC_SUCCESS); 137 } 138 139 PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted) 140 { 141 PetscFunctionBegin; 142 PetscValidLogicalCollectiveIntComm(tsh->comm, n, 2); 143 PetscCheck(n >= 0, tsh->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot request a negative size for history storage"); 144 if (n) PetscValidRealPointer(hist, 3); 145 PetscCall(PetscFree(tsh->hist)); 146 PetscCall(PetscFree(tsh->hist_id)); 147 tsh->n = (size_t)n; 148 tsh->c = (size_t)n; 149 PetscCall(PetscMalloc1(tsh->n, &tsh->hist)); 150 PetscCall(PetscMalloc1(tsh->n, &tsh->hist_id)); 151 for (PetscInt i = 0; i < (PetscInt)tsh->n; i++) { 152 tsh->hist[i] = hist[i]; 153 tsh->hist_id[i] = hist_id ? hist_id[i] : i; 154 } 155 if (!sorted) PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id)); 156 tsh->sorted = PETSC_TRUE; 157 PetscFunctionReturn(PETSC_SUCCESS); 158 } 159 160 PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal *hist[], const PetscInt *hist_id[], PetscBool *sorted) 161 { 162 PetscFunctionBegin; 163 if (n) *n = tsh->n; 164 if (hist) *hist = tsh->hist; 165 if (hist_id) *hist_id = tsh->hist_id; 166 if (sorted) *sorted = tsh->sorted; 167 PetscFunctionReturn(PETSC_SUCCESS); 168 } 169 170 PetscErrorCode TSHistoryDestroy(TSHistory *tsh) 171 { 172 PetscFunctionBegin; 173 if (!*tsh) PetscFunctionReturn(PETSC_SUCCESS); 174 PetscCall(PetscFree((*tsh)->hist)); 175 PetscCall(PetscFree((*tsh)->hist_id)); 176 PetscCall(PetscCommDestroy(&((*tsh)->comm))); 177 PetscCall(PetscFree((*tsh))); 178 *tsh = NULL; 179 PetscFunctionReturn(PETSC_SUCCESS); 180 } 181 182 PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst) 183 { 184 TSHistory tsh; 185 186 PetscFunctionBegin; 187 PetscValidPointer(hst, 2); 188 *hst = NULL; 189 PetscCall(PetscNew(&tsh)); 190 PetscCall(PetscCommDuplicate(comm, &tsh->comm, NULL)); 191 192 tsh->c = 1024; /* capacity */ 193 tsh->s = 1024; /* reallocation size */ 194 tsh->sorted = PETSC_TRUE; 195 196 PetscCall(PetscMalloc1(tsh->c, &tsh->hist)); 197 PetscCall(PetscMalloc1(tsh->c, &tsh->hist_id)); 198 *hst = tsh; 199 PetscFunctionReturn(PETSC_SUCCESS); 200 } 201