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