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; b1[1] = b;                                              \
 10:     MPIU_Allreduce(b1,b2,2,MPIU_INT,MPI_MAX,a); \
 12:   } while (0)

 15:   do {                                                                  \
 16:     PetscMPIInt b1[2],b2[2];                                            \
 17:     b1[0] = -(PetscMPIInt)b; b1[1] = (PetscMPIInt)b;                    \
 18:     MPIU_Allreduce(b1,b2,2,MPI_INT,MPI_MAX,a); \
 20:   } while (0)

 23:   do {                                                                  \
 24:     PetscReal b1[3],b2[3];                                              \
 25:     if (PetscIsNanReal(b)) {b1[2] = 1;} else {b1[2] = 0;};              \
 26:     b1[0] = -b; b1[1] = b;                                              \
 27:     MPI_Allreduce(b1,b2,3,MPIU_REAL,MPIU_MAX,a); \
 29:   } while (0)

 31: #else


 37: #endif

 39: struct _n_TSHistory {
 40:   MPI_Comm  comm;     /* used for runtime collective checks */
 41:   PetscReal *hist;    /* time history */
 42:   PetscInt  *hist_id; /* stores the stepid in time history */
 43:   size_t    n;        /* current number of steps registered */
 44:   PetscBool sorted;   /* if the history is sorted in ascending order */
 45:   size_t    c;        /* current capacity of history */
 46:   size_t    s;        /* reallocation size */
 47: };

 49: PetscErrorCode TSHistoryGetNumSteps(TSHistory tsh, PetscInt *n)
 50: {
 52:   *n = tsh->n;
 53:   return 0;
 54: }

 56: PetscErrorCode TSHistoryUpdate(TSHistory tsh, PetscInt id, PetscReal time)
 57: {
 58:   if (tsh->n == tsh->c) { /* reallocation */
 59:     tsh->c += tsh->s;
 60:     PetscRealloc(tsh->c*sizeof(*tsh->hist),&tsh->hist);
 61:     PetscRealloc(tsh->c*sizeof(*tsh->hist_id),&tsh->hist_id);
 62:   }
 63:   tsh->sorted = (PetscBool)(tsh->sorted && (tsh->n ? time >= tsh->hist[tsh->n-1] : PETSC_TRUE));
 64: #if defined(PETSC_USE_DEBUG)
 65:   if (tsh->n) { /* id should be unique */
 66:     PetscInt loc,*ids;

 68:     PetscMalloc1(tsh->n,&ids);
 69:     PetscArraycpy(ids,tsh->hist_id,tsh->n);
 70:     PetscSortInt(tsh->n,ids);
 71:     PetscFindInt(id,tsh->n,ids,&loc);
 72:     PetscFree(ids);
 74:   }
 75: #endif
 76:   tsh->hist[tsh->n]    = time;
 77:   tsh->hist_id[tsh->n] = id;
 78:   tsh->n += 1;
 79:   return 0;
 80: }

 82: PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t)
 83: {
 84:   if (!t) return 0;
 86:   if (!tsh->sorted) {

 88:     PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);
 89:     tsh->sorted = PETSC_TRUE;
 90:   }
 92:   if (!backward) *t = tsh->hist[step];
 93:   else           *t = tsh->hist[tsh->n-step-1];
 94:   return 0;
 95: }

 97: PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt)
 98: {
 99:   if (!dt) return 0;
101:   if (!tsh->sorted) {

103:     PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);
104:     tsh->sorted = PETSC_TRUE;
105:   }
107:   if (!backward) *dt = tsh->hist[PetscMin(step+1,(PetscInt)tsh->n-1)] - tsh->hist[PetscMin(step,(PetscInt)tsh->n-1)];
108:   else           *dt = tsh->hist[PetscMax((PetscInt)tsh->n-step-1,0)] - tsh->hist[PetscMax((PetscInt)tsh->n-step-2,0)];
109:   return 0;
110: }

112: PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc)
113: {
115:   if (!tsh->sorted) {
116:     PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);
117:     tsh->sorted = PETSC_TRUE;
118:   }
119:   PetscFindReal(time,tsh->n,tsh->hist,PETSC_SMALL,loc);
120:   return 0;
121: }

123: PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted)
124: {
128:   PetscFree(tsh->hist);
129:   PetscFree(tsh->hist_id);
130:   tsh->n = (size_t) n;
131:   tsh->c = (size_t) n;
132:   PetscMalloc1(tsh->n,&tsh->hist);
133:   PetscMalloc1(tsh->n,&tsh->hist_id);
134:   for (PetscInt i = 0; i < (PetscInt)tsh->n; i++) {
135:     tsh->hist[i]    = hist[i];
136:     tsh->hist_id[i] = hist_id ? hist_id[i] : i;
137:   }
138:   if (!sorted) PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);
139:   tsh->sorted = PETSC_TRUE;
140:   return 0;
141: }

143: PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal* hist[], const PetscInt* hist_id[], PetscBool *sorted)
144: {
145:   if (n)             *n = tsh->n;
146:   if (hist)       *hist = tsh->hist;
147:   if (hist_id) *hist_id = tsh->hist_id;
148:   if (sorted)   *sorted = tsh->sorted;
149:   return 0;
150: }

152: PetscErrorCode TSHistoryDestroy(TSHistory *tsh)
153: {
154:   if (!*tsh) return 0;
155:   PetscFree((*tsh)->hist);
156:   PetscFree((*tsh)->hist_id);
157:   PetscCommDestroy(&((*tsh)->comm));
158:   PetscFree((*tsh));
159:   *tsh = NULL;
160:   return 0;
161: }

163: PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst)
164: {
165:   TSHistory      tsh;

168:   *hst = NULL;
169:   PetscNew(&tsh);
170:   PetscCommDuplicate(comm,&tsh->comm,NULL);

172:   tsh->c      = 1024; /* capacity */
173:   tsh->s      = 1024; /* reallocation size */
174:   tsh->sorted = PETSC_TRUE;

176:   PetscMalloc1(tsh->c,&tsh->hist);
177:   PetscMalloc1(tsh->c,&tsh->hist_id);
178:   *hst = tsh;
179:   return 0;
180: }