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