1 static char help[] = "Test point location in DA using DMSwarm\n"; 2 3 #include <petscdmda.h> 4 #include <petscdmswarm.h> 5 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 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