1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "DMDAGetLogicalCoordinate" 10 /*@C 11 DMDAGetLogicalCoordinate - Returns a the i,j,k logical coordinate for the closest mesh point to a x,y,z point in the coordinates of the DMDA 12 13 Collective on DMDA 14 15 Input Parameters: 16 + da - the distributed array 17 . x,y,z - the physical coordinates 18 - 19 20 Output Parameters: 21 + II, JJ, KK - the logical coordinate (-1 on processes that do not contain that point) 22 - X, Y, Z, - (optional) the coordinates of the located grid point 23 24 Level: advanced 25 26 Notes: 27 All processors that share the DMDA must call this with the same coordinate value 28 29 .keywords: distributed array, get, processor subset 30 @*/ 31 PetscErrorCode DMDAGetLogicalCoordinate(DM da,PetscScalar x,PetscScalar y,PetscScalar z,PetscInt *II,PetscInt *JJ,PetscInt *KK,PetscScalar *X,PetscScalar *Y,PetscScalar *Z) 32 { 33 DM_DA *dd = (DM_DA*)da->data; 34 PetscErrorCode ierr; 35 Vec coors; 36 DM dacoors; 37 DMDACoor2d **c; 38 PetscInt i,j,xs,xm,ys,ym; 39 PetscReal d,D = PETSC_MAX_REAL; 40 PetscMPIInt size; 41 42 PetscFunctionBegin; 43 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);CHKERRQ(ierr); 44 if (size > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently not in parallel"); 45 if (dd->dim == 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot get point from 1d DMDA"); 46 if (dd->dim == 3) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot get point from 3d DMDA"); 47 48 *II = -1; 49 if (JJ) *JJ = -1; 50 if (KK) *KK = -1; 51 52 ierr = DMGetCoordinateDM(da,&dacoors);CHKERRQ(ierr); 53 ierr = DMDAGetCorners(dacoors,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 54 ierr = DMGetCoordinates(da,&coors);CHKERRQ(ierr); 55 ierr = DMDAVecGetArray(dacoors,coors,&c);CHKERRQ(ierr); 56 for (j=ys; j<ys+ym; j++) { 57 for (i=xs; i<xs+xm; i++) { 58 d = PetscSqrtReal(PetscRealPart( (c[j][i].x - x)*(c[j][i].x - x) + (c[j][i].y - y)*(c[j][i].y - y) )); 59 if (d < D) { 60 D = d; 61 *II = i; 62 *JJ = j; 63 } 64 } 65 } 66 if (X) *X = c[*JJ][*II].x; 67 if (Y) *Y = c[*JJ][*II].y; 68 ierr = DMDAVecRestoreArray(dacoors,coors,&c);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "DMDAGetRay" 74 /*@C 75 DMDAGetRay - Returns a vector on process zero that contains a row or column of the values in a DMDA vector 76 77 Collective on DMDA 78 79 Input Parameters: 80 + da - the distributed array 81 . vec - the vector 82 . dir - Cartesian direction, either DMDA_X, DMDA_Y, or DMDA_Z 83 - gp - global grid point number in this direction 84 85 Output Parameters: 86 + newvec - the new vector that can hold the values (size zero on all processes except process 0) 87 - scatter - the VecScatter that will map from the original vector to the slice 88 89 Level: advanced 90 91 Notes: 92 All processors that share the DMDA must call this with the same gp value 93 94 .keywords: distributed array, get, processor subset 95 @*/ 96 PetscErrorCode DMDAGetRay(DM da,DMDADirection dir,PetscInt gp,Vec *newvec,VecScatter *scatter) 97 { 98 PetscMPIInt rank; 99 DM_DA *dd = (DM_DA*)da->data; 100 PetscErrorCode ierr; 101 IS is; 102 AO ao; 103 Vec vec; 104 PetscInt *indices,i,j; 105 106 PetscFunctionBegin; 107 if (dd->dim == 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot get slice from 1d DMDA"); 108 if (dd->dim == 3) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot get slice from 3d DMDA"); 109 ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 110 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);CHKERRQ(ierr); 111 if (!rank) { 112 if (dir == DMDA_Y) { 113 ierr = PetscMalloc(dd->w*dd->M*sizeof(PetscInt),&indices);CHKERRQ(ierr); 114 indices[0] = gp*dd->M*dd->w; 115 for (i=1; i<dd->M*dd->w; i++) indices[i] = indices[i-1] + 1; 116 117 ierr = AOApplicationToPetsc(ao,dd->M*dd->w,indices);CHKERRQ(ierr); 118 ierr = VecCreate(PETSC_COMM_SELF,newvec);CHKERRQ(ierr); 119 ierr = VecSetBlockSize(*newvec,dd->w);CHKERRQ(ierr); 120 ierr = VecSetSizes(*newvec,dd->M*dd->w,PETSC_DETERMINE);CHKERRQ(ierr); 121 ierr = VecSetType(*newvec,VECSEQ);CHKERRQ(ierr); 122 ierr = ISCreateGeneral(PETSC_COMM_SELF,dd->w*dd->M,indices,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 123 } else if (dir == DMDA_X) { 124 ierr = PetscMalloc(dd->w*dd->N*sizeof(PetscInt),&indices);CHKERRQ(ierr); 125 indices[0] = dd->w*gp; 126 for (j=1; j<dd->w; j++) indices[j] = indices[j-1] + 1; 127 for (i=1; i<dd->N; i++) { 128 indices[i*dd->w] = indices[i*dd->w-1] + dd->w*dd->M - dd->w + 1; 129 for (j=1; j<dd->w; j++) indices[i*dd->w + j] = indices[i*dd->w + j - 1] + 1; 130 } 131 ierr = AOApplicationToPetsc(ao,dd->w*dd->N,indices);CHKERRQ(ierr); 132 ierr = VecCreate(PETSC_COMM_SELF,newvec);CHKERRQ(ierr); 133 ierr = VecSetBlockSize(*newvec,dd->w);CHKERRQ(ierr); 134 ierr = VecSetSizes(*newvec,dd->N*dd->w,PETSC_DETERMINE);CHKERRQ(ierr); 135 ierr = VecSetType(*newvec,VECSEQ);CHKERRQ(ierr); 136 ierr = ISCreateGeneral(PETSC_COMM_SELF,dd->w*dd->N,indices,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 137 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown DMDADirection"); 138 } else { 139 ierr = VecCreateSeq(PETSC_COMM_SELF,0,newvec);CHKERRQ(ierr); 140 ierr = ISCreateGeneral(PETSC_COMM_SELF,0,0,PETSC_COPY_VALUES,&is);CHKERRQ(ierr); 141 } 142 ierr = DMGetGlobalVector(da,&vec);CHKERRQ(ierr); 143 ierr = VecScatterCreate(vec,is,*newvec,NULL,scatter);CHKERRQ(ierr); 144 ierr = DMRestoreGlobalVector(da,&vec);CHKERRQ(ierr); 145 ierr = ISDestroy(&is);CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNCT__ 150 #define __FUNCT__ "DMDAGetProcessorSubset" 151 /*@C 152 DMDAGetProcessorSubset - Returns a communicator consisting only of the 153 processors in a DMDA that own a particular global x, y, or z grid point 154 (corresponding to a logical plane in a 3D grid or a line in a 2D grid). 155 156 Collective on DMDA 157 158 Input Parameters: 159 + da - the distributed array 160 . dir - Cartesian direction, either DMDA_X, DMDA_Y, or DMDA_Z 161 - gp - global grid point number in this direction 162 163 Output Parameters: 164 . comm - new communicator 165 166 Level: advanced 167 168 Notes: 169 All processors that share the DMDA must call this with the same gp value 170 171 This routine is particularly useful to compute boundary conditions 172 or other application-specific calculations that require manipulating 173 sets of data throughout a logical plane of grid points. 174 175 .keywords: distributed array, get, processor subset 176 @*/ 177 PetscErrorCode DMDAGetProcessorSubset(DM da,DMDADirection dir,PetscInt gp,MPI_Comm *comm) 178 { 179 MPI_Group group,subgroup; 180 PetscErrorCode ierr; 181 PetscInt i,ict,flag,*owners,xs,xm,ys,ym,zs,zm; 182 PetscMPIInt size,*ranks = NULL; 183 DM_DA *dd = (DM_DA*)da->data; 184 185 PetscFunctionBegin; 186 PetscValidHeaderSpecific(da,DM_CLASSID,1); 187 flag = 0; 188 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 189 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);CHKERRQ(ierr); 190 if (dir == DMDA_Z) { 191 if (dd->dim < 3) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"DMDA_Z invalid for DMDA dim < 3"); 192 if (gp < 0 || gp > dd->P) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"invalid grid point"); 193 if (gp >= zs && gp < zs+zm) flag = 1; 194 } else if (dir == DMDA_Y) { 195 if (dd->dim == 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"DMDA_Y invalid for DMDA dim = 1"); 196 if (gp < 0 || gp > dd->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"invalid grid point"); 197 if (gp >= ys && gp < ys+ym) flag = 1; 198 } else if (dir == DMDA_X) { 199 if (gp < 0 || gp > dd->M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"invalid grid point"); 200 if (gp >= xs && gp < xs+xm) flag = 1; 201 } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"Invalid direction"); 202 203 ierr = PetscMalloc2(size,PetscInt,&owners,size,PetscMPIInt,&ranks);CHKERRQ(ierr); 204 ierr = MPI_Allgather(&flag,1,MPIU_INT,owners,1,MPIU_INT,PetscObjectComm((PetscObject)da));CHKERRQ(ierr); 205 ict = 0; 206 ierr = PetscInfo2(da,"DMDAGetProcessorSubset: dim=%D, direction=%d, procs: ",dd->dim,(int)dir);CHKERRQ(ierr); 207 for (i=0; i<size; i++) { 208 if (owners[i]) { 209 ranks[ict] = i; ict++; 210 ierr = PetscInfo1(da,"%D ",i);CHKERRQ(ierr); 211 } 212 } 213 ierr = PetscInfo(da,"\n");CHKERRQ(ierr); 214 ierr = MPI_Comm_group(PetscObjectComm((PetscObject)da),&group);CHKERRQ(ierr); 215 ierr = MPI_Group_incl(group,ict,ranks,&subgroup);CHKERRQ(ierr); 216 ierr = MPI_Comm_create(PetscObjectComm((PetscObject)da),subgroup,comm);CHKERRQ(ierr); 217 ierr = MPI_Group_free(&subgroup);CHKERRQ(ierr); 218 ierr = MPI_Group_free(&group);CHKERRQ(ierr); 219 ierr = PetscFree2(owners,ranks);CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "DMDAGetProcessorSubsets" 225 /*@C 226 DMDAGetProcessorSubsets - Returns communicators consisting only of the 227 processors in a DMDA adjacent in a particular dimension, 228 corresponding to a logical plane in a 3D grid or a line in a 2D grid. 229 230 Collective on DMDA 231 232 Input Parameters: 233 + da - the distributed array 234 - dir - Cartesian direction, either DMDA_X, DMDA_Y, or DMDA_Z 235 236 Output Parameters: 237 . subcomm - new communicator 238 239 Level: advanced 240 241 Notes: 242 This routine is useful for distributing one-dimensional data in a tensor product grid. 243 244 .keywords: distributed array, get, processor subset 245 @*/ 246 PetscErrorCode DMDAGetProcessorSubsets(DM da, DMDADirection dir, MPI_Comm *subcomm) 247 { 248 MPI_Comm comm; 249 MPI_Group group, subgroup; 250 PetscInt subgroupSize = 0; 251 PetscInt *firstPoints; 252 PetscMPIInt size, *subgroupRanks = NULL; 253 PetscInt xs, xm, ys, ym, zs, zm, firstPoint, p; 254 PetscErrorCode ierr; 255 DM_DA *dd = (DM_DA*)da->data; 256 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(da, DM_CLASSID, 1); 259 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 260 ierr = DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);CHKERRQ(ierr); 261 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 262 if (dir == DMDA_Z) { 263 if (dd->dim < 3) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"DMDA_Z invalid for DMDA dim < 3"); 264 firstPoint = zs; 265 } else if (dir == DMDA_Y) { 266 if (dd->dim == 1) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"DMDA_Y invalid for DMDA dim = 1"); 267 firstPoint = ys; 268 } else if (dir == DMDA_X) { 269 firstPoint = xs; 270 } else SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid direction"); 271 272 ierr = PetscMalloc2(size, PetscInt, &firstPoints, size, PetscMPIInt, &subgroupRanks);CHKERRQ(ierr); 273 ierr = MPI_Allgather(&firstPoint, 1, MPIU_INT, firstPoints, 1, MPIU_INT, comm);CHKERRQ(ierr); 274 ierr = PetscInfo2(da,"DMDAGetProcessorSubset: dim=%D, direction=%d, procs: ",dd->dim,(int)dir);CHKERRQ(ierr); 275 for (p = 0; p < size; ++p) { 276 if (firstPoints[p] == firstPoint) { 277 subgroupRanks[subgroupSize++] = p; 278 ierr = PetscInfo1(da, "%D ", p);CHKERRQ(ierr); 279 } 280 } 281 ierr = PetscInfo(da, "\n");CHKERRQ(ierr); 282 ierr = MPI_Comm_group(comm, &group);CHKERRQ(ierr); 283 ierr = MPI_Group_incl(group, subgroupSize, subgroupRanks, &subgroup);CHKERRQ(ierr); 284 ierr = MPI_Comm_create(comm, subgroup, subcomm);CHKERRQ(ierr); 285 ierr = MPI_Group_free(&subgroup);CHKERRQ(ierr); 286 ierr = MPI_Group_free(&group);CHKERRQ(ierr); 287 ierr = PetscFree2(firstPoints, subgroupRanks);CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290