Actual source code: tshistory.c

      1: #include <petsc/private/tshistoryimpl.h>

      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: };

     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: }

     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;

     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: }

     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: }

     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: }

     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: }

     90: PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted)
     91: {
     92:   PetscFunctionBegin;
     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: }

    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: }

    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: }

    133: PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst)
    134: {
    135:   TSHistory tsh;

    137:   PetscFunctionBegin;
    138:   PetscAssertPointer(hst, 2);
    139:   PetscCall(PetscNew(&tsh));
    140:   PetscCall(PetscCommDuplicate(comm, &tsh->comm, NULL));

    142:   tsh->c      = 1024; /* capacity */
    143:   tsh->s      = 1024; /* reallocation size */
    144:   tsh->sorted = PETSC_TRUE;

    146:   PetscCall(PetscMalloc1(tsh->c, &tsh->hist));
    147:   PetscCall(PetscMalloc1(tsh->c, &tsh->hist_id));
    148:   *hst = tsh;
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }