Actual source code: ex21.c
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";
9: #include <petsc.h>
10: #include <petscdm.h>
11: #include <petscdmda.h>
12: #include <petscdmswarm.h>
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;
27: PetscOptionsGetBool(NULL,NULL,"-view",&view,NULL);
28: 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: PetscPrintf(PETSC_COMM_WORLD,"Mesh type: DMDA\n");
37: nxy = 33;
38: 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: DMDASetElementType(celldm,DMDA_ELEMENT_Q1);
42: DMSetFromOptions(celldm);
44: DMSetUp(celldm);
46: DMDASetUniformCoordinates(celldm,0.0,1.0,0.0,1.0,0.0,1.5);
48: minradius = 1.0/((PetscReal)(nxy-1));
50: 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: PetscPrintf(PETSC_COMM_WORLD,"Mesh type: DMPLEX\n");
64: DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, faces, NULL, NULL, PETSC_TRUE, &celldm);
66: /* Distribute mesh over processes */
67: DMPlexDistribute(celldm,0,NULL,&distributedMesh);
68: if (distributedMesh) {
69: DMDestroy(&celldm);
70: celldm = distributedMesh;
71: }
73: DMSetFromOptions(celldm);
75: DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,§ion);
76: DMSetLocalSection(celldm,section);
78: DMSetUp(celldm);
80: /* Calling DMPlexComputeGeometryFVM() generates the value returned by DMPlexGetMinRadius() */
81: DMPlexComputeGeometryFVM(celldm,&cellgeom,&facegeom);
82: DMPlexGetMinRadius(celldm,&minradius);
83: PetscPrintf(PETSC_COMM_WORLD,"PLEX(minradius) %1.4e\n",(double)minradius);
84: VecDestroy(&cellgeom);
85: VecDestroy(&facegeom);
86: PetscSectionDestroy(§ion);
87: }
89: /* Create the DMSwarm */
90: DMCreate(PETSC_COMM_WORLD,&swarm);
91: DMSetType(swarm,DMSWARM);
92: DMSetDimension(swarm,dim);
94: /* Configure swarm to be of type PIC */
95: DMSwarmSetType(swarm,DMSWARM_PIC);
96: DMSwarmSetCellDM(swarm,celldm);
98: /* Register two scalar fields within the DMSwarm */
99: DMSwarmRegisterPetscDatatypeField(swarm,"phi",1,PETSC_REAL);
100: DMSwarmRegisterPetscDatatypeField(swarm,"region",1,PETSC_REAL);
101: DMSwarmFinalizeFieldRegister(swarm);
103: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
104: DMSwarmSetLocalSizes(swarm,4,0);
106: /* Insert swarm coordinates cell-wise */
107: /*DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,ppcell);*/
108: 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: DMSwarmGetLocalSize(swarm,&npoints);
116: DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);
117: DMSwarmGetField(swarm,"phi",NULL,NULL,(void**)&s_phi);
118: 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: DMSwarmRestoreField(swarm,"region",NULL,NULL,(void**)&s_region);
128: DMSwarmRestoreField(swarm,"phi",NULL,NULL,(void**)&s_phi);
129: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);
130: }
132: /* Project initial value of phi onto the mesh */
133: DMSwarmProjectFields(swarm,1,fieldnames,&pfields,PETSC_FALSE);
135: if (view) {
136: /* View swarm all swarm fields using data type PETSC_REAL */
137: DMSwarmViewXDMF(swarm,"ic_dms.xmf");
139: /* View projected swarm field "phi" */
140: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
141: PetscViewerSetType(viewer,PETSCVIEWERVTK);
142: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
143: if (meshtype == 0) { /* DA */
144: PetscViewerFileSetName(viewer,"ic_dmda.vts");
145: VecView(pfields[0],viewer);
146: }
147: if (meshtype == 1) { /* PLEX */
148: PetscViewerFileSetName(viewer,"ic_dmplex.vtk");
149: DMView(celldm,viewer);
150: VecView(pfields[0],viewer);
151: }
152: PetscViewerDestroy(&viewer);
153: }
155: DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);
156: DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);
158: dt = 0.5 * minradius / PetscSqrtReal(vel[0]*vel[0] + vel[1]*vel[1]);
159: for (tk=1; tk<=nt; tk++) {
160: PetscReal *s_coor;
161: PetscInt npoints,p;
163: PetscPrintf(PETSC_COMM_WORLD,"[step %D]\n",tk);
164: /* advect with analytic prescribed (constant) velocity field */
165: DMSwarmGetLocalSize(swarm,&npoints);
166: DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);
167: for (p=0; p<npoints; p++) {
168: s_coor[2*p+0] += dt * vel[0];
169: s_coor[2*p+1] += dt * vel[1];
170: }
171: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);
173: DMSwarmMigrate(swarm,PETSC_TRUE);
175: /* Ad-hoc cell filling algorithm */
176: /*
177: The injection frequency is chosen for default DA case.
178: They will likely not be appropriate for the general case using an unstructure PLEX mesh.
179: */
180: if (tk%10 == 0) {
181: PetscReal dx = 1.0/32.0;
182: PetscInt npoints_dir_x[] = { 32, 1 };
183: PetscReal min[2],max[2];
185: min[0] = 0.5 * dx; max[0] = 0.5 * dx + 31.0 * dx;
186: min[1] = 0.5 * dx; max[1] = 0.5 * dx;
187: DMSwarmSetPointsUniformCoordinates(swarm,min,max,npoints_dir_x,ADD_VALUES);
188: }
189: if (tk%2 == 0) {
190: PetscReal dx = 1.0/32.0;
191: PetscInt npoints_dir_y[] = { 2, 31 };
192: PetscReal min[2],max[2];
194: min[0] = 0.05 * dx; max[0] = 0.5 * dx;
195: min[1] = 0.5 * dx; max[1] = 0.5 * dx + 31.0 * dx;
196: DMSwarmSetPointsUniformCoordinates(swarm,min,max,npoints_dir_y,ADD_VALUES);
197: }
199: /* Project swarm field "phi" onto the cell DM */
200: DMSwarmProjectFields(swarm,1,fieldnames,&pfields,PETSC_TRUE);
202: if (view) {
203: PetscViewer viewer;
204: char fname[PETSC_MAX_PATH_LEN];
206: /* View swarm fields */
207: PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dms.xmf",tk);
208: DMSwarmViewXDMF(swarm,fname);
210: /* View projected field */
211: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
212: PetscViewerSetType(viewer,PETSCVIEWERVTK);
213: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
215: if (meshtype == 0) { /* DA */
216: PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dmda.vts",tk);
217: PetscViewerFileSetName(viewer,fname);
218: VecView(pfields[0],viewer);
219: }
220: if (meshtype == 1) { /* PLEX */
221: PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dmplex.vtk",tk);
222: PetscViewerFileSetName(viewer,fname);
223: DMView(celldm,viewer);
224: VecView(pfields[0],viewer);
225: }
226: PetscViewerDestroy(&viewer);
227: }
229: }
230: VecDestroy(&pfields[0]);
231: PetscFree(pfields);
232: DMDestroy(&celldm);
233: DMDestroy(&swarm);
235: return 0;
236: }
238: int main(int argc,char **args)
239: {
240: PetscInt ppcell = 1;
241: PetscInt meshtype = 0;
243: PetscInitialize(&argc,&args,(char*)0,help);
244: PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL);
245: PetscOptionsGetInt(NULL,NULL,"-meshtype",&meshtype,NULL);
248: pic_advect(ppcell,meshtype);
250: PetscFinalize();
251: return 0;
252: }