Actual source code: ex21.c

      1: static char help[] = "DMSwarm-PIC demonstrator of advecting points within cell DM defined by a DA or PLEX object \n\
      2: Options: \n\
      3: -ppcell   : Number of times to sub-divide the reference cell when layout the initial particle coordinates \n\
      4: -meshtype : 0 ==> DA , 1 ==> PLEX \n\
      5: -nt       : Number of timestep to perform \n\
      6: -view     : Write out initial condition and time dependent data \n";

      8: #include <petsc.h>
      9: #include <petscdm.h>
     10: #include <petscdmda.h>
     11: #include <petscdmswarm.h>

     13: PetscErrorCode pic_advect(PetscInt ppcell, PetscInt meshtype)
     14: {
     15:   const PetscInt dim = 2;
     16:   DM             celldm, swarm;
     17:   PetscInt       tk, nt = 200;
     18:   PetscBool      view = PETSC_FALSE;
     19:   Vec            pfields[1];
     20:   PetscReal      minradius;
     21:   PetscReal      dt;
     22:   PetscReal      vel[]        = {1.0, 0.16};
     23:   const char    *fieldnames[] = {"phi"};
     24:   PetscViewer    viewer;

     26:   PetscFunctionBegin;
     27:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL));
     28:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL));

     30:   /* Create the background cell DM */
     31:   if (meshtype == 0) { /* DA */
     32:     PetscInt nxy;
     33:     PetscInt dof           = 1;
     34:     PetscInt stencil_width = 1;

     36:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mesh type: DMDA\n"));
     37:     nxy = 33;
     38:     PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, nxy, nxy, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &celldm));

     40:     PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1));

     42:     PetscCall(DMSetFromOptions(celldm));

     44:     PetscCall(DMSetUp(celldm));

     46:     PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.5));

     48:     minradius = 1.0 / ((PetscReal)(nxy - 1));

     50:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "DA(minradius) %1.4e\n", (double)minradius));
     51:   }

     53:   if (meshtype == 1) { /* PLEX */
     54:     DM           distributedMesh = NULL;
     55:     PetscInt     numComp[]       = {1};
     56:     PetscInt     numDof[]        = {1, 0, 0}; /* vert, edge, cell */
     57:     PetscInt     faces[]         = {1, 1, 1};
     58:     PetscInt     numBC           = 0;
     59:     PetscSection section;
     60:     Vec          cellgeom = NULL;
     61:     Vec          facegeom = NULL;

     63:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mesh type: DMPLEX\n"));
     64:     PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, faces, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &celldm));

     66:     /* Distribute mesh over processes */
     67:     PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh));
     68:     if (distributedMesh) {
     69:       PetscCall(DMDestroy(&celldm));
     70:       celldm = distributedMesh;
     71:     }

     73:     PetscCall(DMSetFromOptions(celldm));

     75:     PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &section));
     76:     PetscCall(DMSetLocalSection(celldm, section));

     78:     PetscCall(DMSetUp(celldm));

     80:     /* Calling DMPlexComputeGeometryFVM() generates the value returned by DMPlexGetMinRadius() */
     81:     PetscCall(DMPlexComputeGeometryFVM(celldm, &cellgeom, &facegeom));
     82:     PetscCall(DMPlexGetMinRadius(celldm, &minradius));
     83:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "PLEX(minradius) %1.4e\n", (double)minradius));
     84:     PetscCall(VecDestroy(&cellgeom));
     85:     PetscCall(VecDestroy(&facegeom));
     86:     PetscCall(PetscSectionDestroy(&section));
     87:   }

     89:   /* Create the DMSwarm */
     90:   PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
     91:   PetscCall(DMSetType(swarm, DMSWARM));
     92:   PetscCall(DMSetDimension(swarm, dim));

     94:   /* Configure swarm to be of type PIC */
     95:   PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
     96:   PetscCall(DMSwarmSetCellDM(swarm, celldm));

     98:   /* Register two scalar fields within the DMSwarm */
     99:   PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "phi", 1, PETSC_REAL));
    100:   PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "region", 1, PETSC_REAL));
    101:   PetscCall(DMSwarmFinalizeFieldRegister(swarm));

    103:   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
    104:   PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));

    106:   /* Insert swarm coordinates cell-wise */
    107:   /*PetscCall(DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,ppcell));*/
    108:   PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell));

    110:   /* Define initial conditions for th swarm fields "phi" and "region" */
    111:   {
    112:     PetscReal *s_coor, *s_phi, *s_region;
    113:     PetscInt   npoints, p;

    115:     PetscCall(DMSwarmGetLocalSize(swarm, &npoints));
    116:     PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor));
    117:     PetscCall(DMSwarmGetField(swarm, "phi", NULL, NULL, (void **)&s_phi));
    118:     PetscCall(DMSwarmGetField(swarm, "region", NULL, NULL, (void **)&s_region));
    119:     for (p = 0; p < npoints; p++) {
    120:       PetscReal pos[2];
    121:       pos[0] = s_coor[2 * p + 0];
    122:       pos[1] = s_coor[2 * p + 1];

    124:       s_region[p] = 1.0;
    125:       s_phi[p]    = 1.0 + PetscExpReal(-200.0 * ((pos[0] - 0.5) * (pos[0] - 0.5) + (pos[1] - 0.5) * (pos[1] - 0.5)));
    126:     }
    127:     PetscCall(DMSwarmRestoreField(swarm, "region", NULL, NULL, (void **)&s_region));
    128:     PetscCall(DMSwarmRestoreField(swarm, "phi", NULL, NULL, (void **)&s_phi));
    129:     PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor));
    130:   }

    132:   /* Project initial value of phi onto the mesh */
    133:   PetscCall(DMCreateGlobalVector(celldm, &pfields[0]));
    134:   PetscCall(DMSwarmProjectFields(swarm, NULL, 1, fieldnames, pfields, SCATTER_FORWARD));

    136:   if (view) {
    137:     /* View swarm all swarm fields using data type PETSC_REAL */
    138:     PetscCall(DMSwarmViewXDMF(swarm, "ic_dms.xmf"));

    140:     /* View projected swarm field "phi" */
    141:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
    142:     PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
    143:     PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
    144:     if (meshtype == 0) { /* DA */
    145:       PetscCall(PetscViewerFileSetName(viewer, "ic_dmda.vts"));
    146:       PetscCall(VecView(pfields[0], viewer));
    147:     }
    148:     if (meshtype == 1) { /* PLEX */
    149:       PetscCall(PetscViewerFileSetName(viewer, "ic_dmplex.vtk"));
    150:       PetscCall(DMView(celldm, viewer));
    151:       PetscCall(VecView(pfields[0], viewer));
    152:     }
    153:     PetscCall(PetscViewerDestroy(&viewer));
    154:   }

    156:   PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD));
    157:   PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD));

    159:   dt = 0.5 * minradius / PetscSqrtReal(vel[0] * vel[0] + vel[1] * vel[1]);
    160:   for (tk = 1; tk <= nt; tk++) {
    161:     PetscReal *s_coor;
    162:     PetscInt   npoints, p;

    164:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[step %" PetscInt_FMT "]\n", tk));
    165:     /* advect with analytic prescribed (constant) velocity field */
    166:     PetscCall(DMSwarmGetLocalSize(swarm, &npoints));
    167:     PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor));
    168:     for (p = 0; p < npoints; p++) {
    169:       s_coor[2 * p + 0] += dt * vel[0];
    170:       s_coor[2 * p + 1] += dt * vel[1];
    171:     }
    172:     PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor));

    174:     PetscCall(DMSwarmMigrate(swarm, PETSC_TRUE));

    176:     /* Ad-hoc cell filling algorithm */
    177:     /*
    178:        The injection frequency is chosen for default DA case.
    179:        They will likely not be appropriate for the general case using an unstructure PLEX mesh.
    180:     */
    181:     if (tk % 10 == 0) {
    182:       PetscReal dx              = 1.0 / 32.0;
    183:       PetscInt  npoints_dir_x[] = {32, 1};
    184:       PetscReal min[2], max[2];

    186:       min[0] = 0.5 * dx;
    187:       max[0] = 0.5 * dx + 31.0 * dx;
    188:       min[1] = 0.5 * dx;
    189:       max[1] = 0.5 * dx;
    190:       PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, npoints_dir_x, ADD_VALUES));
    191:     }
    192:     if (tk % 2 == 0) {
    193:       PetscReal dx              = 1.0 / 32.0;
    194:       PetscInt  npoints_dir_y[] = {2, 31};
    195:       PetscReal min[2], max[2];

    197:       min[0] = 0.05 * dx;
    198:       max[0] = 0.5 * dx;
    199:       min[1] = 0.5 * dx;
    200:       max[1] = 0.5 * dx + 31.0 * dx;
    201:       PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, npoints_dir_y, ADD_VALUES));
    202:     }

    204:     /* Project swarm field "phi" onto the cell DM */
    205:     PetscCall(DMSwarmProjectFields(swarm, NULL, 1, fieldnames, pfields, SCATTER_FORWARD));

    207:     if (view) {
    208:       PetscViewer viewer;
    209:       char        fname[PETSC_MAX_PATH_LEN];

    211:       /* View swarm fields */
    212:       PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dms.xmf", tk));
    213:       PetscCall(DMSwarmViewXDMF(swarm, fname));

    215:       /* View projected field */
    216:       PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
    217:       PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
    218:       PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));

    220:       if (meshtype == 0) { /* DA */
    221:         PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dmda.vts", tk));
    222:         PetscCall(PetscViewerFileSetName(viewer, fname));
    223:         PetscCall(VecView(pfields[0], viewer));
    224:       }
    225:       if (meshtype == 1) { /* PLEX */
    226:         PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dmplex.vtk", tk));
    227:         PetscCall(PetscViewerFileSetName(viewer, fname));
    228:         PetscCall(DMView(celldm, viewer));
    229:         PetscCall(VecView(pfields[0], viewer));
    230:       }
    231:       PetscCall(PetscViewerDestroy(&viewer));
    232:     }
    233:   }
    234:   PetscCall(VecDestroy(&pfields[0]));
    235:   PetscCall(DMDestroy(&celldm));
    236:   PetscCall(DMDestroy(&swarm));
    237:   PetscFunctionReturn(PETSC_SUCCESS);
    238: }

    240: int main(int argc, char **args)
    241: {
    242:   PetscInt ppcell   = 1;
    243:   PetscInt meshtype = 0;

    245:   PetscFunctionBeginUser;
    246:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
    247:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL));
    248:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-meshtype", &meshtype, NULL));
    249:   PetscCheck(meshtype <= 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "-meshtype <value> must be 0 or 1");

    251:   PetscCall(pic_advect(ppcell, meshtype));

    253:   PetscCall(PetscFinalize());
    254:   return 0;
    255: }