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