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: {
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;
29: PetscOptionsGetBool(NULL,NULL,"-view",&view,NULL);
30: PetscOptionsGetInt(NULL,NULL,"-nt",&nt,NULL);
32: /* Create the background cell DM */
33: if (meshtype == 0) { /* DA */
34: PetscInt nxy;
35: PetscInt dof = 1;
36: PetscInt stencil_width = 1;
38: PetscPrintf(PETSC_COMM_WORLD,"Mesh type: DMDA\n");
39: nxy = 33;
40: 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);
42: DMDASetElementType(celldm,DMDA_ELEMENT_Q1);
44: DMSetFromOptions(celldm);
46: DMSetUp(celldm);
48: DMDASetUniformCoordinates(celldm,0.0,1.0,0.0,1.0,0.0,1.5);
50: minradius = 1.0/((PetscReal)(nxy-1));
52: PetscPrintf(PETSC_COMM_WORLD,"DA(minradius) %1.4e\n",(double)minradius);
53: }
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;
65: PetscPrintf(PETSC_COMM_WORLD,"Mesh type: DMPLEX\n");
66: DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, faces, NULL, NULL, PETSC_TRUE, &celldm);
68: /* Distribute mesh over processes */
69: DMPlexDistribute(celldm,0,NULL,&distributedMesh);
70: if (distributedMesh) {
71: DMDestroy(&celldm);
72: celldm = distributedMesh;
73: }
75: DMSetFromOptions(celldm);
77: DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,§ion);
78: DMSetLocalSection(celldm,section);
80: DMSetUp(celldm);
82: /* Calling DMPlexComputeGeometryFVM() generates the value returned by DMPlexGetMinRadius() */
83: DMPlexComputeGeometryFVM(celldm,&cellgeom,&facegeom);
84: DMPlexGetMinRadius(celldm,&minradius);
85: PetscPrintf(PETSC_COMM_WORLD,"PLEX(minradius) %1.4e\n",(double)minradius);
86: VecDestroy(&cellgeom);
87: VecDestroy(&facegeom);
88: PetscSectionDestroy(§ion);
89: }
91: /* Create the DMSwarm */
92: DMCreate(PETSC_COMM_WORLD,&swarm);
93: DMSetType(swarm,DMSWARM);
94: DMSetDimension(swarm,dim);
96: /* Configure swarm to be of type PIC */
97: DMSwarmSetType(swarm,DMSWARM_PIC);
98: DMSwarmSetCellDM(swarm,celldm);
100: /* Register two scalar fields within the DMSwarm */
101: DMSwarmRegisterPetscDatatypeField(swarm,"phi",1,PETSC_REAL);
102: DMSwarmRegisterPetscDatatypeField(swarm,"region",1,PETSC_REAL);
103: DMSwarmFinalizeFieldRegister(swarm);
105: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
106: DMSwarmSetLocalSizes(swarm,4,0);
108: /* Insert swarm coordinates cell-wise */
109: /*DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,ppcell);*/
110: DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_SUBDIVISION,ppcell);
112: /* Define initial conditions for th swarm fields "phi" and "region" */
113: {
114: PetscReal *s_coor,*s_phi,*s_region;
115: PetscInt npoints,p;
117: DMSwarmGetLocalSize(swarm,&npoints);
118: DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);
119: DMSwarmGetField(swarm,"phi",NULL,NULL,(void**)&s_phi);
120: 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];
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: DMSwarmRestoreField(swarm,"region",NULL,NULL,(void**)&s_region);
130: DMSwarmRestoreField(swarm,"phi",NULL,NULL,(void**)&s_phi);
131: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);
132: }
134: /* Project initial value of phi onto the mesh */
135: DMSwarmProjectFields(swarm,1,fieldnames,&pfields,PETSC_FALSE);
137: if (view) {
138: /* View swarm all swarm fields using data type PETSC_REAL */
139: DMSwarmViewXDMF(swarm,"ic_dms.xmf");
141: /* View projected swarm field "phi" */
142: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
143: PetscViewerSetType(viewer,PETSCVIEWERVTK);
144: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
145: if (meshtype == 0) { /* DA */
146: PetscViewerFileSetName(viewer,"ic_dmda.vts");
147: VecView(pfields[0],viewer);
148: }
149: if (meshtype == 1) { /* PLEX */
150: PetscViewerFileSetName(viewer,"ic_dmplex.vtk");
151: DMView(celldm,viewer);
152: VecView(pfields[0],viewer);
153: }
154: PetscViewerDestroy(&viewer);
155: }
157: DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);
158: DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);
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;
165: PetscPrintf(PETSC_COMM_WORLD,"[step %D]\n",tk);
166: /* advect with analytic prescribed (constant) velocity field */
167: DMSwarmGetLocalSize(swarm,&npoints);
168: 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: DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);
175: DMSwarmMigrate(swarm,PETSC_TRUE);
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];
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: 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];
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: DMSwarmSetPointsUniformCoordinates(swarm,min,max,npoints_dir_y,ADD_VALUES);
199: }
201: /* Project swarm field "phi" onto the cell DM */
202: DMSwarmProjectFields(swarm,1,fieldnames,&pfields,PETSC_TRUE);
204: if (view) {
205: PetscViewer viewer;
206: char fname[PETSC_MAX_PATH_LEN];
208: /* View swarm fields */
209: PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dms.xmf",tk);
210: DMSwarmViewXDMF(swarm,fname);
212: /* View projected field */
213: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
214: PetscViewerSetType(viewer,PETSCVIEWERVTK);
215: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
217: if (meshtype == 0) { /* DA */
218: PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dmda.vts",tk);
219: PetscViewerFileSetName(viewer,fname);
220: VecView(pfields[0],viewer);
221: }
222: if (meshtype == 1) { /* PLEX */
223: PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dmplex.vtk",tk);
224: PetscViewerFileSetName(viewer,fname);
225: DMView(celldm,viewer);
226: VecView(pfields[0],viewer);
227: }
228: PetscViewerDestroy(&viewer);
229: }
231: }
232: VecDestroy(&pfields[0]);
233: PetscFree(pfields);
234: DMDestroy(&celldm);
235: DMDestroy(&swarm);
237: return(0);
238: }
240: int main(int argc,char **args)
241: {
243: PetscInt ppcell = 1;
244: PetscInt meshtype = 0;
246: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
247: PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL);
248: PetscOptionsGetInt(NULL,NULL,"-meshtype",&meshtype,NULL);
249: if (meshtype > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"-meshtype <value> must be 0 or 1");
251: pic_advect(ppcell,meshtype);
253: PetscFinalize();
254: return ierr;
255: }