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);CHKERRQ(_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 = PetscMemcpy(ids,tsh->hist_id,tsh->n*sizeof(PetscInt));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 TSHistoryGetTimeStep(TSHistory tsh, PetscBool backward, PetscInt step, PetscReal *dt) 92 { 93 PetscFunctionBegin; 94 PetscValidLogicalCollectiveBoolComm(tsh->comm,backward,2); 95 PetscValidLogicalCollectiveIntComm(tsh->comm,step,3); 96 PetscValidRealPointer(dt,4); 97 if (!tsh->sorted) { 98 PetscErrorCode ierr; 99 100 ierr = PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);CHKERRQ(ierr); 101 tsh->sorted = PETSC_TRUE; 102 } 103 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); 104 if (!backward) *dt = tsh->hist[PetscMin(step+1,tsh->n-1)] - tsh->hist[PetscMin(step,tsh->n-1)]; 105 else *dt = tsh->hist[PetscMax(tsh->n-step-1,0)] - tsh->hist[PetscMax(tsh->n-step-2,0)]; 106 PetscFunctionReturn(0); 107 } 108 109 PetscErrorCode TSHistoryGetLocFromTime(TSHistory tsh, PetscReal time, PetscInt *loc) 110 { 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 PetscValidLogicalCollectiveRealComm(tsh->comm,time,2); 115 PetscValidIntPointer(loc,3); 116 if (!tsh->sorted) { 117 ierr = PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);CHKERRQ(ierr); 118 tsh->sorted = PETSC_TRUE; 119 } 120 ierr = PetscFindReal(time,tsh->n,tsh->hist,PETSC_SMALL,loc);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 PetscErrorCode TSHistorySetHistory(TSHistory tsh, PetscInt n, PetscReal hist[], PetscInt hist_id[], PetscBool sorted) 125 { 126 PetscInt i; 127 PetscErrorCode ierr; 128 129 PetscFunctionBegin; 130 PetscValidLogicalCollectiveIntComm(tsh->comm,n,2); 131 if (n) PetscValidRealPointer(hist,3); 132 ierr = PetscFree(tsh->hist);CHKERRQ(ierr); 133 ierr = PetscFree(tsh->hist_id);CHKERRQ(ierr); 134 tsh->n = n; 135 tsh->c = n; 136 ierr = PetscMalloc1(tsh->n,&tsh->hist);CHKERRQ(ierr); 137 ierr = PetscMalloc1(tsh->n,&tsh->hist_id);CHKERRQ(ierr); 138 for (i = 0; i < tsh->n; i++) { 139 tsh->hist[i] = hist[i]; 140 tsh->hist_id[i] = hist_id ? hist_id[i] : i; 141 } 142 if (!sorted) { 143 ierr = PetscSortRealWithArrayInt(tsh->n,tsh->hist,tsh->hist_id);CHKERRQ(ierr); 144 } 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 PetscErrorCode ierr; 162 163 PetscFunctionBegin; 164 if (!*tsh) PetscFunctionReturn(0); 165 ierr = PetscFree((*tsh)->hist);CHKERRQ(ierr); 166 ierr = PetscFree((*tsh)->hist_id);CHKERRQ(ierr); 167 ierr = PetscCommDestroy(&((*tsh)->comm));CHKERRQ(ierr); 168 ierr = PetscFree((*tsh));CHKERRQ(ierr); 169 *tsh = NULL; 170 PetscFunctionReturn(0); 171 } 172 173 PetscErrorCode TSHistoryCreate(MPI_Comm comm, TSHistory *hst) 174 { 175 TSHistory tsh; 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 PetscValidPointer(hst,2); 180 *hst = NULL; 181 ierr = PetscNew(&tsh);CHKERRQ(ierr); 182 ierr = PetscCommDuplicate(comm,&tsh->comm,NULL);CHKERRQ(ierr); 183 184 tsh->c = 1024; /* capacity */ 185 tsh->s = 1024; /* reallocation size */ 186 tsh->sorted = PETSC_TRUE; 187 188 ierr = PetscMalloc1(tsh->c,&tsh->hist);CHKERRQ(ierr); 189 ierr = PetscMalloc1(tsh->c,&tsh->hist_id);CHKERRQ(ierr); 190 *hst = tsh; 191 PetscFunctionReturn(0); 192 } 193