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:       MPIU_Allreduce(b1, b2, 2, MPIU_INT, MPI_MAX, a); \
 13:     } while (0)

 16:     do { \
 17:       PetscMPIInt b1[2], b2[2]; \
 18:       b1[0] = -(PetscMPIInt)b; \
 19:       b1[1] = (PetscMPIInt)b; \
 20:       MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, a); \
 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:       MPI_Allreduce(b1, b2, 3, MPIU_REAL, MPIU_MAX, a); \
 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:   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: };

 62: PetscErrorCode TSHistoryGetNumSteps(TSHistory tsh, PetscInt *n)
 63: {
 65:   *n = tsh->n;
 66:   return 0;
 67: }

 69: PetscErrorCode TSHistoryUpdate(TSHistory tsh, PetscInt id, PetscReal time)
 70: {
 71:   if (tsh->n == tsh->c) { /* reallocation */
 72:     tsh->c += tsh->s;
 73:     PetscRealloc(tsh->c * sizeof(*tsh->hist), &tsh->hist);
 74:     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;

 81:     PetscMalloc1(tsh->n, &ids);
 82:     PetscArraycpy(ids, tsh->hist_id, tsh->n);
 83:     PetscSortInt(tsh->n, ids);
 84:     PetscFindInt(id, tsh->n, ids, &loc);
 85:     PetscFree(ids);
 87:   }
 88: #endif
 89:   tsh->hist[tsh->n]    = time;
 90:   tsh->hist_id[tsh->n] = id;
 91:   tsh->n += 1;
 92:   return 0;
 93: }

 95: PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t)
 96: {
 97:   if (!t) return 0;
 99:   if (!tsh->sorted) {
100:     PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id);
101:     tsh->sorted = PETSC_TRUE;
102:   }
104:   if (!backward) *t = tsh->hist[step];
105:   else *t = tsh->hist[tsh->n - step - 1];
106:   return 0;
107: }

109: PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt)
110: {
111:   if (!dt) return 0;
113:   if (!tsh->sorted) {
114:     PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id);
115:     tsh->sorted = PETSC_TRUE;
116:   }
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:   return 0;
121: }

123: PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc)
124: {
126:   if (!tsh->sorted) {
127:     PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id);
128:     tsh->sorted = PETSC_TRUE;
129:   }
130:   PetscFindReal(time, tsh->n, tsh->hist, PETSC_SMALL, loc);
131:   return 0;
132: }

134: PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted)
135: {
139:   PetscFree(tsh->hist);
140:   PetscFree(tsh->hist_id);
141:   tsh->n = (size_t)n;
142:   tsh->c = (size_t)n;
143:   PetscMalloc1(tsh->n, &tsh->hist);
144:   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) PetscSortRealWithArrayInt(tsh->n, tsh->hist, tsh->hist_id);
150:   tsh->sorted = PETSC_TRUE;
151:   return 0;
152: }

154: PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal *hist[], const PetscInt *hist_id[], PetscBool *sorted)
155: {
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:   return 0;
161: }

163: PetscErrorCode TSHistoryDestroy(TSHistory *tsh)
164: {
165:   if (!*tsh) return 0;
166:   PetscFree((*tsh)->hist);
167:   PetscFree((*tsh)->hist_id);
168:   PetscCommDestroy(&((*tsh)->comm));
169:   PetscFree((*tsh));
170:   *tsh = NULL;
171:   return 0;
172: }

174: PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst)
175: {
176:   TSHistory tsh;

179:   *hst = NULL;
180:   PetscNew(&tsh);
181:   PetscCommDuplicate(comm, &tsh->comm, NULL);

183:   tsh->c      = 1024; /* capacity */
184:   tsh->s      = 1024; /* reallocation size */
185:   tsh->sorted = PETSC_TRUE;

187:   PetscMalloc1(tsh->c, &tsh->hist);
188:   PetscMalloc1(tsh->c, &tsh->hist_id);
189:   *hst = tsh;
190:   return 0;
191: }