1 static const char help[] = "Test initialization and migration with swarm.\n"; 2 3 #include <petscdmplex.h> 4 #include <petscdmswarm.h> 5 6 typedef struct { 7 PetscReal L[3]; /* Dimensions of the mesh bounding box */ 8 PetscBool setClosurePermutation; 9 } AppCtx; 10 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 } 20 21 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user) 22 { 23 PetscReal low[3], high[3]; 24 PetscInt cdim, d; 25 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; 33 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 } 45 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 PetscInt particleInitSize = 10; 53 PetscReal *coords, upper[3], lower[3]; 54 PetscInt *cellid, Np, dim; 55 PetscMPIInt rank, size; 56 MPI_Comm comm; 57 58 PetscFunctionBegin; 59 comm = PETSC_COMM_WORLD; 60 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 61 PetscCallMPI(MPI_Comm_size(comm, &size)); 62 PetscCall(DMGetBoundingBox(dm, lower, upper)); 63 PetscCall(DMCreate(PETSC_COMM_WORLD, sw)); 64 PetscCall(DMGetDimension(dm, &dim)); 65 PetscCall(DMSetType(*sw, DMSWARM)); 66 PetscCall(DMSetDimension(*sw, dim)); 67 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC)); 68 PetscCall(DMSwarmSetCellDM(*sw, dm)); 69 PetscCall(DMSwarmFinalizeFieldRegister(*sw)); 70 PetscCall(DMSwarmSetLocalSizes(*sw, rank == 0 ? particleInitSize : 0, 0)); 71 PetscCall(DMSetFromOptions(*sw)); 72 PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 73 PetscCall(DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid)); 74 PetscCall(DMSwarmGetLocalSize(*sw, &Np)); 75 for (PetscInt p = 0; p < Np; ++p) { 76 for (PetscInt d = 0; d < dim; ++d) { coords[p * dim + d] = 0.5 * (upper[d] - lower[d]); } 77 coords[p * dim + 1] = (upper[1] - lower[1]) / particleInitSize * p + lower[1]; 78 cellid[p] = 0; 79 } 80 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 81 PetscCall(DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid)); 82 PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view")); 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 86 /* 87 Configure the swarm on rank 0 with all particles 88 located there, then migrate where they need to be 89 */ 90 static PetscErrorCode CheckMigrate(DM sw) 91 { 92 Vec preMigrate, postMigrate, tmp; 93 PetscInt preSize, postSize; 94 PetscReal prenorm, postnorm; 95 96 PetscFunctionBeginUser; 97 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &tmp)); 98 PetscCall(VecDuplicate(tmp, &preMigrate)); 99 PetscCall(VecCopy(tmp, preMigrate)); 100 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &tmp)); 101 PetscCall(DMSwarmMigrate(sw, PETSC_TRUE)); 102 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &tmp)); 103 PetscCall(VecDuplicate(tmp, &postMigrate)); 104 PetscCall(VecCopy(tmp, postMigrate)); 105 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &tmp)); 106 PetscCall(VecGetSize(preMigrate, &preSize)); 107 PetscCall(VecGetSize(postMigrate, &postSize)); 108 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); 109 PetscCall(VecNorm(preMigrate, NORM_2, &prenorm)); 110 PetscCall(VecNorm(postMigrate, NORM_2, &postnorm)); 111 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)); 112 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Migrate check passes.\n")); 113 PetscCall(VecDestroy(&preMigrate)); 114 PetscCall(VecDestroy(&postMigrate)); 115 PetscFunctionReturn(PETSC_SUCCESS); 116 } 117 118 /* 119 Checks added points persist on migrate 120 */ 121 static PetscErrorCode CheckPointInsertion(DM sw) 122 { 123 PetscInt Np_pre, Np_post; 124 PetscMPIInt rank, size; 125 MPI_Comm comm; 126 127 PetscFunctionBeginUser; 128 comm = PETSC_COMM_WORLD; 129 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 130 PetscCallMPI(MPI_Comm_size(comm, &size)); 131 PetscCall(PetscPrintf(comm, "Basic point insertion check...\n")); 132 PetscCall(DMSwarmGetSize(sw, &Np_pre)); 133 if (rank == 0) PetscCall(DMSwarmAddPoint(sw)); 134 PetscCall(DMSwarmGetSize(sw, &Np_post)); 135 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); 136 PetscCall(CheckMigrate(sw)); 137 PetscCall(PetscPrintf(comm, "Basic point insertion check passes.\n")); 138 PetscFunctionReturn(PETSC_SUCCESS); 139 } 140 141 /* 142 Checks tie breaking works properly when a particle 143 is located at a shared boundary. The higher rank should 144 receive the particle while the lower rank deletes it. 145 146 TODO: Currently only works for 2 procs. 147 */ 148 static PetscErrorCode CheckPointInsertion_Boundary(DM sw) 149 { 150 PetscInt Np_loc_pre, Np_loc_post, dim; 151 PetscMPIInt rank, size; 152 PetscReal lbox_low[3], lbox_high[3], gbox_low[3], gbox_high[3]; 153 MPI_Comm comm; 154 DM cdm; 155 156 PetscFunctionBeginUser; 157 comm = PETSC_COMM_WORLD; 158 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 159 PetscCallMPI(MPI_Comm_size(comm, &size)); 160 PetscCall(PetscPrintf(comm, "Rank boundary point insertion check...\n")); 161 PetscCall(DMSwarmGetCellDM(sw, &cdm)); 162 PetscCall(DMGetDimension(cdm, &dim)); 163 PetscCall(DMGetBoundingBox(cdm, gbox_low, gbox_high)); 164 if (rank == 0) { 165 PetscReal *coords; 166 PetscInt adjacentdim = 0, Np; 167 168 PetscCall(DMGetLocalBoundingBox(cdm, lbox_low, lbox_high)); 169 // find a face that belongs to the neighbor. 170 for (PetscInt d = 0; d < dim; ++d) { 171 if ((gbox_high[d] - lbox_high[d]) != 0.) adjacentdim = d; 172 } 173 PetscCall(DMSwarmAddPoint(sw)); 174 PetscCall(DMSwarmGetLocalSize(sw, &Np)); 175 PetscCall(DMSwarmGetField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 176 for (PetscInt d = 0; d < dim; ++d) coords[(Np - 1) * dim + d] = 0.5 * (lbox_high[d] - lbox_low[d]); 177 coords[(Np - 1) * dim + adjacentdim] = lbox_high[adjacentdim]; 178 PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **)&coords)); 179 } 180 PetscCall(DMSwarmGetLocalSize(sw, &Np_loc_pre)); 181 PetscCall(CheckMigrate(sw)); 182 PetscCall(DMSwarmGetLocalSize(sw, &Np_loc_post)); 183 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); 184 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); 185 PetscCall(PetscPrintf(comm, "Rank boundary point insertion check passes.\n")); 186 PetscFunctionReturn(PETSC_SUCCESS); 187 } 188 189 int main(int argc, char **argv) 190 { 191 DM dm, sw; 192 MPI_Comm comm; 193 AppCtx user; 194 195 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 196 comm = PETSC_COMM_WORLD; 197 PetscCall(ProcessOptions(comm, &user)); 198 PetscCall(CreateMesh(comm, &dm, &user)); 199 PetscCall(CreateSwarm(dm, &sw, &user)); 200 PetscCall(CheckMigrate(sw)); 201 PetscCall(CheckPointInsertion(sw)); 202 PetscCall(CheckPointInsertion_Boundary(sw)); 203 PetscCall(DMDestroy(&sw)); 204 PetscCall(DMDestroy(&dm)); 205 PetscCall(PetscFinalize()); 206 return 0; 207 } 208 209 /*TEST 210 211 # Swarm does not handle complex or quad 212 build: 213 requires: !complex double 214 # swarm_migrate_hash and swarm_migrate_scan test swarm migration against point location types 215 # with a distributed mesh where ranks overlap by 1. Points in the shared boundary should 216 # be sent to the process which has the highest rank that has that portion of the domain. 217 test: 218 suffix: swarm_migrate_hash 219 nsize: 2 220 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_distribute_overlap 1 -dm_plex_box_faces 10,10,10\ 221 -dm_plex_box_lower 0.,0.,0. -dm_plex_box_upper 1.,1.,10. -dm_plex_box_bd none,none,none\ 222 -dm_plex_hash_location true 223 filter: grep -v marker | grep -v atomic | grep -v usage 224 test: 225 suffix: swarm_migrate_hash_tensor_permutation 226 nsize: 2 227 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_distribute_overlap 1 -dm_plex_box_faces 10,10,10\ 228 -dm_plex_box_lower 0.,0.,0. -dm_plex_box_upper 1.,1.,10. -dm_plex_box_bd none,none,none\ 229 -dm_plex_hash_location true -set_closure_permutation 230 filter: grep -v marker | grep -v atomic | grep -v usage 231 test: 232 suffix: swarm_migrate_scan 233 nsize: 2 234 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_distribute_overlap 1 -dm_plex_box_faces 10,10,10\ 235 -dm_plex_box_lower 0.,0.,0. -dm_plex_box_upper 1.,1.,10. -dm_plex_box_bd none,none,none\ 236 -dm_plex_hash_location false 237 filter: grep -v marker | grep -v atomic | grep -v usage 238 test: 239 suffix: swarm_migrate_scan_tensor_permutation 240 nsize: 2 241 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_distribute_overlap 1 -dm_plex_box_faces 10,10,10\ 242 -dm_plex_box_lower 0.,0.,0. -dm_plex_box_upper 1.,1.,10. -dm_plex_box_bd none,none,none\ 243 -dm_plex_hash_location false -set_closure_permutation 244 filter: grep -v marker | grep -v atomic | grep -v usage 245 TEST*/ 246