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:     PetscErrorCode _7_ierr;                                             \
      9:     PetscInt b1[2],b2[2];                                               \
     10:     b1[0] = -b; b1[1] = b;                                              \
     11:     _7_MPIU_Allreduce(b1,b2,2,MPIU_INT,MPI_MAX,a);CHKERRQ(_7_ierr); \
     12:     if (-b2[0] != b2[1]) SETERRQ1(a,PETSC_ERR_ARG_WRONG,"Int value must be same on all processes, argument # %d",c); \
     13:   } while (0)

     16:   do {                                                                  \
     17:     PetscErrorCode _7_ierr;                                             \
     18:     PetscMPIInt b1[2],b2[2];                                            \
     19:     b1[0] = -(PetscMPIInt)b; b1[1] = (PetscMPIInt)b;                    \
     20:     _7_MPIU_Allreduce(b1,b2,2,MPI_INT,MPI_MAX,a);CHKERRQ(_7_ierr); \
     21:     if (-b2[0] != b2[1]) SETERRQ1(a,PETSC_ERR_ARG_WRONG,"Bool value must be same on all processes, argument # %d",c); \
     22:   } while (0)

     25:   do {                                                                  \
     26:     PetscErrorCode _7_ierr;                                             \
     27:     PetscReal b1[3],b2[3];                                              \
     28:     if (PetscIsNanReal(b)) {b1[2] = 1;} else {b1[2] = 0;};              \
     29:     b1[0] = -b; b1[1] = b;                                              \
     30:     _7_MPI_Allreduce(b1,b2,3,MPIU_REAL,MPIU_MAX,a);CHKERRMPI(_7_ierr); \
     31:     if (!(b2[2] > 0) && !PetscEqualReal(-b2[0],b2[1])) SETERRQ1(a,PETSC_ERR_ARG_WRONG,"Real value must be same on all processes, argument # %d",c); \
     32:   } while (0)

     34: #else

     40: #endif

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

     52: PetscErrorCode TSHistoryGetNumSteps(TSHistory tsh, PetscInt *n)
     53: {
     56:   *n = tsh->n;
     57:   return(0);
     58: }

     60: PetscErrorCode TSHistoryUpdate(TSHistory tsh, PetscInt id, PetscReal time)
     61: {

     67:   if (tsh->n == tsh->c) { /* reallocation */
     68:     tsh->c += tsh->s;
     69:     PetscRealloc(tsh->c*sizeof(*tsh->hist),&tsh->hist);
     70:     PetscRealloc(tsh->c*sizeof(*tsh->hist_id),&tsh->hist_id);
     71:   }
     72:   tsh->sorted = (PetscBool)(tsh->sorted && (tsh->n ? time >= tsh->hist[tsh->n-1] : PETSC_TRUE));
     73: #if defined(PETSC_USE_DEBUG)
     74:   if (tsh->n) { /* id should be unique */
     75:     PetscInt loc,*ids;

     77:     PetscMalloc1(tsh->n,&ids);
     78:     PetscArraycpy(ids,tsh->hist_id,tsh->n);
     79:     PetscSortInt(tsh->n,ids);
     80:     PetscFindInt(id,tsh->n,ids,&loc);
     81:     PetscFree(ids);
     82:     if (loc >=0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"History id should be unique");
     83:   }
     84: #endif
     85:   tsh->hist[tsh->n]    = time;
     86:   tsh->hist_id[tsh->n] = id;
     87:   tsh->n += 1;
     88:   return(0);
     89: }

     91: PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t)
     92: {
     96:   if (!t) return(0);
     98:   if (!tsh->sorted) {

    101:     PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);
    102:     tsh->sorted = PETSC_TRUE;
    103:   }
    104:   if (step < 0 || step >= tsh->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Given time step %D does not match any in history [0,%D]",step,tsh->n);
    105:   if (!backward) *t = tsh->hist[step];
    106:   else           *t = tsh->hist[tsh->n-step-1];
    107:   return(0);
    108: }

    110: PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt)
    111: {
    115:   if (!dt) return(0);
    117:   if (!tsh->sorted) {

    120:     PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);
    121:     tsh->sorted = PETSC_TRUE;
    122:   }
    123:   if (step < 0 || step > tsh->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Given time step %D does not match any in history [0,%D]",step,tsh->n);
    124:   if (!backward) *dt = tsh->hist[PetscMin(step+1,tsh->n-1)] - tsh->hist[PetscMin(step,tsh->n-1)];
    125:   else           *dt = tsh->hist[PetscMax(tsh->n-step-1,0)] - tsh->hist[PetscMax(tsh->n-step-2,0)];
    126:   return(0);
    127: }

    129: PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc)
    130: {

    136:   if (!tsh->sorted) {
    137:     PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);
    138:     tsh->sorted = PETSC_TRUE;
    139:   }
    140:   PetscFindReal(time,tsh->n,tsh->hist,PETSC_SMALL,loc);
    141:   return(0);
    142: }

    144: PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted)
    145: {
    146:   PetscInt       i;

    152:   PetscFree(tsh->hist);
    153:   PetscFree(tsh->hist_id);
    154:   tsh->n = n;
    155:   tsh->c = n;
    156:   PetscMalloc1(tsh->n,&tsh->hist);
    157:   PetscMalloc1(tsh->n,&tsh->hist_id);
    158:   for (i = 0; i < tsh->n; i++) {
    159:     tsh->hist[i]    = hist[i];
    160:     tsh->hist_id[i] = hist_id ? hist_id[i] : i;
    161:   }
    162:   if (!sorted) {
    163:     PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);
    164:   }
    165:   tsh->sorted = PETSC_TRUE;
    166:   return(0);
    167: }

    169: PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal* hist[], const PetscInt* hist_id[], PetscBool *sorted)
    170: {
    172:   if (n)             *n = tsh->n;
    173:   if (hist)       *hist = tsh->hist;
    174:   if (hist_id) *hist_id = tsh->hist_id;
    175:   if (sorted)   *sorted = tsh->sorted;
    176:   return(0);
    177: }

    179: PetscErrorCode TSHistoryDestroy(TSHistory *tsh)
    180: {

    184:   if (!*tsh) return(0);
    185:   PetscFree((*tsh)->hist);
    186:   PetscFree((*tsh)->hist_id);
    187:   PetscCommDestroy(&((*tsh)->comm));
    188:   PetscFree((*tsh));
    189:   *tsh = NULL;
    190:   return(0);
    191: }

    193: PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst)
    194: {
    195:   TSHistory      tsh;

    200:   *hst = NULL;
    201:   PetscNew(&tsh);
    202:   PetscCommDuplicate(comm,&tsh->comm,NULL);

    204:   tsh->c      = 1024; /* capacity */
    205:   tsh->s      = 1024; /* reallocation size */
    206:   tsh->sorted = PETSC_TRUE;

    208:   PetscMalloc1(tsh->c,&tsh->hist);
    209:   PetscMalloc1(tsh->c,&tsh->hist_id);
    210:   *hst = tsh;
    211:   return(0);
    212: }