particlesInDm.c

static const char help[] = "Simple Tests";

#include <petsc.h>

int main(int argc, char **argv)
{
    PetscCall(PetscInitialize(&argc, &argv, NULL, help));

    PetscInt dimensions  = 3;
    PetscInt faces[3] = {10,10, 10};
    PetscReal lower[3] = {0.0, 0.0, 0.0};
    PetscReal upper[3] = {1.0, 1.0, 10.0};
    DMBoundaryType bc[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};

    DM dm;
    PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dimensions,  PETSC_FALSE, faces, lower, upper, bc, PETSC_TRUE, &dm));
    PetscCall(DMSetFromOptions(dm));

    // get the rank
    int rank, size;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
    MPI_Comm_size(PETSC_COMM_WORLD, &size);

    // Set the max number of particles
    PetscInt particleInitSize = 10;


    // Create a dm swarm
    DM swarmDm;
    PetscCall(DMCreate(PETSC_COMM_WORLD, &swarmDm));
    PetscCall(DMSetFromOptions(swarmDm));
    PetscCall(DMSetType(swarmDm, DMSWARM));

    PetscCall(DMSetDimension(swarmDm, dimensions));

    PetscCall(DMSwarmSetType(swarmDm, DMSWARM_PIC));
    PetscCall(DMSwarmSetCellDM(swarmDm, dm));

    PetscCall(DMSwarmFinalizeFieldRegister(swarmDm));

    PetscCall(DMSwarmSetLocalSizes(swarmDm, rank == 0? particleInitSize : 0, 1));

    PetscReal* coords;
    PetscInt *cellid;
    PetscCall(DMSwarmGetField(swarmDm, DMSwarmPICField_coor, NULL, NULL, (void**)&coords));
    PetscCall(DMSwarmGetField(swarmDm, DMSwarmPICField_cellid, NULL, NULL, (void**)&cellid));

    PetscInt np;
    PetscCall(DMSwarmGetLocalSize(swarmDm, &np));
    for (PetscInt p = 0; p < np; ++p) {
        for (PetscInt d = 0; d < dimensions; ++d) {
            coords[p * dimensions + d] = 0.5*(upper[d] - lower[d]);
        }
        coords[p * dimensions + 1] = (upper[1] - lower[1])/particleInitSize* p + lower[1];
        cellid[p] = 0;
    }
    PetscCall(DMSwarmRestoreField(swarmDm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
    PetscCall(DMSwarmRestoreField(swarmDm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));


//    Vec coordVec;
//    PetscCall(DMSwarmCreateGlobalVectorFromField(swarmDm, DMSwarmPICField_coor, &coordVec));
//    VecView(coordVec, PETSC_VIEWER_STDOUT_WORLD);
//    PetscCall(DMSwarmDestroyGlobalVectorFromField(swarmDm, DMSwarmPICField_coor, &coordVec));

    DMSwarmMigrate(swarmDm, PETSC_TRUE);

    // while there are particles
    PetscInt npGlobal, npLocal;
    PetscCall(DMSwarmGetSize(swarmDm, &npGlobal));
    PetscCall(DMSwarmGetLocalSize(swarmDm, &npLocal));

    PetscPrintf(PETSC_COMM_WORLD, "Global Particles: %" PetscInt_FMT "\n", npGlobal);
    PetscPrintf(PETSC_COMM_WORLD, "########################################################################\n");
    PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Local/Global Particles %" PetscInt_FMT "/%" PetscInt_FMT " on Rank %"  PetscInt_FMT "\n",npLocal, npGlobal, rank);
    PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);

    PetscCall(DMDestroy(&swarmDm));
    PetscCall(DMDestroy(&dm));

    PetscFinalize();
    return 0;
}
/**
* Grows particles:
 * -n 4 ./particlesInDm -dm_distribute_overlap 1 -dm_plex_hash_location true
* Does not grow particles:
 * -n 4 ./particlesInDm -dm_distribute_overlap 1 -dm_plex_hash_location false

*/