1 #include <petsc/private/tshistoryimpl.h>
2
3 struct _n_TSHistory {
4 MPI_Comm comm; /* used for runtime collective checks */
5 PetscReal *hist; /* time history */
6 PetscInt *hist_id; /* stores the stepid in time history */
7 PetscCount n; /* current number of steps registered */
8 PetscBool sorted; /* if the history is sorted in ascending order */
9 PetscCount c; /* current capacity of history */
10 PetscCount s; /* reallocation size */
11 };
12
TSHistoryGetNumSteps(TSHistory tsh,PetscInt * n)13 PetscErrorCode TSHistoryGetNumSteps(TSHistory tsh, PetscInt *n)
14 {
15 PetscFunctionBegin;
16 PetscAssertPointer(n, 2);
17 PetscCall(PetscIntCast(tsh->n, n));
18 PetscFunctionReturn(PETSC_SUCCESS);
19 }
20
TSHistoryUpdate(TSHistory tsh,PetscInt id,PetscReal time)21 PetscErrorCode TSHistoryUpdate(TSHistory tsh, PetscInt id, PetscReal time)
22 {
23 PetscFunctionBegin;
24 if (tsh->n == tsh->c) { /* reallocation */
25 tsh->c += tsh->s;
26 PetscCall(PetscRealloc(tsh->c * sizeof(*tsh->hist), &tsh->hist));
27 PetscCall(PetscRealloc(tsh->c * sizeof(*tsh->hist_id), &tsh->hist_id));
28 }
29 tsh->sorted = (PetscBool)(tsh->sorted && (tsh->n ? (PetscBool)(time >= tsh->hist[tsh->n - 1]) : PETSC_TRUE));
30 #if defined(PETSC_USE_DEBUG)
31 if (tsh->n) { /* id should be unique */
32 PetscInt loc, *ids;
33
34 PetscCall(PetscMalloc1(tsh->n, &ids));
35 PetscCall(PetscArraycpy(ids, tsh->hist_id, tsh->n));
36 PetscCall(PetscSortInt(tsh->n, ids));
37 PetscCall(PetscFindInt(id, tsh->n, ids, &loc));
38 PetscCall(PetscFree(ids));
39 PetscCheck(loc < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "History id should be unique");
40 }
41 #endif
42 tsh->hist[tsh->n] = time;
43 tsh->hist_id[tsh->n] = id;
44 tsh->n += 1;
45 PetscFunctionReturn(PETSC_SUCCESS);
46 }
47
TSHistoryGetTime(TSHistory tsh,PetscBool backward,PetscInt step,PetscReal * t)48 PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t)
49 {
50 PetscFunctionBegin;
51 if (!t) PetscFunctionReturn(PETSC_SUCCESS);
52 PetscAssertPointer(t, 4);
53 if (!tsh->sorted) {
54 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id));
55 tsh->sorted = PETSC_TRUE;
56 }
57 PetscCheck(step >= 0 && step < tsh->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Given time step %" PetscInt_FMT " does not match any in history [0,%" PetscCount_FMT "]", step, tsh->n);
58 if (!backward) *t = tsh->hist[step];
59 else *t = tsh->hist[tsh->n - step - 1];
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
TSHistoryGetTimeStep(TSHistory tsh,PetscBool backward,PetscInt step,PetscReal * dt)63 PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt)
64 {
65 PetscFunctionBegin;
66 if (!dt) PetscFunctionReturn(PETSC_SUCCESS);
67 PetscAssertPointer(dt, 4);
68 if (!tsh->sorted) {
69 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id));
70 tsh->sorted = PETSC_TRUE;
71 }
72 PetscCheck(step >= 0 && step <= tsh->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Given time step %" PetscInt_FMT " does not match any in history [0,%" PetscCount_FMT "]", step, tsh->n);
73 if (!backward) *dt = tsh->hist[PetscMin(step + 1, tsh->n - 1)] - tsh->hist[PetscMin(step, tsh->n - 1)];
74 else *dt = tsh->hist[PetscMax(tsh->n - step - 1, 0)] - tsh->hist[PetscMax(tsh->n - step - 2, 0)];
75 PetscFunctionReturn(PETSC_SUCCESS);
76 }
77
TSHistoryGetLocFromTime(TSHistory tsh,PetscReal time,PetscInt * loc)78 PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc)
79 {
80 PetscFunctionBegin;
81 PetscAssertPointer(loc, 3);
82 if (!tsh->sorted) {
83 PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id));
84 tsh->sorted = PETSC_TRUE;
85 }
86 PetscCall(PetscFindReal(time, tsh->n, tsh->hist, PETSC_SMALL, loc));
87 PetscFunctionReturn(PETSC_SUCCESS);
88 }
89
TSHistorySetHistory(TSHistory tsh,PetscInt n,PetscReal hist[],PetscInt hist_id[],PetscBool sorted)90 PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted)
91 {
92 PetscFunctionBegin;
93 PetscValidLogicalCollectiveIntComm(tsh->comm, n, 2);
94 PetscCheck(n >= 0, tsh->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot request a negative size for history storage");
95 if (n) PetscAssertPointer(hist, 3);
96 PetscCall(PetscFree(tsh->hist));
97 PetscCall(PetscFree(tsh->hist_id));
98 tsh->n = (size_t)n;
99 tsh->c = (size_t)n;
100 PetscCall(PetscMalloc1(tsh->n, &tsh->hist));
101 PetscCall(PetscMalloc1(tsh->n, &tsh->hist_id));
102 for (PetscInt i = 0; i < n; i++) {
103 tsh->hist[i] = hist[i];
104 tsh->hist_id[i] = hist_id ? hist_id[i] : i;
105 }
106 if (!sorted) PetscCall(PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id));
107 tsh->sorted = PETSC_TRUE;
108 PetscFunctionReturn(PETSC_SUCCESS);
109 }
110
TSHistoryGetHistory(TSHistory tsh,PetscInt * n,const PetscReal * hist[],const PetscInt * hist_id[],PetscBool * sorted)111 PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal *hist[], const PetscInt *hist_id[], PetscBool *sorted)
112 {
113 PetscFunctionBegin;
114 if (n) PetscCall(PetscIntCast(tsh->n, n));
115 if (hist) *hist = tsh->hist;
116 if (hist_id) *hist_id = tsh->hist_id;
117 if (sorted) *sorted = tsh->sorted;
118 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120
TSHistoryDestroy(TSHistory * tsh)121 PetscErrorCode TSHistoryDestroy(TSHistory *tsh)
122 {
123 PetscFunctionBegin;
124 if (!*tsh) PetscFunctionReturn(PETSC_SUCCESS);
125 PetscCall(PetscFree((*tsh)->hist));
126 PetscCall(PetscFree((*tsh)->hist_id));
127 PetscCall(PetscCommDestroy(&(*tsh)->comm));
128 PetscCall(PetscFree(*tsh));
129 *tsh = NULL;
130 PetscFunctionReturn(PETSC_SUCCESS);
131 }
132
TSHistoryCreate(MPI_Comm comm,TSHistory * hst)133 PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst)
134 {
135 TSHistory tsh;
136
137 PetscFunctionBegin;
138 PetscAssertPointer(hst, 2);
139 PetscCall(PetscNew(&tsh));
140 PetscCall(PetscCommDuplicate(comm, &tsh->comm, NULL));
141
142 tsh->c = 1024; /* capacity */
143 tsh->s = 1024; /* reallocation size */
144 tsh->sorted = PETSC_TRUE;
145
146 PetscCall(PetscMalloc1(tsh->c, &tsh->hist));
147 PetscCall(PetscMalloc1(tsh->c, &tsh->hist_id));
148 *hst = tsh;
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151