1 static char help[] = "Test point location in DA using DMSwarm\n";
2
3 #include <petscdmda.h>
4 #include <petscdmswarm.h>
5
DMSwarmPrint(DM sw)6 PetscErrorCode DMSwarmPrint(DM sw)
7 {
8 DMSwarmCellDM celldm;
9 PetscReal *array;
10 PetscInt *pidArray, *cidArray;
11 PetscInt Np, bs, Nfc;
12 PetscMPIInt rank;
13 const char **coordFields, *cellid;
14
15 PetscFunctionBeginUser;
16 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
17 PetscCall(DMSwarmGetLocalSize(sw, &Np));
18 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
19 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
20 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
21 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
22 PetscCall(DMSwarmGetField(sw, coordFields[0], &bs, NULL, (void **)&array));
23 PetscCall(DMSwarmGetField(sw, DMSwarmField_pid, &bs, NULL, (void **)&pidArray));
24 PetscCall(DMSwarmGetField(sw, cellid, &bs, NULL, (void **)&cidArray));
25 for (PetscInt p = 0; p < Np; ++p) {
26 const PetscReal th = PetscAtan2Real(array[2 * p + 1], array[2 * p]) / PETSC_PI;
27 const PetscReal r = PetscSqrtReal(array[2 * p + 1] * array[2 * p + 1] + array[2 * p] * array[2 * p]);
28 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] p %" PetscInt_FMT " (%+1.4f,%+1.4f) r=%1.2f th=%1.3f*pi cellid=%" PetscInt_FMT "\n", rank, pidArray[p], (double)array[2 * p], (double)array[2 * p + 1], (double)r, (double)th, cidArray[p]));
29 }
30 PetscCall(DMSwarmRestoreField(sw, coordFields[0], &bs, NULL, (void **)&array));
31 PetscCall(DMSwarmRestoreField(sw, DMSwarmField_pid, &bs, NULL, (void **)&pidArray));
32 PetscCall(DMSwarmRestoreField(sw, cellid, &bs, NULL, (void **)&pidArray));
33 PetscFunctionReturn(PETSC_SUCCESS);
34 }
35
main(int argc,char ** argv)36 int main(int argc, char **argv)
37 {
38 DM dm, sw;
39 PetscDataType dtype;
40 PetscReal *coords, r, dr;
41 PetscInt Nx = 4, Ny = 4, Np = 8, bs;
42 PetscMPIInt rank, size;
43
44 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
45 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
46 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
47
48 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DMDA_STENCIL_BOX, Nx, Ny, PETSC_DECIDE, PETSC_DECIDE, 1, 2, NULL, NULL, &dm));
49 PetscCall(DMSetFromOptions(dm));
50 PetscCall(DMSetUp(dm));
51 PetscCall(DMDASetUniformCoordinates(dm, -1., 1., -1., 1., -1., 1.));
52 PetscCall(DMViewFromOptions(dm, NULL, "-da_view"));
53
54 PetscCall(DMCreate(PETSC_COMM_WORLD, &sw));
55 PetscCall(PetscObjectSetName((PetscObject)sw, "Particle Grid"));
56 PetscCall(DMSetType(sw, DMSWARM));
57 PetscCall(DMSetDimension(sw, 2));
58 PetscCall(DMSwarmSetType(sw, DMSWARM_PIC));
59 PetscCall(DMSetFromOptions(sw));
60 PetscCall(DMSwarmSetCellDM(sw, dm));
61 PetscCall(DMSwarmInitializeFieldRegister(sw));
62 PetscCall(DMSwarmRegisterPetscDatatypeField(sw, "u", 1, PETSC_SCALAR));
63 PetscCall(DMSwarmFinalizeFieldRegister(sw));
64 PetscCall(DMSwarmSetLocalSizes(sw, Np, 2));
65 PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
66 dr = 1.0 / (size + 1);
67 r = (rank + 1) * dr;
68 for (PetscInt p = 0; p < Np; ++p) {
69 const PetscReal th = (p + 0.5) * 2. * PETSC_PI / Np;
70
71 coords[p * 2 + 0] = r * PetscCosReal(th);
72 coords[p * 2 + 1] = r * PetscSinReal(th);
73 }
74 PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
75 PetscCall(DMViewFromOptions(sw, NULL, "-swarm_view"));
76 PetscCall(DMSwarmPrint(sw));
77
78 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n... calling DMSwarmMigrate ...\n\n"));
79 PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
80 PetscCall(DMViewFromOptions(sw, NULL, "-swarm_view"));
81 PetscCall(DMSwarmPrint(sw));
82
83 PetscCall(DMDestroy(&sw));
84 PetscCall(DMDestroy(&dm));
85 PetscCall(PetscFinalize());
86 return 0;
87 }
88
89 /*TEST
90
91 # Swarm does not handle complex or quad
92 build:
93 requires: !complex double
94
95 test:
96 suffix: 0
97
98 TEST*/
99