Actual source code: ex1.c

  1: static char help[] = "Test point location in DA using DMSwarm\n";

  3: #include <petscdmda.h>
  4: #include <petscdmswarm.h>

  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;

 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: }

 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;

 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));

 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"));

 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;

 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));

 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));

 83:   PetscCall(DMDestroy(&sw));
 84:   PetscCall(DMDestroy(&dm));
 85:   PetscCall(PetscFinalize());
 86:   return 0;
 87: }

 89: /*TEST

 91:   # Swarm does not handle complex or quad
 92:   build:
 93:     requires: !complex double

 95:   test:
 96:     suffix: 0

 98: TEST*/