1 2 static char help[] = "Tests DMSwarm with DMShell\n\n"; 3 4 #include <petscsf.h> 5 #include <petscdm.h> 6 #include <petscdmshell.h> 7 #include <petscdmda.h> 8 #include <petscdmswarm.h> 9 #include <petsc/private/dmimpl.h> 10 11 PetscErrorCode _DMLocatePoints_DMDARegular_IS(DM dm,Vec pos,IS *iscell) 12 { 13 PetscInt p,n,bs,npoints,si,sj,milocal,mjlocal,mx,my; 14 DM dmregular; 15 PetscInt *cellidx; 16 const PetscScalar *coor; 17 PetscReal dx,dy; 18 PetscMPIInt rank; 19 20 PetscFunctionBegin; 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 22 PetscCall(VecGetLocalSize(pos,&n)); 23 PetscCall(VecGetBlockSize(pos,&bs)); 24 npoints = n/bs; 25 26 PetscCall(PetscMalloc1(npoints,&cellidx)); 27 PetscCall(DMGetApplicationContext(dm,&dmregular)); 28 PetscCall(DMDAGetCorners(dmregular,&si,&sj,NULL,&milocal,&mjlocal,NULL)); 29 PetscCall(DMDAGetInfo(dmregular,NULL,&mx,&my,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 30 31 dx = 2.0/((PetscReal)mx); 32 dy = 2.0/((PetscReal)my); 33 34 PetscCall(VecGetArrayRead(pos,&coor)); 35 for (p=0; p<npoints; p++) { 36 PetscReal coorx,coory; 37 PetscInt mi,mj; 38 39 coorx = PetscRealPart(coor[2*p]); 40 coory = PetscRealPart(coor[2*p+1]); 41 42 mi = (PetscInt)( (coorx - (-1.0))/dx); 43 mj = (PetscInt)( (coory - (-1.0))/dy); 44 45 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 46 47 if ((mj >= sj) && (mj < sj + mjlocal)) { 48 if ((mi >= si) && (mi < si + milocal)) { 49 cellidx[p] = (mi-si) + (mj-sj) * milocal; 50 } 51 } 52 if (coorx < -1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 53 if (coorx > 1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 54 if (coory < -1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 55 if (coory > 1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 56 } 57 PetscCall(VecRestoreArrayRead(pos,&coor)); 58 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell)); 59 PetscFunctionReturn(0); 60 } 61 62 PetscErrorCode DMLocatePoints_DMDARegular(DM dm,Vec pos,DMPointLocationType ltype, PetscSF cellSF) 63 { 64 IS iscell; 65 PetscSFNode *cells; 66 PetscInt p,bs,npoints,nfound; 67 const PetscInt *boxCells; 68 69 PetscFunctionBegin; 70 PetscCall(_DMLocatePoints_DMDARegular_IS(dm,pos,&iscell)); 71 PetscCall(VecGetLocalSize(pos,&npoints)); 72 PetscCall(VecGetBlockSize(pos,&bs)); 73 npoints = npoints / bs; 74 75 PetscCall(PetscMalloc1(npoints, &cells)); 76 PetscCall(ISGetIndices(iscell, &boxCells)); 77 78 for (p=0; p<npoints; p++) { 79 cells[p].rank = 0; 80 cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND; 81 cells[p].index = boxCells[p]; 82 } 83 PetscCall(ISRestoreIndices(iscell, &boxCells)); 84 PetscCall(ISDestroy(&iscell)); 85 nfound = npoints; 86 PetscCall(PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER)); 87 PetscCall(ISDestroy(&iscell)); 88 PetscFunctionReturn(0); 89 } 90 91 PetscErrorCode DMGetNeighbors_DMDARegular(DM dm,PetscInt *nneighbors,const PetscMPIInt **neighbors) 92 { 93 DM dmregular; 94 95 PetscFunctionBegin; 96 PetscCall(DMGetApplicationContext(dm,&dmregular)); 97 PetscCall(DMGetNeighbors(dmregular,nneighbors,neighbors)); 98 PetscFunctionReturn(0); 99 } 100 101 PetscErrorCode SwarmViewGP(DM dms,const char prefix[]) 102 { 103 PetscReal *array; 104 PetscInt *iarray; 105 PetscInt npoints,p,bs; 106 FILE *fp; 107 char name[PETSC_MAX_PATH_LEN]; 108 PetscMPIInt rank; 109 110 PetscFunctionBegin; 111 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 112 PetscCall(PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"%s-rank%d.gp",prefix,rank)); 113 fp = fopen(name,"w"); 114 PetscCheck(fp,PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file %s",name); 115 PetscCall(DMSwarmGetLocalSize(dms,&npoints)); 116 PetscCall(DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array)); 117 PetscCall(DMSwarmGetField(dms,"itag",NULL,NULL,(void**)&iarray)); 118 for (p=0; p<npoints; p++) { 119 fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array[2*p],array[2*p+1],(double)iarray[p]); 120 } 121 PetscCall(DMSwarmRestoreField(dms,"itag",NULL,NULL,(void**)&iarray)); 122 PetscCall(DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array)); 123 fclose(fp); 124 PetscFunctionReturn(0); 125 } 126 127 /* 128 Create a DMShell and attach a regularly spaced DMDA for point location 129 Override methods for point location 130 */ 131 PetscErrorCode ex3_1(void) 132 { 133 DM dms,dmcell,dmregular; 134 PetscMPIInt rank; 135 PetscInt p,bs,nlocal,overlap,mx,tk; 136 PetscReal dx; 137 PetscReal *array,dt; 138 PetscInt *iarray; 139 PetscRandom rand; 140 141 PetscFunctionBegin; 142 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 143 144 /* Create a regularly spaced DMDA */ 145 mx = 40; 146 overlap = 0; 147 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,mx,PETSC_DECIDE,PETSC_DECIDE,1,overlap,NULL,NULL,&dmregular)); 148 PetscCall(DMSetFromOptions(dmregular)); 149 PetscCall(DMSetUp(dmregular)); 150 151 dx = 2.0/((PetscReal)mx); 152 PetscCall(DMDASetUniformCoordinates(dmregular,-1.0+0.5*dx,1.0-0.5*dx,-1.0+0.5*dx,1.0-0.5*dx,-1.0,1.0)); 153 154 /* Create a DMShell for point location purposes */ 155 PetscCall(DMShellCreate(PETSC_COMM_WORLD,&dmcell)); 156 PetscCall(DMSetApplicationContext(dmcell,dmregular)); 157 dmcell->ops->locatepoints = DMLocatePoints_DMDARegular; 158 dmcell->ops->getneighbors = DMGetNeighbors_DMDARegular; 159 160 /* Create the swarm */ 161 PetscCall(DMCreate(PETSC_COMM_WORLD,&dms)); 162 PetscCall(DMSetType(dms,DMSWARM)); 163 PetscCall(DMSetDimension(dms,2)); 164 PetscCall(PetscObjectSetName((PetscObject) dms, "Particles")); 165 166 PetscCall(DMSwarmSetType(dms,DMSWARM_PIC)); 167 PetscCall(DMSwarmSetCellDM(dms,dmcell)); 168 169 PetscCall(DMSwarmRegisterPetscDatatypeField(dms,"itag",1,PETSC_INT)); 170 PetscCall(DMSwarmFinalizeFieldRegister(dms)); 171 { 172 PetscInt si,sj,milocal,mjlocal; 173 const PetscScalar *LA_coors; 174 Vec coors; 175 PetscInt cnt; 176 177 PetscCall(DMDAGetCorners(dmregular,&si,&sj,NULL,&milocal,&mjlocal,NULL)); 178 PetscCall(DMGetCoordinates(dmregular,&coors)); 179 PetscCall(VecGetArrayRead(coors,&LA_coors)); 180 PetscCall(DMSwarmSetLocalSizes(dms,milocal*mjlocal,4)); 181 PetscCall(DMSwarmGetLocalSize(dms,&nlocal)); 182 PetscCall(DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array)); 183 cnt = 0; 184 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand)); 185 PetscCall(PetscRandomSetInterval(rand,-1.0,1.0)); 186 for (p=0; p<nlocal; p++) { 187 PetscReal px,py,rx,ry,r2; 188 189 PetscCall(PetscRandomGetValueReal(rand,&rx)); 190 PetscCall(PetscRandomGetValueReal(rand,&ry)); 191 192 px = PetscRealPart(LA_coors[2*p+0]) + 0.1*rx*dx; 193 py = PetscRealPart(LA_coors[2*p+1]) + 0.1*ry*dx; 194 195 r2 = px*px + py*py; 196 if (r2 < 0.75*0.75) { 197 array[bs*cnt+0] = px; 198 array[bs*cnt+1] = py; 199 cnt++; 200 } 201 } 202 PetscCall(PetscRandomDestroy(&rand)); 203 PetscCall(DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array)); 204 PetscCall(VecRestoreArrayRead(coors,&LA_coors)); 205 PetscCall(DMSwarmSetLocalSizes(dms,cnt,4)); 206 207 PetscCall(DMSwarmGetLocalSize(dms,&nlocal)); 208 PetscCall(DMSwarmGetField(dms,"itag",&bs,NULL,(void**)&iarray)); 209 for (p=0; p<nlocal; p++) { 210 iarray[p] = (PetscInt)rank; 211 } 212 PetscCall(DMSwarmRestoreField(dms,"itag",&bs,NULL,(void**)&iarray)); 213 } 214 215 PetscCall(DMView(dms,PETSC_VIEWER_STDOUT_WORLD)); 216 PetscCall(SwarmViewGP(dms,"step0")); 217 218 dt = 0.1; 219 for (tk=1; tk<20; tk++) { 220 char prefix[PETSC_MAX_PATH_LEN]; 221 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Step %" PetscInt_FMT " \n",tk)); 222 /* push points */ 223 PetscCall(DMSwarmGetLocalSize(dms,&nlocal)); 224 PetscCall(DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array)); 225 for (p=0; p<nlocal; p++) { 226 PetscReal cx,cy,vx,vy; 227 228 cx = array[2*p]; 229 cy = array[2*p+1]; 230 vx = cy; 231 vy = -cx; 232 233 array[2*p ] += dt * vx; 234 array[2*p+1] += dt * vy; 235 } 236 PetscCall(DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array)); 237 238 /* migrate points */ 239 PetscCall(DMSwarmMigrate(dms,PETSC_TRUE)); 240 /* view points */ 241 PetscCall(PetscSNPrintf(prefix,PETSC_MAX_PATH_LEN-1,"step%" PetscInt_FMT,tk)); 242 /* should use the regular SwarmView() api, not one for a particular type */ 243 PetscCall(SwarmViewGP(dms,prefix)); 244 } 245 PetscCall(DMDestroy(&dmregular)); 246 PetscCall(DMDestroy(&dmcell)); 247 PetscCall(DMDestroy(&dms)); 248 PetscFunctionReturn(0); 249 } 250 251 int main(int argc,char **argv) 252 { 253 254 PetscFunctionBeginUser; 255 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 256 PetscCall(ex3_1()); 257 PetscCall(PetscFinalize()); 258 return 0; 259 } 260 261 /*TEST 262 263 build: 264 requires: double !complex 265 266 test: 267 filter: grep -v atomic 268 filter_output: grep -v atomic 269 TEST*/ 270