Actual source code: tshistory.c
1: #include <petsc/private/tshistoryimpl.h>
3: /* These macros can be moved to petscimpl.h eventually */
4: #if defined(PETSC_USE_DEBUG)
7: do { \
8: PetscInt b1[2], b2[2]; \
9: b1[0] = -b; \
10: b1[1] = b; \
11: PetscCallMPI(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)
16: do { \
17: PetscMPIInt b1[2], b2[2]; \
18: b1[0] = -(PetscMPIInt)b; \
19: b1[1] = (PetscMPIInt)b; \
20: PetscCallMPI(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)
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(MPIU_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)
38: #else
41: do { \
42: } while (0)
44: do { \
45: } while (0)
47: do { \
48: } while (0)
50: #endif
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: PetscCount n; /* current number of steps registered */
57: PetscBool sorted; /* if the history is sorted in ascending order */
58: PetscCount c; /* current capacity of history */
59: PetscCount s; /* reallocation size */
60: };
62: PetscErrorCode TSHistoryGetNumSteps(TSHistory tsh, PetscInt *n)
63: {
64: PetscFunctionBegin;
65: PetscAssertPointer(n, 2);
66: PetscCall(PetscIntCast(tsh->n, n));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
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 ? (PetscBool)(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;
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: }
97: PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t)
98: {
99: PetscFunctionBegin;
100: if (!t) PetscFunctionReturn(PETSC_SUCCESS);
101: PetscAssertPointer(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: }
112: PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt)
113: {
114: PetscFunctionBegin;
115: if (!dt) PetscFunctionReturn(PETSC_SUCCESS);
116: PetscAssertPointer(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: }
127: PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc)
128: {
129: PetscFunctionBegin;
130: PetscAssertPointer(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: }
139: PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted)
140: {
141: PetscFunctionBegin;
143: PetscCheck(n >= 0, tsh->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot request a negative size for history storage");
144: if (n) PetscAssertPointer(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: }
160: PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal *hist[], const PetscInt *hist_id[], PetscBool *sorted)
161: {
162: PetscFunctionBegin;
163: if (n) PetscCall(PetscIntCast(tsh->n, 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: }
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: }
182: PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst)
183: {
184: TSHistory tsh;
186: PetscFunctionBegin;
187: PetscAssertPointer(hst, 2);
188: PetscCall(PetscNew(&tsh));
189: PetscCall(PetscCommDuplicate(comm, &tsh->comm, NULL));
191: tsh->c = 1024; /* capacity */
192: tsh->s = 1024; /* reallocation size */
193: tsh->sorted = PETSC_TRUE;
195: PetscCall(PetscMalloc1(tsh->c, &tsh->hist));
196: PetscCall(PetscMalloc1(tsh->c, &tsh->hist_id));
197: *hst = tsh;
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }