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 const PetscInt dim = 2; 17 DM celldm,swarm; 18 PetscInt tk,nt = 200; 19 PetscBool view = PETSC_FALSE; 20 Vec *pfields; 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(DMSwarmProjectFields(swarm,1,fieldnames,&pfields,PETSC_FALSE)); 135 136 if (view) { 137 /* View swarm all swarm fields using data type PETSC_REAL */ 138 PetscCall(DMSwarmViewXDMF(swarm,"ic_dms.xmf")); 139 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 } 155 156 PetscCall(DMView(celldm,PETSC_VIEWER_STDOUT_WORLD)); 157 PetscCall(DMView(swarm,PETSC_VIEWER_STDOUT_WORLD)); 158 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; 163 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)); 173 174 PetscCall(DMSwarmMigrate(swarm,PETSC_TRUE)); 175 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]; 185 186 min[0] = 0.5 * dx; max[0] = 0.5 * dx + 31.0 * dx; 187 min[1] = 0.5 * dx; max[1] = 0.5 * dx; 188 PetscCall(DMSwarmSetPointsUniformCoordinates(swarm,min,max,npoints_dir_x,ADD_VALUES)); 189 } 190 if (tk%2 == 0) { 191 PetscReal dx = 1.0/32.0; 192 PetscInt npoints_dir_y[] = { 2, 31 }; 193 PetscReal min[2],max[2]; 194 195 min[0] = 0.05 * dx; max[0] = 0.5 * dx; 196 min[1] = 0.5 * dx; max[1] = 0.5 * dx + 31.0 * dx; 197 PetscCall(DMSwarmSetPointsUniformCoordinates(swarm,min,max,npoints_dir_y,ADD_VALUES)); 198 } 199 200 /* Project swarm field "phi" onto the cell DM */ 201 PetscCall(DMSwarmProjectFields(swarm,1,fieldnames,&pfields,PETSC_TRUE)); 202 203 if (view) { 204 PetscViewer viewer; 205 char fname[PETSC_MAX_PATH_LEN]; 206 207 /* View swarm fields */ 208 PetscCall(PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4" PetscInt_FMT "_dms.xmf",tk)); 209 PetscCall(DMSwarmViewXDMF(swarm,fname)); 210 211 /* View projected field */ 212 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD,&viewer)); 213 PetscCall(PetscViewerSetType(viewer,PETSCVIEWERVTK)); 214 PetscCall(PetscViewerFileSetMode(viewer,FILE_MODE_WRITE)); 215 216 if (meshtype == 0) { /* DA */ 217 PetscCall(PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4" PetscInt_FMT "_dmda.vts",tk)); 218 PetscCall(PetscViewerFileSetName(viewer,fname)); 219 PetscCall(VecView(pfields[0],viewer)); 220 } 221 if (meshtype == 1) { /* PLEX */ 222 PetscCall(PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4" PetscInt_FMT "_dmplex.vtk",tk)); 223 PetscCall(PetscViewerFileSetName(viewer,fname)); 224 PetscCall(DMView(celldm,viewer)); 225 PetscCall(VecView(pfields[0],viewer)); 226 } 227 PetscCall(PetscViewerDestroy(&viewer)); 228 } 229 230 } 231 PetscCall(VecDestroy(&pfields[0])); 232 PetscCall(PetscFree(pfields)); 233 PetscCall(DMDestroy(&celldm)); 234 PetscCall(DMDestroy(&swarm)); 235 236 PetscFunctionReturn(0); 237 } 238 239 int main(int argc,char **args) 240 { 241 PetscInt ppcell = 1; 242 PetscInt meshtype = 0; 243 244 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 245 PetscCall(PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL)); 246 PetscCall(PetscOptionsGetInt(NULL,NULL,"-meshtype",&meshtype,NULL)); 247 PetscCheck(meshtype <= 1,PETSC_COMM_WORLD,PETSC_ERR_USER,"-meshtype <value> must be 0 or 1"); 248 249 PetscCall(pic_advect(ppcell,meshtype)); 250 251 PetscCall(PetscFinalize()); 252 return 0; 253 } 254