1 2 static char help[] = "Tests DMSwarm\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 #include <petscdmswarm.h> 7 8 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM,PetscErrorCode (*)(DM,void*,PetscInt*,PetscInt**),size_t,void*,PetscInt*); 9 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM,PetscInt*); 10 11 PetscErrorCode ex1_1(void) 12 { 13 DM dms; 14 PetscErrorCode ierr; 15 Vec x; 16 PetscMPIInt rank,size; 17 PetscInt p; 18 19 PetscFunctionBegin; 20 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 21 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 22 if ((size > 1) && (size != 4)) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Must be run wuth 4 MPI ranks"); 23 24 ierr = DMCreate(PETSC_COMM_WORLD,&dms);CHKERRQ(ierr); 25 ierr = DMSetType(dms,DMSWARM);CHKERRQ(ierr); 26 27 ierr = DMSwarmInitializeFieldRegister(dms);CHKERRQ(ierr); 28 ierr = DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL);CHKERRQ(ierr); 29 ierr = DMSwarmRegisterPetscDatatypeField(dms,"strain",1,PETSC_REAL);CHKERRQ(ierr); 30 ierr = DMSwarmFinalizeFieldRegister(dms);CHKERRQ(ierr); 31 ierr = DMSwarmSetLocalSizes(dms,5+rank,4);CHKERRQ(ierr); 32 ierr = DMView(dms,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 33 34 { 35 PetscReal *array; 36 ierr = DMSwarmGetField(dms,"viscosity",NULL,NULL,(void**)&array);CHKERRQ(ierr); 37 for (p=0; p<5+rank; p++) { 38 array[p] = 11.1 + p*0.1 + rank*100.0; 39 } 40 ierr = DMSwarmRestoreField(dms,"viscosity",NULL,NULL,(void**)&array);CHKERRQ(ierr); 41 } 42 43 { 44 PetscReal *array; 45 ierr = DMSwarmGetField(dms,"strain",NULL,NULL,(void**)&array);CHKERRQ(ierr); 46 for (p=0; p<5+rank; p++) { 47 array[p] = 2.0e-2 + p*0.001 + rank*1.0; 48 } 49 ierr = DMSwarmRestoreField(dms,"strain",NULL,NULL,(void**)&array);CHKERRQ(ierr); 50 } 51 52 ierr = DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr); 53 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr); 54 55 ierr = DMSwarmVectorDefineField(dms,"strain");CHKERRQ(ierr); 56 ierr = DMCreateGlobalVector(dms,&x);CHKERRQ(ierr); 57 ierr = VecDestroy(&x);CHKERRQ(ierr); 58 59 { 60 PetscInt *rankval; 61 PetscInt npoints[2],npoints_orig[2]; 62 63 ierr = DMSwarmGetLocalSize(dms,&npoints_orig[0]);CHKERRQ(ierr); 64 ierr = DMSwarmGetSize(dms,&npoints_orig[1]);CHKERRQ(ierr); 65 ierr = DMSwarmGetField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 66 if ((rank == 0) && (size > 1)) { 67 rankval[0] = 1; 68 rankval[3] = 1; 69 } 70 if (rank == 3) { 71 rankval[2] = 1; 72 } 73 ierr = DMSwarmRestoreField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 74 ierr = DMSwarmMigrate(dms,PETSC_TRUE);CHKERRQ(ierr); 75 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr); 76 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr); 77 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 78 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%D,%D) after(%D,%D)\n",rank,npoints_orig[0],npoints_orig[1],npoints[0],npoints[1]);CHKERRQ(ierr); 79 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 80 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 81 } 82 { 83 ierr = DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr); 84 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 85 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr); 86 } 87 { 88 ierr = DMSwarmCreateGlobalVectorFromField(dms,"strain",&x);CHKERRQ(ierr); 89 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 90 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"strain",&x);CHKERRQ(ierr); 91 } 92 93 ierr = DMDestroy(&dms);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 97 PetscErrorCode ex1_2(void) 98 { 99 DM dms; 100 PetscErrorCode ierr; 101 Vec x; 102 PetscMPIInt rank,size; 103 PetscInt p; 104 105 PetscFunctionBegin; 106 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 107 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 108 ierr = DMCreate(PETSC_COMM_WORLD,&dms);CHKERRQ(ierr); 109 ierr = DMSetType(dms,DMSWARM);CHKERRQ(ierr); 110 ierr = DMSwarmInitializeFieldRegister(dms);CHKERRQ(ierr); 111 112 ierr = DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL);CHKERRQ(ierr); 113 ierr = DMSwarmRegisterPetscDatatypeField(dms,"strain",1,PETSC_REAL);CHKERRQ(ierr); 114 ierr = DMSwarmFinalizeFieldRegister(dms);CHKERRQ(ierr); 115 ierr = DMSwarmSetLocalSizes(dms,5+rank,4);CHKERRQ(ierr); 116 ierr = DMView(dms,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 117 { 118 PetscReal *array; 119 ierr = DMSwarmGetField(dms,"viscosity",NULL,NULL,(void**)&array);CHKERRQ(ierr); 120 for (p=0; p<5+rank; p++) { 121 array[p] = 11.1 + p*0.1 + rank*100.0; 122 } 123 ierr = DMSwarmRestoreField(dms,"viscosity",NULL,NULL,(void**)&array);CHKERRQ(ierr); 124 } 125 { 126 PetscReal *array; 127 ierr = DMSwarmGetField(dms,"strain",NULL,NULL,(void**)&array);CHKERRQ(ierr); 128 for (p=0; p<5+rank; p++) { 129 array[p] = 2.0e-2 + p*0.001 + rank*1.0; 130 } 131 ierr = DMSwarmRestoreField(dms,"strain",NULL,NULL,(void**)&array);CHKERRQ(ierr); 132 } 133 { 134 PetscInt *rankval; 135 PetscInt npoints[2],npoints_orig[2]; 136 137 ierr = DMSwarmGetLocalSize(dms,&npoints_orig[0]);CHKERRQ(ierr); 138 ierr = DMSwarmGetSize(dms,&npoints_orig[1]);CHKERRQ(ierr); 139 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 140 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%D,%D)\n",rank,npoints_orig[0],npoints_orig[1]);CHKERRQ(ierr); 141 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 142 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 143 144 ierr = DMSwarmGetField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 145 146 if (rank == 1) { 147 rankval[0] = -1; 148 } 149 if (rank == 2) { 150 rankval[1] = -1; 151 } 152 if (rank == 3) { 153 rankval[3] = -1; 154 rankval[4] = -1; 155 } 156 ierr = DMSwarmRestoreField(dms,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 157 ierr = DMSwarmCollectViewCreate(dms);CHKERRQ(ierr); 158 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr); 159 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr); 160 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 161 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%D,%D)\n",rank,npoints[0],npoints[1]);CHKERRQ(ierr); 162 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 163 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 164 165 ierr = DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr); 166 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 167 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr); 168 169 ierr = DMSwarmCollectViewDestroy(dms);CHKERRQ(ierr); 170 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr); 171 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr); 172 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 173 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after_v(%D,%D)\n",rank,npoints[0],npoints[1]);CHKERRQ(ierr); 174 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 175 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 176 177 ierr = DMSwarmCreateGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr); 178 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 179 ierr = DMSwarmDestroyGlobalVectorFromField(dms,"viscosity",&x);CHKERRQ(ierr); 180 } 181 ierr = DMDestroy(&dms);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 /* 186 splot "c-rank0.gp","c-rank1.gp","c-rank2.gp","c-rank3.gp" 187 */ 188 PetscErrorCode ex1_3(void) 189 { 190 DM dms; 191 PetscErrorCode ierr; 192 PetscMPIInt rank,size; 193 PetscInt is,js,ni,nj,overlap; 194 DM dmcell; 195 196 PetscFunctionBegin; 197 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 198 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 199 overlap = 2; 200 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,13,13,PETSC_DECIDE,PETSC_DECIDE,1,overlap,NULL,NULL,&dmcell);CHKERRQ(ierr); 201 ierr = DMSetFromOptions(dmcell);CHKERRQ(ierr); 202 ierr = DMSetUp(dmcell);CHKERRQ(ierr); 203 ierr = DMDASetUniformCoordinates(dmcell,-1.0,1.0,-1.0,1.0,-1.0,1.0);CHKERRQ(ierr); 204 ierr = DMDAGetCorners(dmcell,&is,&js,NULL,&ni,&nj,NULL);CHKERRQ(ierr); 205 ierr = DMCreate(PETSC_COMM_WORLD,&dms);CHKERRQ(ierr); 206 ierr = DMSetType(dms,DMSWARM);CHKERRQ(ierr); 207 ierr = DMSwarmSetCellDM(dms,dmcell);CHKERRQ(ierr); 208 209 /* load in data types */ 210 ierr = DMSwarmInitializeFieldRegister(dms);CHKERRQ(ierr); 211 ierr = DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL);CHKERRQ(ierr); 212 ierr = DMSwarmRegisterPetscDatatypeField(dms,"coorx",1,PETSC_REAL);CHKERRQ(ierr); 213 ierr = DMSwarmRegisterPetscDatatypeField(dms,"coory",1,PETSC_REAL);CHKERRQ(ierr); 214 ierr = DMSwarmFinalizeFieldRegister(dms);CHKERRQ(ierr); 215 ierr = DMSwarmSetLocalSizes(dms,ni*nj*4,4);CHKERRQ(ierr); 216 ierr = DMView(dms,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 217 218 /* set values within the swarm */ 219 { 220 PetscReal *array_x,*array_y; 221 PetscInt npoints,i,j,cnt; 222 DMDACoor2d **LA_coor; 223 Vec coor; 224 DM dmcellcdm; 225 226 ierr = DMGetCoordinateDM(dmcell,&dmcellcdm);CHKERRQ(ierr); 227 ierr = DMGetCoordinates(dmcell,&coor);CHKERRQ(ierr); 228 ierr = DMDAVecGetArray(dmcellcdm,coor,&LA_coor);CHKERRQ(ierr); 229 ierr = DMSwarmGetLocalSize(dms,&npoints);CHKERRQ(ierr); 230 ierr = DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 231 ierr = DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 232 cnt = 0; 233 for (j=js; j<js+nj; j++) { 234 for (i=is; i<is+ni; i++) { 235 PetscReal xp,yp; 236 xp = PetscRealPart(LA_coor[j][i].x); 237 yp = PetscRealPart(LA_coor[j][i].y); 238 array_x[4*cnt+0] = xp - 0.05; if (array_x[4*cnt+0] < -1.0) { array_x[4*cnt+0] = -1.0+1.0e-12; } 239 array_x[4*cnt+1] = xp + 0.05; if (array_x[4*cnt+1] > 1.0) { array_x[4*cnt+1] = 1.0-1.0e-12; } 240 array_x[4*cnt+2] = xp - 0.05; if (array_x[4*cnt+2] < -1.0) { array_x[4*cnt+2] = -1.0+1.0e-12; } 241 array_x[4*cnt+3] = xp + 0.05; if (array_x[4*cnt+3] > 1.0) { array_x[4*cnt+3] = 1.0-1.0e-12; } 242 243 array_y[4*cnt+0] = yp - 0.05; if (array_y[4*cnt+0] < -1.0) { array_y[4*cnt+0] = -1.0+1.0e-12; } 244 array_y[4*cnt+1] = yp - 0.05; if (array_y[4*cnt+1] < -1.0) { array_y[4*cnt+1] = -1.0+1.0e-12; } 245 array_y[4*cnt+2] = yp + 0.05; if (array_y[4*cnt+2] > 1.0) { array_y[4*cnt+2] = 1.0-1.0e-12; } 246 array_y[4*cnt+3] = yp + 0.05; if (array_y[4*cnt+3] > 1.0) { array_y[4*cnt+3] = 1.0-1.0e-12; } 247 cnt++; 248 } 249 } 250 ierr = DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 251 ierr = DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 252 ierr = DMDAVecRestoreArray(dmcellcdm,coor,&LA_coor);CHKERRQ(ierr); 253 } 254 { 255 PetscInt npoints[2],npoints_orig[2],ng; 256 257 ierr = DMSwarmGetLocalSize(dms,&npoints_orig[0]);CHKERRQ(ierr); 258 ierr = DMSwarmGetSize(dms,&npoints_orig[1]);CHKERRQ(ierr); 259 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 260 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%D,%D)\n",rank,npoints_orig[0],npoints_orig[1]);CHKERRQ(ierr); 261 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 262 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 263 ierr = DMSwarmCollect_DMDABoundingBox(dms,&ng);CHKERRQ(ierr); 264 265 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr); 266 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr); 267 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 268 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%D,%D)\n",rank,npoints[0],npoints[1]);CHKERRQ(ierr); 269 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 270 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 271 } 272 { 273 PetscReal *array_x,*array_y; 274 PetscInt npoints,p; 275 FILE *fp = NULL; 276 char name[PETSC_MAX_PATH_LEN]; 277 278 ierr = PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"c-rank%d.gp",rank);CHKERRQ(ierr); 279 fp = fopen(name,"w"); 280 if (!fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file %s",name); 281 ierr = DMSwarmGetLocalSize(dms,&npoints);CHKERRQ(ierr); 282 ierr = DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 283 ierr = DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 284 for (p=0; p<npoints; p++) { 285 fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array_x[p],array_y[p],(double)rank); 286 } 287 ierr = DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 288 ierr = DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 289 fclose(fp); 290 } 291 ierr = DMDestroy(&dmcell);CHKERRQ(ierr); 292 ierr = DMDestroy(&dms);CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 typedef struct { 297 PetscReal cx[2]; 298 PetscReal radius; 299 } CollectZoneCtx; 300 301 PetscErrorCode collect_zone(DM dm,void *ctx,PetscInt *nfound,PetscInt **foundlist) 302 { 303 CollectZoneCtx *zone = (CollectZoneCtx*)ctx; 304 PetscInt p,npoints; 305 PetscReal *array_x,*array_y,r2; 306 PetscInt p2collect,*plist; 307 PetscMPIInt rank; 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 312 ierr = DMSwarmGetLocalSize(dm,&npoints);CHKERRQ(ierr); 313 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 314 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 315 316 r2 = zone->radius * zone->radius; 317 p2collect = 0; 318 for (p=0; p<npoints; p++) { 319 PetscReal sep2; 320 321 sep2 = (array_x[p] - zone->cx[0])*(array_x[p] - zone->cx[0]); 322 sep2 += (array_y[p] - zone->cx[1])*(array_y[p] - zone->cx[1]); 323 if (sep2 < r2) { 324 p2collect++; 325 } 326 } 327 328 ierr = PetscMalloc1(p2collect+1,&plist);CHKERRQ(ierr); 329 p2collect = 0; 330 for (p=0; p<npoints; p++) { 331 PetscReal sep2; 332 333 sep2 = (array_x[p] - zone->cx[0])*(array_x[p] - zone->cx[0]); 334 sep2 += (array_y[p] - zone->cx[1])*(array_y[p] - zone->cx[1]); 335 if (sep2 < r2) { 336 plist[p2collect] = p; 337 p2collect++; 338 } 339 } 340 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 341 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 342 343 *nfound = p2collect; 344 *foundlist = plist; 345 PetscFunctionReturn(0); 346 } 347 348 PetscErrorCode ex1_4(void) 349 { 350 DM dms; 351 PetscErrorCode ierr; 352 PetscMPIInt rank,size; 353 PetscInt is,js,ni,nj,overlap,nn; 354 DM dmcell; 355 CollectZoneCtx *zone; 356 PetscReal dx; 357 358 PetscFunctionBegin; 359 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 360 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 361 nn = 101; 362 dx = 2.0/ (PetscReal)(nn-1); 363 overlap = 0; 364 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nn,nn,PETSC_DECIDE,PETSC_DECIDE,1,overlap,NULL,NULL,&dmcell);CHKERRQ(ierr); 365 ierr = DMSetFromOptions(dmcell);CHKERRQ(ierr); 366 ierr = DMSetUp(dmcell);CHKERRQ(ierr); 367 ierr = DMDASetUniformCoordinates(dmcell,-1.0,1.0,-1.0,1.0,-1.0,1.0);CHKERRQ(ierr); 368 ierr = DMDAGetCorners(dmcell,&is,&js,NULL,&ni,&nj,NULL);CHKERRQ(ierr); 369 ierr = DMCreate(PETSC_COMM_WORLD,&dms);CHKERRQ(ierr); 370 ierr = DMSetType(dms,DMSWARM);CHKERRQ(ierr); 371 372 /* load in data types */ 373 ierr = DMSwarmInitializeFieldRegister(dms);CHKERRQ(ierr); 374 ierr = DMSwarmRegisterPetscDatatypeField(dms,"viscosity",1,PETSC_REAL);CHKERRQ(ierr); 375 ierr = DMSwarmRegisterPetscDatatypeField(dms,"coorx",1,PETSC_REAL);CHKERRQ(ierr); 376 ierr = DMSwarmRegisterPetscDatatypeField(dms,"coory",1,PETSC_REAL);CHKERRQ(ierr); 377 ierr = DMSwarmFinalizeFieldRegister(dms);CHKERRQ(ierr); 378 ierr = DMSwarmSetLocalSizes(dms,ni*nj*4,4);CHKERRQ(ierr); 379 ierr = DMView(dms,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 380 381 /* set values within the swarm */ 382 { 383 PetscReal *array_x,*array_y; 384 PetscInt npoints,i,j,cnt; 385 DMDACoor2d **LA_coor; 386 Vec coor; 387 DM dmcellcdm; 388 389 ierr = DMGetCoordinateDM(dmcell,&dmcellcdm);CHKERRQ(ierr); 390 ierr = DMGetCoordinates(dmcell,&coor);CHKERRQ(ierr); 391 ierr = DMDAVecGetArray(dmcellcdm,coor,&LA_coor);CHKERRQ(ierr); 392 ierr = DMSwarmGetLocalSize(dms,&npoints);CHKERRQ(ierr); 393 ierr = DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 394 ierr = DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 395 cnt = 0; 396 for (j=js; j<js+nj; j++) { 397 for (i=is; i<is+ni; i++) { 398 PetscReal xp,yp; 399 400 xp = PetscRealPart(LA_coor[j][i].x); 401 yp = PetscRealPart(LA_coor[j][i].y); 402 array_x[4*cnt+0] = xp - dx*0.1; /*if (array_x[4*cnt+0] < -1.0) { array_x[4*cnt+0] = -1.0+1.0e-12; }*/ 403 array_x[4*cnt+1] = xp + dx*0.1; /*if (array_x[4*cnt+1] > 1.0) { array_x[4*cnt+1] = 1.0-1.0e-12; }*/ 404 array_x[4*cnt+2] = xp - dx*0.1; /*if (array_x[4*cnt+2] < -1.0) { array_x[4*cnt+2] = -1.0+1.0e-12; }*/ 405 array_x[4*cnt+3] = xp + dx*0.1; /*if (array_x[4*cnt+3] > 1.0) { array_x[4*cnt+3] = 1.0-1.0e-12; }*/ 406 array_y[4*cnt+0] = yp - dx*0.1; /*if (array_y[4*cnt+0] < -1.0) { array_y[4*cnt+0] = -1.0+1.0e-12; }*/ 407 array_y[4*cnt+1] = yp - dx*0.1; /*if (array_y[4*cnt+1] < -1.0) { array_y[4*cnt+1] = -1.0+1.0e-12; }*/ 408 array_y[4*cnt+2] = yp + dx*0.1; /*if (array_y[4*cnt+2] > 1.0) { array_y[4*cnt+2] = 1.0-1.0e-12; }*/ 409 array_y[4*cnt+3] = yp + dx*0.1; /*if (array_y[4*cnt+3] > 1.0) { array_y[4*cnt+3] = 1.0-1.0e-12; }*/ 410 cnt++; 411 } 412 } 413 ierr = DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 414 ierr = DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 415 ierr = DMDAVecRestoreArray(dmcellcdm,coor,&LA_coor);CHKERRQ(ierr); 416 } 417 ierr = PetscMalloc1(1,&zone);CHKERRQ(ierr); 418 if (size == 4) { 419 if (rank == 0) { 420 zone->cx[0] = 0.5; 421 zone->cx[1] = 0.5; 422 zone->radius = 0.3; 423 } 424 if (rank == 1) { 425 zone->cx[0] = -0.5; 426 zone->cx[1] = 0.5; 427 zone->radius = 0.25; 428 } 429 if (rank == 2) { 430 zone->cx[0] = 0.5; 431 zone->cx[1] = -0.5; 432 zone->radius = 0.2; 433 } 434 if (rank == 3) { 435 zone->cx[0] = -0.5; 436 zone->cx[1] = -0.5; 437 zone->radius = 0.1; 438 } 439 } else { 440 if (rank == 0) { 441 zone->cx[0] = 0.5; 442 zone->cx[1] = 0.5; 443 zone->radius = 0.8; 444 } else { 445 zone->cx[0] = 10.0; 446 zone->cx[1] = 10.0; 447 zone->radius = 0.0; 448 } 449 } 450 { 451 PetscInt npoints[2],npoints_orig[2],ng; 452 453 ierr = DMSwarmGetLocalSize(dms,&npoints_orig[0]);CHKERRQ(ierr); 454 ierr = DMSwarmGetSize(dms,&npoints_orig[1]);CHKERRQ(ierr); 455 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 456 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] before(%D,%D)\n",rank,npoints_orig[0],npoints_orig[1]);CHKERRQ(ierr); 457 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 458 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 459 ierr = DMSwarmCollect_General(dms,collect_zone,sizeof(CollectZoneCtx),zone,&ng);CHKERRQ(ierr); 460 ierr = DMSwarmGetLocalSize(dms,&npoints[0]);CHKERRQ(ierr); 461 ierr = DMSwarmGetSize(dms,&npoints[1]);CHKERRQ(ierr); 462 ierr = PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 463 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"rank[%d] after(%D,%D)\n",rank,npoints[0],npoints[1]);CHKERRQ(ierr); 464 ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 465 ierr = PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 466 } 467 { 468 PetscReal *array_x,*array_y; 469 PetscInt npoints,p; 470 FILE *fp = NULL; 471 char name[PETSC_MAX_PATH_LEN]; 472 473 ierr = PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"c-rank%d.gp",rank);CHKERRQ(ierr); 474 fp = fopen(name,"w"); 475 if (!fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file %s",name); 476 ierr = DMSwarmGetLocalSize(dms,&npoints);CHKERRQ(ierr); 477 ierr = DMSwarmGetField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 478 ierr = DMSwarmGetField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 479 for (p=0; p<npoints; p++) { 480 fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array_x[p],array_y[p],(double)rank); 481 } 482 ierr = DMSwarmRestoreField(dms,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 483 ierr = DMSwarmRestoreField(dms,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 484 fclose(fp); 485 } 486 ierr = DMDestroy(&dmcell);CHKERRQ(ierr); 487 ierr = DMDestroy(&dms);CHKERRQ(ierr); 488 ierr = PetscFree(zone);CHKERRQ(ierr); 489 PetscFunctionReturn(0); 490 } 491 492 int main(int argc,char **argv) 493 { 494 PetscErrorCode ierr; 495 PetscInt test_mode = 4; 496 497 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 498 ierr = PetscOptionsGetInt(NULL,NULL,"-test_mode",&test_mode,NULL);CHKERRQ(ierr); 499 if (test_mode == 1) { 500 ierr = ex1_1();CHKERRQ(ierr); 501 } else if (test_mode == 2) { 502 ierr = ex1_2();CHKERRQ(ierr); 503 } else if (test_mode == 3) { 504 ierr = ex1_3();CHKERRQ(ierr); 505 } else if (test_mode == 4) { 506 ierr = ex1_4();CHKERRQ(ierr); 507 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown test_mode value, should be 1,2,3,4"); 508 ierr = PetscFinalize(); 509 return ierr; 510 } 511 512 /*TEST 513 514 build: 515 requires: !complex double 516 517 test: 518 args: -test_mode 1 519 filter: grep -v atomic 520 filter_output: grep -v atomic 521 522 test: 523 suffix: 2 524 args: -test_mode 2 525 filter: grep -v atomic 526 filter_output: grep -v atomic 527 528 test: 529 suffix: 3 530 args: -test_mode 3 531 filter: grep -v atomic 532 filter_output: grep -v atomic 533 TODO: broken 534 535 test: 536 suffix: 4 537 args: -test_mode 4 538 filter: grep -v atomic 539 filter_output: grep -v atomic 540 541 test: 542 suffix: 5 543 nsize: 4 544 args: -test_mode 1 545 filter: grep -v atomic 546 filter_output: grep -v atomic 547 548 test: 549 suffix: 6 550 nsize: 4 551 args: -test_mode 2 552 filter: grep -v atomic 553 filter_output: grep -v atomic 554 555 test: 556 suffix: 7 557 nsize: 4 558 args: -test_mode 3 559 filter: grep -v atomic 560 filter_output: grep -v atomic 561 TODO: broken 562 563 test: 564 suffix: 8 565 nsize: 4 566 args: -test_mode 4 567 filter: grep -v atomic 568 filter_output: grep -v atomic 569 570 TEST*/ 571