1 #include <petsc/private/tshistoryimpl.h> 2 3 /* These macros can be moved to petscimpl.h eventually */ 4 #if defined(PETSC_USE_DEBUG) 5 6 #define PetscValidLogicalCollectiveIntComm(a,b,c) \ 7 do { \ 8 PetscInt b1[2],b2[2]; \ 9 b1[0] = -b; b1[1] = b; \ 10 PetscCallMPI(MPIU_Allreduce(b1,b2,2,MPIU_INT,MPI_MAX,a)); \ 11 PetscCheck(-b2[0] == b2[1],a,PETSC_ERR_ARG_WRONG,"Int value must be same on all processes, argument # %d",c); \ 12 } while (0) 13 14 #define PetscValidLogicalCollectiveBoolComm(a,b,c) \ 15 do { \ 16 PetscMPIInt b1[2],b2[2]; \ 17 b1[0] = -(PetscMPIInt)b; b1[1] = (PetscMPIInt)b; \ 18 PetscCallMPI(MPIU_Allreduce(b1,b2,2,MPI_INT,MPI_MAX,a)); \ 19 PetscCheck(-b2[0] == b2[1],a,PETSC_ERR_ARG_WRONG,"Bool value must be same on all processes, argument # %d",c); \ 20 } while (0) 21 22 #define PetscValidLogicalCollectiveRealComm(a,b,c) \ 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 PetscCallMPI(MPI_Allreduce(b1,b2,3,MPIU_REAL,MPIU_MAX,a)); \ 28 PetscCheck((b2[2] == 1) || PetscEqualReal(-b2[0],b2[1]),a,PETSC_ERR_ARG_WRONG,"Real value must be same on all processes, argument # %d",c); \ 29 } while (0) 30 31 #else 32 33 #define PetscValidLogicalCollectiveRealComm(a,b,c) do {} while (0) 34 #define PetscValidLogicalCollectiveIntComm(a,b,c) do {} while (0) 35 #define PetscValidLogicalCollectiveBoolComm(a,b,c) do {} while (0) 36 37 #endif 38 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 }; 48 49 PetscErrorCode TSHistoryGetNumSteps(TSHistory tsh, PetscInt *n) 50 { 51 PetscFunctionBegin; 52 PetscValidIntPointer(n,2); 53 *n = tsh->n; 54 PetscFunctionReturn(0); 55 } 56 57 PetscErrorCode TSHistoryUpdate(TSHistory tsh, PetscInt id, PetscReal time) 58 { 59 PetscFunctionBegin; 60 if (tsh->n == tsh->c) { /* reallocation */ 61 tsh->c += tsh->s; 62 PetscCall(PetscRealloc(tsh->c*sizeof(*tsh->hist),&tsh->hist)); 63 PetscCall(PetscRealloc(tsh->c*sizeof(*tsh->hist_id),&tsh->hist_id)); 64 } 65 tsh->sorted = (PetscBool)(tsh->sorted && (tsh->n ? time >= tsh->hist[tsh->n-1] : PETSC_TRUE)); 66 #if defined(PETSC_USE_DEBUG) 67 if (tsh->n) { /* id should be unique */ 68 PetscInt loc,*ids; 69 70 PetscCall(PetscMalloc1(tsh->n,&ids)); 71 PetscCall(PetscArraycpy(ids,tsh->hist_id,tsh->n)); 72 PetscCall(PetscSortInt(tsh->n,ids)); 73 PetscCall(PetscFindInt(id,tsh->n,ids,&loc)); 74 PetscCall(PetscFree(ids)); 75 PetscCheck(loc < 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"History id should be unique"); 76 } 77 #endif 78 tsh->hist[tsh->n] = time; 79 tsh->hist_id[tsh->n] = id; 80 tsh->n += 1; 81 PetscFunctionReturn(0); 82 } 83 84 PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t) 85 { 86 PetscFunctionBegin; 87 if (!t) PetscFunctionReturn(0); 88 PetscValidRealPointer(t,4); 89 if (!tsh->sorted) { 90 91 PetscCall(PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id)); 92 tsh->sorted = PETSC_TRUE; 93 } 94 PetscCheck(step >= 0 && step < (PetscInt)tsh->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Given time step %D does not match any in history [0,%D]",step,(PetscInt)tsh->n); 95 if (!backward) *t = tsh->hist[step]; 96 else *t = tsh->hist[tsh->n-step-1]; 97 PetscFunctionReturn(0); 98 } 99 100 PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt) 101 { 102 PetscFunctionBegin; 103 if (!dt) PetscFunctionReturn(0); 104 PetscValidRealPointer(dt,4); 105 if (!tsh->sorted) { 106 107 PetscCall(PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id)); 108 tsh->sorted = PETSC_TRUE; 109 } 110 PetscCheck(step >= 0 && step <= (PetscInt)tsh->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Given time step %D does not match any in history [0,%D]",step,(PetscInt)tsh->n); 111 if (!backward) *dt = tsh->hist[PetscMin(step+1,(PetscInt)tsh->n-1)] - tsh->hist[PetscMin(step,(PetscInt)tsh->n-1)]; 112 else *dt = tsh->hist[PetscMax((PetscInt)tsh->n-step-1,0)] - tsh->hist[PetscMax((PetscInt)tsh->n-step-2,0)]; 113 PetscFunctionReturn(0); 114 } 115 116 PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc) 117 { 118 PetscFunctionBegin; 119 PetscValidIntPointer(loc,3); 120 if (!tsh->sorted) { 121 PetscCall(PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id)); 122 tsh->sorted = PETSC_TRUE; 123 } 124 PetscCall(PetscFindReal(time,tsh->n,tsh->hist,PETSC_SMALL,loc)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted) 129 { 130 PetscFunctionBegin; 131 PetscValidLogicalCollectiveIntComm(tsh->comm,n,2); 132 PetscCheck(n >= 0,tsh->comm,PETSC_ERR_ARG_OUTOFRANGE,"Cannot request a negative size for history storage"); 133 if (n) PetscValidRealPointer(hist,3); 134 PetscCall(PetscFree(tsh->hist)); 135 PetscCall(PetscFree(tsh->hist_id)); 136 tsh->n = (size_t) n; 137 tsh->c = (size_t) n; 138 PetscCall(PetscMalloc1(tsh->n,&tsh->hist)); 139 PetscCall(PetscMalloc1(tsh->n,&tsh->hist_id)); 140 for (PetscInt i = 0; i < (PetscInt)tsh->n; i++) { 141 tsh->hist[i] = hist[i]; 142 tsh->hist_id[i] = hist_id ? hist_id[i] : i; 143 } 144 if (!sorted) PetscCall(PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id)); 145 tsh->sorted = PETSC_TRUE; 146 PetscFunctionReturn(0); 147 } 148 149 PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal* hist[], const PetscInt* hist_id[], PetscBool *sorted) 150 { 151 PetscFunctionBegin; 152 if (n) *n = tsh->n; 153 if (hist) *hist = tsh->hist; 154 if (hist_id) *hist_id = tsh->hist_id; 155 if (sorted) *sorted = tsh->sorted; 156 PetscFunctionReturn(0); 157 } 158 159 PetscErrorCode TSHistoryDestroy(TSHistory *tsh) 160 { 161 PetscFunctionBegin; 162 if (!*tsh) PetscFunctionReturn(0); 163 PetscCall(PetscFree((*tsh)->hist)); 164 PetscCall(PetscFree((*tsh)->hist_id)); 165 PetscCall(PetscCommDestroy(&((*tsh)->comm))); 166 PetscCall(PetscFree((*tsh))); 167 *tsh = NULL; 168 PetscFunctionReturn(0); 169 } 170 171 PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst) 172 { 173 TSHistory tsh; 174 175 PetscFunctionBegin; 176 PetscValidPointer(hst,2); 177 *hst = NULL; 178 PetscCall(PetscNew(&tsh)); 179 PetscCall(PetscCommDuplicate(comm,&tsh->comm,NULL)); 180 181 tsh->c = 1024; /* capacity */ 182 tsh->s = 1024; /* reallocation size */ 183 tsh->sorted = PETSC_TRUE; 184 185 PetscCall(PetscMalloc1(tsh->c,&tsh->hist)); 186 PetscCall(PetscMalloc1(tsh->c,&tsh->hist_id)); 187 *hst = tsh; 188 PetscFunctionReturn(0); 189 } 190