Actual source code: ex1.c
1: static const char help[] = "Test initialization and migration with swarm.\n";
3: #include <petscdmplex.h>
4: #include <petscdmswarm.h>
6: typedef struct {
7: PetscReal L[3]; /* Dimensions of the mesh bounding box */
8: PetscBool setClosurePermutation;
9: } AppCtx;
11: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
12: {
13: PetscFunctionBeginUser;
14: PetscOptionsBegin(comm, "", "Swarm configuration options", "DMSWARM");
15: options->setClosurePermutation = PETSC_FALSE;
16: PetscCall(PetscOptionsBool("-set_closure_permutation", "Set the closure permutation to tensor ordering", NULL, options->setClosurePermutation, &options->setClosurePermutation, NULL));
17: PetscOptionsEnd();
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
21: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
22: {
23: PetscReal low[3], high[3];
24: PetscInt cdim, d;
26: PetscFunctionBeginUser;
27: PetscCall(DMCreate(comm, dm));
28: PetscCall(DMSetType(*dm, DMPLEX));
29: PetscCall(DMSetFromOptions(*dm));
30: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
31: if (user->setClosurePermutation) {
32: DM cdm;
34: // -- Set tensor permutation
35: PetscCall(DMGetCoordinateDM(*dm, &cdm));
36: PetscCall(DMPlexSetClosurePermutationTensor(*dm, PETSC_DETERMINE, NULL));
37: PetscCall(DMPlexSetClosurePermutationTensor(cdm, PETSC_DETERMINE, NULL));
38: }
39: PetscCall(DMGetCoordinateDim(*dm, &cdim));
40: PetscCall(DMGetBoundingBox(*dm, low, high));
41: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "dim: %" PetscInt_FMT "\n", cdim));
42: for (d = 0; d < cdim; ++d) user->L[d] = high[d] - low[d];
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*
47: This function initializes all particles on rank 0.
48: They are sent to other ranks to test migration across non nearest neighbors
49: */
50: static PetscErrorCode CreateSwarm(DM dm, DM *sw, AppCtx *user)
51: {
52: DMSwarmCellDM celldm;
53: PetscInt particleInitSize = 10;
54: PetscReal *coords, upper[3], lower[3];
55: PetscInt *swarm_cellid, Np, dim, Nfc;
56: PetscMPIInt rank, size;
57: MPI_Comm comm;
58: const char **coordFields, *cellid;
60: PetscFunctionBegin;
61: comm = PETSC_COMM_WORLD;
62: PetscCallMPI(MPI_Comm_rank(comm, &rank));
63: PetscCallMPI(MPI_Comm_size(comm, &size));
64: PetscCall(DMGetBoundingBox(dm, lower, upper));
65: PetscCall(DMCreate(PETSC_COMM_WORLD, sw));
66: PetscCall(DMGetDimension(dm, &dim));
67: PetscCall(DMSetType(*sw, DMSWARM));
68: PetscCall(DMSetDimension(*sw, dim));
69: PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
70: PetscCall(DMSwarmSetCellDM(*sw, dm));
71: PetscCall(DMSwarmFinalizeFieldRegister(*sw));
72: PetscCall(DMSwarmSetLocalSizes(*sw, rank == 0 ? particleInitSize : 0, 0));
73: PetscCall(DMSetFromOptions(*sw));
74: PetscCall(DMSwarmGetCellDMActive(*sw, &celldm));
75: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
76: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
77: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
78: PetscCall(DMSwarmGetField(*sw, coordFields[0], NULL, NULL, (void **)&coords));
79: PetscCall(DMSwarmGetField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
80: PetscCall(DMSwarmGetLocalSize(*sw, &Np));
81: for (PetscInt p = 0; p < Np; ++p) {
82: for (PetscInt d = 0; d < dim; ++d) { coords[p * dim + d] = 0.5 * (upper[d] - lower[d]); }
83: coords[p * dim + 1] = (upper[1] - lower[1]) / particleInitSize * p + lower[1];
84: swarm_cellid[p] = 0;
85: }
86: PetscCall(DMSwarmRestoreField(*sw, coordFields[0], NULL, NULL, (void **)&coords));
87: PetscCall(DMSwarmRestoreField(*sw, cellid, NULL, NULL, (void **)&swarm_cellid));
88: PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: /*
93: Configure the swarm on rank 0 with all particles
94: located there, then migrate where they need to be
95: */
96: static PetscErrorCode CheckMigrate(DM sw)
97: {
98: Vec preMigrate, postMigrate, tmp;
99: PetscInt preSize, postSize;
100: PetscReal prenorm, postnorm;
102: PetscFunctionBeginUser;
103: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &tmp));
104: PetscCall(VecDuplicate(tmp, &preMigrate));
105: PetscCall(VecCopy(tmp, preMigrate));
106: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &tmp));
107: PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
108: PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &tmp));
109: PetscCall(VecDuplicate(tmp, &postMigrate));
110: PetscCall(VecCopy(tmp, postMigrate));
111: PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &tmp));
112: PetscCall(VecGetSize(preMigrate, &preSize));
113: PetscCall(VecGetSize(postMigrate, &postSize));
114: PetscCheck(preSize == postSize, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Particles either lost or duplicated. Pre migrate global size %" PetscInt_FMT " != Post migrate size %" PetscInt_FMT, preSize, postSize);
115: PetscCall(VecNorm(preMigrate, NORM_2, &prenorm));
116: PetscCall(VecNorm(postMigrate, NORM_2, &postnorm));
117: PetscCheck(PetscAbsReal(prenorm - postnorm) < 100. * PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_COR, "Particle coordinates corrupted in migrate with abs(norm(pre) - norm(post)) = %.16g", PetscAbsReal(prenorm - postnorm));
118: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Migrate check passes.\n"));
119: PetscCall(VecDestroy(&preMigrate));
120: PetscCall(VecDestroy(&postMigrate));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: /*
125: Checks added points persist on migrate
126: */
127: static PetscErrorCode CheckPointInsertion(DM sw)
128: {
129: PetscInt Np_pre, Np_post;
130: PetscMPIInt rank, size;
131: MPI_Comm comm;
133: PetscFunctionBeginUser;
134: comm = PETSC_COMM_WORLD;
135: PetscCallMPI(MPI_Comm_rank(comm, &rank));
136: PetscCallMPI(MPI_Comm_size(comm, &size));
137: PetscCall(PetscPrintf(comm, "Basic point insertion check...\n"));
138: PetscCall(DMSwarmGetSize(sw, &Np_pre));
139: if (rank == 0) PetscCall(DMSwarmAddPoint(sw));
140: PetscCall(DMSwarmGetSize(sw, &Np_post));
141: PetscCheck(Np_post == (Np_pre + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Particle insertion failed. Global size pre insertion: %" PetscInt_FMT " global size post insertion: %" PetscInt_FMT, Np_pre, Np_post);
142: PetscCall(CheckMigrate(sw));
143: PetscCall(PetscPrintf(comm, "Basic point insertion check passes.\n"));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*
148: Checks tie breaking works properly when a particle
149: is located at a shared boundary. The higher rank should
150: receive the particle while the lower rank deletes it.
152: TODO: Currently only works for 2 procs.
153: */
154: static PetscErrorCode CheckPointInsertion_Boundary(DM sw)
155: {
156: PetscInt Np_loc_pre, Np_loc_post, dim;
157: PetscMPIInt rank, size;
158: PetscReal lbox_low[3], lbox_high[3], gbox_low[3], gbox_high[3];
159: MPI_Comm comm;
160: DM cdm;
162: PetscFunctionBeginUser;
163: comm = PETSC_COMM_WORLD;
164: PetscCallMPI(MPI_Comm_rank(comm, &rank));
165: PetscCallMPI(MPI_Comm_size(comm, &size));
166: PetscCall(PetscPrintf(comm, "Rank boundary point insertion check...\n"));
167: PetscCall(DMSwarmGetCellDM(sw, &cdm));
168: PetscCall(DMGetDimension(cdm, &dim));
169: PetscCall(DMGetBoundingBox(cdm, gbox_low, gbox_high));
170: if (rank == 0) {
171: PetscReal *coords;
172: PetscInt adjacentdim = 0, Np;
174: PetscCall(DMGetLocalBoundingBox(cdm, lbox_low, lbox_high));
175: // find a face that belongs to the neighbor.
176: for (PetscInt d = 0; d < dim; ++d) {
177: if ((gbox_high[d] - lbox_high[d]) != 0.) adjacentdim = d;
178: }
179: PetscCall(DMSwarmAddPoint(sw));
180: PetscCall(DMSwarmGetLocalSize(sw, &Np));
181: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
182: for (PetscInt d = 0; d < dim; ++d) coords[(Np - 1) * dim + d] = 0.5 * (lbox_high[d] - lbox_low[d]);
183: coords[(Np - 1) * dim + adjacentdim] = lbox_high[adjacentdim];
184: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
185: }
186: PetscCall(DMSwarmGetLocalSize(sw, &Np_loc_pre));
187: PetscCall(CheckMigrate(sw));
188: PetscCall(DMSwarmGetLocalSize(sw, &Np_loc_post));
189: if (rank == 0) PetscCheck(Np_loc_pre == (Np_loc_post + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Migration tie breaking failed on rank %d. Particle on boundary not sent.", rank);
190: if (rank == 1) PetscCheck(Np_loc_pre == (Np_loc_post - 1), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Migration tie breaking failed on rank %d. Particle on boundary not received.", rank);
191: PetscCall(PetscPrintf(comm, "Rank boundary point insertion check passes.\n"));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: int main(int argc, char **argv)
196: {
197: DM dm, sw;
198: MPI_Comm comm;
199: AppCtx user;
201: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
202: comm = PETSC_COMM_WORLD;
203: PetscCall(ProcessOptions(comm, &user));
204: PetscCall(CreateMesh(comm, &dm, &user));
205: PetscCall(CreateSwarm(dm, &sw, &user));
206: PetscCall(CheckMigrate(sw));
207: PetscCall(CheckPointInsertion(sw));
208: PetscCall(CheckPointInsertion_Boundary(sw));
209: PetscCall(DMDestroy(&sw));
210: PetscCall(DMDestroy(&dm));
211: PetscCall(PetscFinalize());
212: return 0;
213: }
215: /*TEST
217: # Swarm does not handle complex or quad
218: build:
219: requires: !complex double
220: # swarm_migrate_hash and swarm_migrate_scan test swarm migration against point location types
221: # with a distributed mesh where ranks overlap by 1. Points in the shared boundary should
222: # be sent to the process which has the highest rank that has that portion of the domain.
223: test:
224: suffix: swarm_migrate_hash
225: nsize: 2
226: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_distribute_overlap 1 -dm_plex_box_faces 10,10,10\
227: -dm_plex_box_lower 0.,0.,0. -dm_plex_box_upper 1.,1.,10. -dm_plex_box_bd none,none,none\
228: -dm_plex_hash_location true
229: filter: grep -v marker | grep -v atomic | grep -v usage
230: test:
231: suffix: swarm_migrate_hash_tensor_permutation
232: nsize: 2
233: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_distribute_overlap 1 -dm_plex_box_faces 10,10,10\
234: -dm_plex_box_lower 0.,0.,0. -dm_plex_box_upper 1.,1.,10. -dm_plex_box_bd none,none,none\
235: -dm_plex_hash_location true -set_closure_permutation
236: filter: grep -v marker | grep -v atomic | grep -v usage
237: test:
238: suffix: swarm_migrate_scan
239: nsize: 2
240: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_distribute_overlap 1 -dm_plex_box_faces 10,10,10\
241: -dm_plex_box_lower 0.,0.,0. -dm_plex_box_upper 1.,1.,10. -dm_plex_box_bd none,none,none\
242: -dm_plex_hash_location false
243: filter: grep -v marker | grep -v atomic | grep -v usage
244: test:
245: suffix: swarm_migrate_scan_tensor_permutation
246: nsize: 2
247: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_distribute_overlap 1 -dm_plex_box_faces 10,10,10\
248: -dm_plex_box_lower 0.,0.,0. -dm_plex_box_upper 1.,1.,10. -dm_plex_box_bd none,none,none\
249: -dm_plex_hash_location false -set_closure_permutation
250: filter: grep -v marker | grep -v atomic | grep -v usage
251: TEST*/