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 %" PetscInt_FMT "]\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;
186: max[0] = 0.5 * dx + 31.0 * dx;
187: min[1] = 0.5 * dx;
188: 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;
197: max[0] = 0.5 * dx;
198: min[1] = 0.5 * dx;
199: max[1] = 0.5 * dx + 31.0 * dx;
200: DMSwarmSetPointsUniformCoordinates(swarm, min, max, npoints_dir_y, ADD_VALUES);
201: }
203: /* Project swarm field "phi" onto the cell DM */
204: DMSwarmProjectFields(swarm, 1, fieldnames, &pfields, PETSC_TRUE);
206: if (view) {
207: PetscViewer viewer;
208: char fname[PETSC_MAX_PATH_LEN];
210: /* View swarm fields */
211: PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dms.xmf", tk);
212: DMSwarmViewXDMF(swarm, fname);
214: /* View projected field */
215: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
216: PetscViewerSetType(viewer, PETSCVIEWERVTK);
217: PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
219: if (meshtype == 0) { /* DA */
220: PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dmda.vts", tk);
221: PetscViewerFileSetName(viewer, fname);
222: VecView(pfields[0], viewer);
223: }
224: if (meshtype == 1) { /* PLEX */
225: PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dmplex.vtk", tk);
226: PetscViewerFileSetName(viewer, fname);
227: DMView(celldm, viewer);
228: VecView(pfields[0], viewer);
229: }
230: PetscViewerDestroy(&viewer);
231: }
232: }
233: VecDestroy(&pfields[0]);
234: PetscFree(pfields);
235: DMDestroy(&celldm);
236: DMDestroy(&swarm);
238: return 0;
239: }
241: int main(int argc, char **args)
242: {
243: PetscInt ppcell = 1;
244: PetscInt meshtype = 0;
247: PetscInitialize(&argc, &args, (char *)0, help);
248: PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL);
249: PetscOptionsGetInt(NULL, NULL, "-meshtype", &meshtype, NULL);
252: pic_advect(ppcell, meshtype);
254: PetscFinalize();
255: return 0;
256: }