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 PetscFunctionBegin; 64 PetscValidIntPointer(n, 2); 65 *n = tsh->n; 66 PetscFunctionReturn(0); 67 } 68 69 PetscErrorCode TSHistoryUpdate(TSHistory tsh, PetscInt id, PetscReal time) { 70 PetscFunctionBegin; 71 if (tsh->n == tsh->c) { /* reallocation */ 72 tsh->c += tsh->s; 73 PetscCall(PetscRealloc(tsh->c * sizeof(*tsh->hist), &tsh->hist)); 74 PetscCall(PetscRealloc(tsh->c * sizeof(*tsh->hist_id), &tsh->hist_id)); 75 } 76 tsh->sorted = (PetscBool)(tsh->sorted && (tsh->n ? time >= tsh->hist[tsh->n - 1] : PETSC_TRUE)); 77 #if defined(PETSC_USE_DEBUG) 78 if (tsh->n) { /* id should be unique */ 79 PetscInt loc, *ids; 80 81 PetscCall(PetscMalloc1(tsh->n, &ids)); 82 PetscCall(PetscArraycpy(ids, tsh->hist_id, tsh->n)); 83 PetscCall(PetscSortInt(tsh->n, ids)); 84 PetscCall(PetscFindInt(id, tsh->n, ids, &loc)); 85 PetscCall(PetscFree(ids)); 86 PetscCheck(loc < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "History id should be unique"); 87 } 88 #endif 89 tsh->hist[tsh->n] = time; 90 tsh->hist_id[tsh->n] = id; 91 tsh->n += 1; 92 PetscFunctionReturn(0); 93 } 94 95 PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t) { 96 PetscFunctionBegin; 97 if (!t) PetscFunctionReturn(0); 98 PetscValidRealPointer(t, 4); 99 if (!tsh->sorted) { 100 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id)); 101 tsh->sorted = PETSC_TRUE; 102 } 103 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); 104 if (!backward) *t = tsh->hist[step]; 105 else *t = tsh->hist[tsh->n - step - 1]; 106 PetscFunctionReturn(0); 107 } 108 109 PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt) { 110 PetscFunctionBegin; 111 if (!dt) PetscFunctionReturn(0); 112 PetscValidRealPointer(dt, 4); 113 if (!tsh->sorted) { 114 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id)); 115 tsh->sorted = PETSC_TRUE; 116 } 117 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); 118 if (!backward) *dt = tsh->hist[PetscMin(step + 1, (PetscInt)tsh->n - 1)] - tsh->hist[PetscMin(step, (PetscInt)tsh->n - 1)]; 119 else *dt = tsh->hist[PetscMax((PetscInt)tsh->n - step - 1, 0)] - tsh->hist[PetscMax((PetscInt)tsh->n - step - 2, 0)]; 120 PetscFunctionReturn(0); 121 } 122 123 PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc) { 124 PetscFunctionBegin; 125 PetscValidIntPointer(loc, 3); 126 if (!tsh->sorted) { 127 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id)); 128 tsh->sorted = PETSC_TRUE; 129 } 130 PetscCall(PetscFindReal(time, tsh->n, tsh->hist, PETSC_SMALL, loc)); 131 PetscFunctionReturn(0); 132 } 133 134 PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted) { 135 PetscFunctionBegin; 136 PetscValidLogicalCollectiveIntComm(tsh->comm, n, 2); 137 PetscCheck(n >= 0, tsh->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot request a negative size for history storage"); 138 if (n) PetscValidRealPointer(hist, 3); 139 PetscCall(PetscFree(tsh->hist)); 140 PetscCall(PetscFree(tsh->hist_id)); 141 tsh->n = (size_t)n; 142 tsh->c = (size_t)n; 143 PetscCall(PetscMalloc1(tsh->n, &tsh->hist)); 144 PetscCall(PetscMalloc1(tsh->n, &tsh->hist_id)); 145 for (PetscInt i = 0; i < (PetscInt)tsh->n; i++) { 146 tsh->hist[i] = hist[i]; 147 tsh->hist_id[i] = hist_id ? hist_id[i] : i; 148 } 149 if (!sorted) PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id)); 150 tsh->sorted = PETSC_TRUE; 151 PetscFunctionReturn(0); 152 } 153 154 PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal *hist[], const PetscInt *hist_id[], PetscBool *sorted) { 155 PetscFunctionBegin; 156 if (n) *n = tsh->n; 157 if (hist) *hist = tsh->hist; 158 if (hist_id) *hist_id = tsh->hist_id; 159 if (sorted) *sorted = tsh->sorted; 160 PetscFunctionReturn(0); 161 } 162 163 PetscErrorCode TSHistoryDestroy(TSHistory *tsh) { 164 PetscFunctionBegin; 165 if (!*tsh) PetscFunctionReturn(0); 166 PetscCall(PetscFree((*tsh)->hist)); 167 PetscCall(PetscFree((*tsh)->hist_id)); 168 PetscCall(PetscCommDestroy(&((*tsh)->comm))); 169 PetscCall(PetscFree((*tsh))); 170 *tsh = NULL; 171 PetscFunctionReturn(0); 172 } 173 174 PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst) { 175 TSHistory tsh; 176 177 PetscFunctionBegin; 178 PetscValidPointer(hst, 2); 179 *hst = NULL; 180 PetscCall(PetscNew(&tsh)); 181 PetscCall(PetscCommDuplicate(comm, &tsh->comm, NULL)); 182 183 tsh->c = 1024; /* capacity */ 184 tsh->s = 1024; /* reallocation size */ 185 tsh->sorted = PETSC_TRUE; 186 187 PetscCall(PetscMalloc1(tsh->c, &tsh->hist)); 188 PetscCall(PetscMalloc1(tsh->c, &tsh->hist_id)); 189 *hst = tsh; 190 PetscFunctionReturn(0); 191 } 192