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 PetscErrorCode _7_ierr; \ 9 PetscInt b1[2],b2[2]; \ 10 b1[0] = -b; b1[1] = b; \ 11 _7_ierr = MPIU_Allreduce(b1,b2,2,MPIU_INT,MPI_MAX,a);CHKERRMPI(_7_ierr); \ 12 PetscCheckFalse(-b2[0] != b2[1],a,PETSC_ERR_ARG_WRONG,"Int value must be same on all processes, argument # %d",c); \ 13 } while (0) 14 15 #define PetscValidLogicalCollectiveBoolComm(a,b,c) \ 16 do { \ 17 PetscErrorCode _7_ierr; \ 18 PetscMPIInt b1[2],b2[2]; \ 19 b1[0] = -(PetscMPIInt)b; b1[1] = (PetscMPIInt)b; \ 20 _7_ierr = MPIU_Allreduce(b1,b2,2,MPI_INT,MPI_MAX,a);CHKERRMPI(_7_ierr); \ 21 PetscCheckFalse(-b2[0] != b2[1],a,PETSC_ERR_ARG_WRONG,"Bool value must be same on all processes, argument # %d",c); \ 22 } while (0) 23 24 #define PetscValidLogicalCollectiveRealComm(a,b,c) \ 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_ierr = MPI_Allreduce(b1,b2,3,MPIU_REAL,MPIU_MAX,a);CHKERRMPI(_7_ierr); \ 31 PetscCheckFalse(!(b2[2] > 0) && !PetscEqualReal(-b2[0],b2[1]),a,PETSC_ERR_ARG_WRONG,"Real value must be same on all processes, argument # %d",c); \ 32 } while (0) 33 34 #else 35 36 #define PetscValidLogicalCollectiveRealComm(a,b,c) do {} while (0) 37 #define PetscValidLogicalCollectiveIntComm(a,b,c) do {} while (0) 38 #define PetscValidLogicalCollectiveBoolComm(a,b,c) do {} while (0) 39 40 #endif 41 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 size_t n; /* current number of steps registered */ 47 PetscBool sorted; /* if the history is sorted in ascending order */ 48 size_t c; /* current capacity of history */ 49 size_t s; /* reallocation size */ 50 }; 51 52 PetscErrorCode TSHistoryGetNumSteps(TSHistory tsh, PetscInt *n) 53 { 54 PetscFunctionBegin; 55 PetscValidPointer(n,2); 56 *n = tsh->n; 57 PetscFunctionReturn(0); 58 } 59 60 PetscErrorCode TSHistoryUpdate(TSHistory tsh, PetscInt id, PetscReal time) 61 { 62 PetscErrorCode ierr; 63 64 PetscFunctionBegin; 65 PetscValidLogicalCollectiveIntComm(tsh->comm,id,2); 66 PetscValidLogicalCollectiveRealComm(tsh->comm,time,3); 67 if (tsh->n == tsh->c) { /* reallocation */ 68 tsh->c += tsh->s; 69 ierr = PetscRealloc(tsh->c*sizeof(*tsh->hist),&tsh->hist);CHKERRQ(ierr); 70 ierr = PetscRealloc(tsh->c*sizeof(*tsh->hist_id),&tsh->hist_id);CHKERRQ(ierr); 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; 76 77 ierr = PetscMalloc1(tsh->n,&ids);CHKERRQ(ierr); 78 ierr = PetscArraycpy(ids,tsh->hist_id,tsh->n);CHKERRQ(ierr); 79 ierr = PetscSortInt(tsh->n,ids);CHKERRQ(ierr); 80 ierr = PetscFindInt(id,tsh->n,ids,&loc);CHKERRQ(ierr); 81 ierr = PetscFree(ids);CHKERRQ(ierr); 82 PetscCheckFalse(loc >=0,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 PetscFunctionReturn(0); 89 } 90 91 PetscErrorCode TSHistoryGetTime(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *t) 92 { 93 PetscFunctionBegin; 94 PetscValidLogicalCollectiveBoolComm(tsh->comm,backward,2); 95 PetscValidLogicalCollectiveIntComm(tsh->comm,step,3); 96 if (!t) PetscFunctionReturn(0); 97 PetscValidRealPointer(t,4); 98 if (!tsh->sorted) { 99 PetscErrorCode ierr; 100 101 ierr = PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);CHKERRQ(ierr); 102 tsh->sorted = PETSC_TRUE; 103 } 104 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); 105 if (!backward) *t = tsh->hist[step]; 106 else *t = tsh->hist[tsh->n-step-1]; 107 PetscFunctionReturn(0); 108 } 109 110 PetscErrorCode TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt) 111 { 112 PetscFunctionBegin; 113 PetscValidLogicalCollectiveBoolComm(tsh->comm,backward,2); 114 PetscValidLogicalCollectiveIntComm(tsh->comm,step,3); 115 if (!dt) PetscFunctionReturn(0); 116 PetscValidRealPointer(dt,4); 117 if (!tsh->sorted) { 118 PetscErrorCode ierr; 119 120 ierr = PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);CHKERRQ(ierr); 121 tsh->sorted = PETSC_TRUE; 122 } 123 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); 124 if (!backward) *dt = tsh->hist[PetscMin(step+1,(PetscInt)tsh->n-1)] - tsh->hist[PetscMin(step,(PetscInt)tsh->n-1)]; 125 else *dt = tsh->hist[PetscMax((PetscInt)tsh->n-step-1,0)] - tsh->hist[PetscMax((PetscInt)tsh->n-step-2,0)]; 126 PetscFunctionReturn(0); 127 } 128 129 PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc) 130 { 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 PetscValidLogicalCollectiveRealComm(tsh->comm,time,2); 135 PetscValidIntPointer(loc,3); 136 if (!tsh->sorted) { 137 ierr = PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);CHKERRQ(ierr); 138 tsh->sorted = PETSC_TRUE; 139 } 140 ierr = PetscFindReal(time,tsh->n,tsh->hist,PETSC_SMALL,loc);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 144 PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted) 145 { 146 PetscInt i; 147 PetscErrorCode ierr; 148 149 PetscFunctionBegin; 150 PetscValidLogicalCollectiveIntComm(tsh->comm,n,2); 151 if (n < 0) SETERRQ(tsh->comm,PETSC_ERR_ARG_OUTOFRANGE,"Cannot request a negative size for history storage"); 152 if (n) PetscValidRealPointer(hist,3); 153 ierr = PetscFree(tsh->hist);CHKERRQ(ierr); 154 ierr = PetscFree(tsh->hist_id);CHKERRQ(ierr); 155 tsh->n = (size_t) n; 156 tsh->c = (size_t) n; 157 ierr = PetscMalloc1(tsh->n,&tsh->hist);CHKERRQ(ierr); 158 ierr = PetscMalloc1(tsh->n,&tsh->hist_id);CHKERRQ(ierr); 159 for (i = 0; i < (PetscInt)tsh->n; i++) { 160 tsh->hist[i] = hist[i]; 161 tsh->hist_id[i] = hist_id ? hist_id[i] : i; 162 } 163 if (!sorted) { 164 ierr = PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);CHKERRQ(ierr); 165 } 166 tsh->sorted = PETSC_TRUE; 167 PetscFunctionReturn(0); 168 } 169 170 PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal* hist[], const PetscInt* hist_id[], PetscBool *sorted) 171 { 172 PetscFunctionBegin; 173 if (n) *n = tsh->n; 174 if (hist) *hist = tsh->hist; 175 if (hist_id) *hist_id = tsh->hist_id; 176 if (sorted) *sorted = tsh->sorted; 177 PetscFunctionReturn(0); 178 } 179 180 PetscErrorCode TSHistoryDestroy(TSHistory *tsh) 181 { 182 PetscErrorCode ierr; 183 184 PetscFunctionBegin; 185 if (!*tsh) PetscFunctionReturn(0); 186 ierr = PetscFree((*tsh)->hist);CHKERRQ(ierr); 187 ierr = PetscFree((*tsh)->hist_id);CHKERRQ(ierr); 188 ierr = PetscCommDestroy(&((*tsh)->comm));CHKERRQ(ierr); 189 ierr = PetscFree((*tsh));CHKERRQ(ierr); 190 *tsh = NULL; 191 PetscFunctionReturn(0); 192 } 193 194 PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst) 195 { 196 TSHistory tsh; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 PetscValidPointer(hst,2); 201 *hst = NULL; 202 ierr = PetscNew(&tsh);CHKERRQ(ierr); 203 ierr = PetscCommDuplicate(comm,&tsh->comm,NULL);CHKERRQ(ierr); 204 205 tsh->c = 1024; /* capacity */ 206 tsh->s = 1024; /* reallocation size */ 207 tsh->sorted = PETSC_TRUE; 208 209 ierr = PetscMalloc1(tsh->c,&tsh->hist);CHKERRQ(ierr); 210 ierr = PetscMalloc1(tsh->c,&tsh->hist_id);CHKERRQ(ierr); 211 *hst = tsh; 212 PetscFunctionReturn(0); 213 } 214