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