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);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) 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);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) 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 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) 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 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 }; 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 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 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 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 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 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 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) PetscValidRealPointer(hist,3); 152 ierr = PetscFree(tsh->hist);CHKERRQ(ierr); 153 ierr = PetscFree(tsh->hist_id);CHKERRQ(ierr); 154 tsh->n = n; 155 tsh->c = n; 156 ierr = PetscMalloc1(tsh->n,&tsh->hist);CHKERRQ(ierr); 157 ierr = PetscMalloc1(tsh->n,&tsh->hist_id);CHKERRQ(ierr); 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 ierr = PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);CHKERRQ(ierr); 164 } 165 tsh->sorted = PETSC_TRUE; 166 PetscFunctionReturn(0); 167 } 168 169 PetscErrorCode TSHistoryGetHistory(TSHistory tsh, PetscInt *n, const PetscReal* hist[], const PetscInt* hist_id[], PetscBool *sorted) 170 { 171 PetscFunctionBegin; 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 PetscFunctionReturn(0); 177 } 178 179 PetscErrorCode TSHistoryDestroy(TSHistory *tsh) 180 { 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 if (!*tsh) PetscFunctionReturn(0); 185 ierr = PetscFree((*tsh)->hist);CHKERRQ(ierr); 186 ierr = PetscFree((*tsh)->hist_id);CHKERRQ(ierr); 187 ierr = PetscCommDestroy(&((*tsh)->comm));CHKERRQ(ierr); 188 ierr = PetscFree((*tsh));CHKERRQ(ierr); 189 *tsh = NULL; 190 PetscFunctionReturn(0); 191 } 192 193 PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst) 194 { 195 TSHistory tsh; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 PetscValidPointer(hst,2); 200 *hst = NULL; 201 ierr = PetscNew(&tsh);CHKERRQ(ierr); 202 ierr = PetscCommDuplicate(comm,&tsh->comm,NULL);CHKERRQ(ierr); 203 204 tsh->c = 1024; /* capacity */ 205 tsh->s = 1024; /* reallocation size */ 206 tsh->sorted = PETSC_TRUE; 207 208 ierr = PetscMalloc1(tsh->c,&tsh->hist);CHKERRQ(ierr); 209 ierr = PetscMalloc1(tsh->c,&tsh->hist_id);CHKERRQ(ierr); 210 *hst = tsh; 211 PetscFunctionReturn(0); 212 } 213