1 2 static char help[] = "DMSwarm-PIC demonstrator of advecting points within cell DM defined by a DA or PLEX object \n\ 3 Options: \n\ 4 -ppcell : Number of times to sub-divide the reference cell when layout the initial particle coordinates \n\ 5 -meshtype : 0 ==> DA , 1 ==> PLEX \n\ 6 -nt : Number of timestep to perform \n\ 7 -view : Write out initial condition and time dependent data \n"; 8 9 #include <petsc.h> 10 #include <petscdm.h> 11 #include <petscdmda.h> 12 #include <petscdmswarm.h> 13 14 PetscErrorCode pic_advect(PetscInt ppcell, PetscInt meshtype) 15 { 16 const PetscInt dim = 2; 17 DM celldm, swarm; 18 PetscInt tk, nt = 200; 19 PetscBool view = PETSC_FALSE; 20 Vec pfields[1]; 21 PetscReal minradius; 22 PetscReal dt; 23 PetscReal vel[] = {1.0, 0.16}; 24 const char *fieldnames[] = {"phi"}; 25 PetscViewer viewer; 26 27 PetscFunctionBegin; 28 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL)); 29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL)); 30 31 /* Create the background cell DM */ 32 if (meshtype == 0) { /* DA */ 33 PetscInt nxy; 34 PetscInt dof = 1; 35 PetscInt stencil_width = 1; 36 37 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mesh type: DMDA\n")); 38 nxy = 33; 39 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 41 PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1)); 42 43 PetscCall(DMSetFromOptions(celldm)); 44 45 PetscCall(DMSetUp(celldm)); 46 47 PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.5)); 48 49 minradius = 1.0 / ((PetscReal)(nxy - 1)); 50 51 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "DA(minradius) %1.4e\n", (double)minradius)); 52 } 53 54 if (meshtype == 1) { /* PLEX */ 55 DM distributedMesh = NULL; 56 PetscInt numComp[] = {1}; 57 PetscInt numDof[] = {1, 0, 0}; /* vert, edge, cell */ 58 PetscInt faces[] = {1, 1, 1}; 59 PetscInt numBC = 0; 60 PetscSection section; 61 Vec cellgeom = NULL; 62 Vec facegeom = NULL; 63 64 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mesh type: DMPLEX\n")); 65 PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, faces, NULL, NULL, PETSC_TRUE, &celldm)); 66 67 /* Distribute mesh over processes */ 68 PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh)); 69 if (distributedMesh) { 70 PetscCall(DMDestroy(&celldm)); 71 celldm = distributedMesh; 72 } 73 74 PetscCall(DMSetFromOptions(celldm)); 75 76 PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion)); 77 PetscCall(DMSetLocalSection(celldm, section)); 78 79 PetscCall(DMSetUp(celldm)); 80 81 /* Calling DMPlexComputeGeometryFVM() generates the value returned by DMPlexGetMinRadius() */ 82 PetscCall(DMPlexComputeGeometryFVM(celldm, &cellgeom, &facegeom)); 83 PetscCall(DMPlexGetMinRadius(celldm, &minradius)); 84 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "PLEX(minradius) %1.4e\n", (double)minradius)); 85 PetscCall(VecDestroy(&cellgeom)); 86 PetscCall(VecDestroy(&facegeom)); 87 PetscCall(PetscSectionDestroy(§ion)); 88 } 89 90 /* Create the DMSwarm */ 91 PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm)); 92 PetscCall(DMSetType(swarm, DMSWARM)); 93 PetscCall(DMSetDimension(swarm, dim)); 94 95 /* Configure swarm to be of type PIC */ 96 PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC)); 97 PetscCall(DMSwarmSetCellDM(swarm, celldm)); 98 99 /* Register two scalar fields within the DMSwarm */ 100 PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "phi", 1, PETSC_REAL)); 101 PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "region", 1, PETSC_REAL)); 102 PetscCall(DMSwarmFinalizeFieldRegister(swarm)); 103 104 /* Set initial local sizes of the DMSwarm with a buffer length of zero */ 105 PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0)); 106 107 /* Insert swarm coordinates cell-wise */ 108 /*PetscCall(DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,ppcell));*/ 109 PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell)); 110 111 /* Define initial conditions for th swarm fields "phi" and "region" */ 112 { 113 PetscReal *s_coor, *s_phi, *s_region; 114 PetscInt npoints, p; 115 116 PetscCall(DMSwarmGetLocalSize(swarm, &npoints)); 117 PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor)); 118 PetscCall(DMSwarmGetField(swarm, "phi", NULL, NULL, (void **)&s_phi)); 119 PetscCall(DMSwarmGetField(swarm, "region", NULL, NULL, (void **)&s_region)); 120 for (p = 0; p < npoints; p++) { 121 PetscReal pos[2]; 122 pos[0] = s_coor[2 * p + 0]; 123 pos[1] = s_coor[2 * p + 1]; 124 125 s_region[p] = 1.0; 126 s_phi[p] = 1.0 + PetscExpReal(-200.0 * ((pos[0] - 0.5) * (pos[0] - 0.5) + (pos[1] - 0.5) * (pos[1] - 0.5))); 127 } 128 PetscCall(DMSwarmRestoreField(swarm, "region", NULL, NULL, (void **)&s_region)); 129 PetscCall(DMSwarmRestoreField(swarm, "phi", NULL, NULL, (void **)&s_phi)); 130 PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor)); 131 } 132 133 /* Project initial value of phi onto the mesh */ 134 PetscCall(DMCreateGlobalVector(celldm, &pfields[0])); 135 PetscCall(DMSwarmProjectFields(swarm, 1, fieldnames, pfields, SCATTER_FORWARD)); 136 137 if (view) { 138 /* View swarm all swarm fields using data type PETSC_REAL */ 139 PetscCall(DMSwarmViewXDMF(swarm, "ic_dms.xmf")); 140 141 /* View projected swarm field "phi" */ 142 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer)); 143 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK)); 144 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE)); 145 if (meshtype == 0) { /* DA */ 146 PetscCall(PetscViewerFileSetName(viewer, "ic_dmda.vts")); 147 PetscCall(VecView(pfields[0], viewer)); 148 } 149 if (meshtype == 1) { /* PLEX */ 150 PetscCall(PetscViewerFileSetName(viewer, "ic_dmplex.vtk")); 151 PetscCall(DMView(celldm, viewer)); 152 PetscCall(VecView(pfields[0], viewer)); 153 } 154 PetscCall(PetscViewerDestroy(&viewer)); 155 } 156 157 PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD)); 158 PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD)); 159 160 dt = 0.5 * minradius / PetscSqrtReal(vel[0] * vel[0] + vel[1] * vel[1]); 161 for (tk = 1; tk <= nt; tk++) { 162 PetscReal *s_coor; 163 PetscInt npoints, p; 164 165 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[step %" PetscInt_FMT "]\n", tk)); 166 /* advect with analytic prescribed (constant) velocity field */ 167 PetscCall(DMSwarmGetLocalSize(swarm, &npoints)); 168 PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor)); 169 for (p = 0; p < npoints; p++) { 170 s_coor[2 * p + 0] += dt * vel[0]; 171 s_coor[2 * p + 1] += dt * vel[1]; 172 } 173 PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor)); 174 175 PetscCall(DMSwarmMigrate(swarm, PETSC_TRUE)); 176 177 /* Ad-hoc cell filling algorithm */ 178 /* 179 The injection frequency is chosen for default DA case. 180 They will likely not be appropriate for the general case using an unstructure PLEX mesh. 181 */ 182 if (tk % 10 == 0) { 183 PetscReal dx = 1.0 / 32.0; 184 PetscInt npoints_dir_x[] = {32, 1}; 185 PetscReal min[2], max[2]; 186 187 min[0] = 0.5 * dx; 188 max[0] = 0.5 * dx + 31.0 * dx; 189 min[1] = 0.5 * dx; 190 max[1] = 0.5 * dx; 191 PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, npoints_dir_x, ADD_VALUES)); 192 } 193 if (tk % 2 == 0) { 194 PetscReal dx = 1.0 / 32.0; 195 PetscInt npoints_dir_y[] = {2, 31}; 196 PetscReal min[2], max[2]; 197 198 min[0] = 0.05 * dx; 199 max[0] = 0.5 * dx; 200 min[1] = 0.5 * dx; 201 max[1] = 0.5 * dx + 31.0 * dx; 202 PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, npoints_dir_y, ADD_VALUES)); 203 } 204 205 /* Project swarm field "phi" onto the cell DM */ 206 PetscCall(DMSwarmProjectFields(swarm, 1, fieldnames, pfields, SCATTER_FORWARD)); 207 208 if (view) { 209 PetscViewer viewer; 210 char fname[PETSC_MAX_PATH_LEN]; 211 212 /* View swarm fields */ 213 PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dms.xmf", tk)); 214 PetscCall(DMSwarmViewXDMF(swarm, fname)); 215 216 /* View projected field */ 217 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer)); 218 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK)); 219 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE)); 220 221 if (meshtype == 0) { /* DA */ 222 PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dmda.vts", tk)); 223 PetscCall(PetscViewerFileSetName(viewer, fname)); 224 PetscCall(VecView(pfields[0], viewer)); 225 } 226 if (meshtype == 1) { /* PLEX */ 227 PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dmplex.vtk", tk)); 228 PetscCall(PetscViewerFileSetName(viewer, fname)); 229 PetscCall(DMView(celldm, viewer)); 230 PetscCall(VecView(pfields[0], viewer)); 231 } 232 PetscCall(PetscViewerDestroy(&viewer)); 233 } 234 } 235 PetscCall(VecDestroy(&pfields[0])); 236 PetscCall(DMDestroy(&celldm)); 237 PetscCall(DMDestroy(&swarm)); 238 239 PetscFunctionReturn(PETSC_SUCCESS); 240 } 241 242 int main(int argc, char **args) 243 { 244 PetscInt ppcell = 1; 245 PetscInt meshtype = 0; 246 247 PetscFunctionBeginUser; 248 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 249 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL)); 250 PetscCall(PetscOptionsGetInt(NULL, NULL, "-meshtype", &meshtype, NULL)); 251 PetscCheck(meshtype <= 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "-meshtype <value> must be 0 or 1"); 252 253 PetscCall(pic_advect(ppcell, meshtype)); 254 255 PetscCall(PetscFinalize()); 256 return 0; 257 } 258