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