Actual source code: ex21.c
1: static char help[] = "DMSwarm-PIC demonstrator of advecting points within cell DM defined by a DA or PLEX object \n\
2: Options: \n\
3: -ppcell : Number of times to sub-divide the reference cell when layout the initial particle coordinates \n\
4: -meshtype : 0 ==> DA , 1 ==> PLEX \n\
5: -nt : Number of timestep to perform \n\
6: -view : Write out initial condition and time dependent data \n";
8: #include <petsc.h>
9: #include <petscdm.h>
10: #include <petscdmda.h>
11: #include <petscdmswarm.h>
13: PetscErrorCode pic_advect(PetscInt ppcell, PetscInt meshtype)
14: {
15: const PetscInt dim = 2;
16: DM celldm, swarm;
17: PetscInt tk, nt = 200;
18: PetscBool view = PETSC_FALSE;
19: Vec pfields[1];
20: PetscReal minradius;
21: PetscReal dt;
22: PetscReal vel[] = {1.0, 0.16};
23: const char *fieldnames[] = {"phi"};
24: PetscViewer viewer;
26: PetscFunctionBegin;
27: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL));
28: PetscCall(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: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mesh type: DMDA\n"));
37: nxy = 33;
38: 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: PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1));
42: PetscCall(DMSetFromOptions(celldm));
44: PetscCall(DMSetUp(celldm));
46: PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.5));
48: minradius = 1.0 / ((PetscReal)(nxy - 1));
50: PetscCall(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: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mesh type: DMPLEX\n"));
64: PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, faces, NULL, NULL, PETSC_TRUE, &celldm));
66: /* Distribute mesh over processes */
67: PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh));
68: if (distributedMesh) {
69: PetscCall(DMDestroy(&celldm));
70: celldm = distributedMesh;
71: }
73: PetscCall(DMSetFromOptions(celldm));
75: PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion));
76: PetscCall(DMSetLocalSection(celldm, section));
78: PetscCall(DMSetUp(celldm));
80: /* Calling DMPlexComputeGeometryFVM() generates the value returned by DMPlexGetMinRadius() */
81: PetscCall(DMPlexComputeGeometryFVM(celldm, &cellgeom, &facegeom));
82: PetscCall(DMPlexGetMinRadius(celldm, &minradius));
83: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "PLEX(minradius) %1.4e\n", (double)minradius));
84: PetscCall(VecDestroy(&cellgeom));
85: PetscCall(VecDestroy(&facegeom));
86: PetscCall(PetscSectionDestroy(§ion));
87: }
89: /* Create the DMSwarm */
90: PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
91: PetscCall(DMSetType(swarm, DMSWARM));
92: PetscCall(DMSetDimension(swarm, dim));
94: /* Configure swarm to be of type PIC */
95: PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
96: PetscCall(DMSwarmSetCellDM(swarm, celldm));
98: /* Register two scalar fields within the DMSwarm */
99: PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "phi", 1, PETSC_REAL));
100: PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "region", 1, PETSC_REAL));
101: PetscCall(DMSwarmFinalizeFieldRegister(swarm));
103: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
104: PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));
106: /* Insert swarm coordinates cell-wise */
107: /*PetscCall(DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,ppcell));*/
108: PetscCall(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: PetscCall(DMSwarmGetLocalSize(swarm, &npoints));
116: PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor));
117: PetscCall(DMSwarmGetField(swarm, "phi", NULL, NULL, (void **)&s_phi));
118: PetscCall(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: PetscCall(DMSwarmRestoreField(swarm, "region", NULL, NULL, (void **)&s_region));
128: PetscCall(DMSwarmRestoreField(swarm, "phi", NULL, NULL, (void **)&s_phi));
129: PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&s_coor));
130: }
132: /* Project initial value of phi onto the mesh */
133: PetscCall(DMCreateGlobalVector(celldm, &pfields[0]));
134: PetscCall(DMSwarmProjectFields(swarm, NULL, 1, fieldnames, pfields, SCATTER_FORWARD));
136: if (view) {
137: /* View swarm all swarm fields using data type PETSC_REAL */
138: PetscCall(DMSwarmViewXDMF(swarm, "ic_dms.xmf"));
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: }
156: PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD));
157: PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD));
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;
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));
174: PetscCall(DMSwarmMigrate(swarm, PETSC_TRUE));
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];
186: min[0] = 0.5 * dx;
187: max[0] = 0.5 * dx + 31.0 * dx;
188: min[1] = 0.5 * dx;
189: max[1] = 0.5 * dx;
190: PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, npoints_dir_x, ADD_VALUES));
191: }
192: if (tk % 2 == 0) {
193: PetscReal dx = 1.0 / 32.0;
194: PetscInt npoints_dir_y[] = {2, 31};
195: PetscReal min[2], max[2];
197: min[0] = 0.05 * dx;
198: max[0] = 0.5 * dx;
199: min[1] = 0.5 * dx;
200: max[1] = 0.5 * dx + 31.0 * dx;
201: PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, npoints_dir_y, ADD_VALUES));
202: }
204: /* Project swarm field "phi" onto the cell DM */
205: PetscCall(DMSwarmProjectFields(swarm, NULL, 1, fieldnames, pfields, SCATTER_FORWARD));
207: if (view) {
208: PetscViewer viewer;
209: char fname[PETSC_MAX_PATH_LEN];
211: /* View swarm fields */
212: PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dms.xmf", tk));
213: PetscCall(DMSwarmViewXDMF(swarm, fname));
215: /* View projected field */
216: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
217: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
218: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
220: if (meshtype == 0) { /* DA */
221: PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dmda.vts", tk));
222: PetscCall(PetscViewerFileSetName(viewer, fname));
223: PetscCall(VecView(pfields[0], viewer));
224: }
225: if (meshtype == 1) { /* PLEX */
226: PetscCall(PetscSNPrintf(fname, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_dmplex.vtk", tk));
227: PetscCall(PetscViewerFileSetName(viewer, fname));
228: PetscCall(DMView(celldm, viewer));
229: PetscCall(VecView(pfields[0], viewer));
230: }
231: PetscCall(PetscViewerDestroy(&viewer));
232: }
233: }
234: PetscCall(VecDestroy(&pfields[0]));
235: PetscCall(DMDestroy(&celldm));
236: PetscCall(DMDestroy(&swarm));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: int main(int argc, char **args)
241: {
242: PetscInt ppcell = 1;
243: PetscInt meshtype = 0;
245: PetscFunctionBeginUser;
246: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
247: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL));
248: PetscCall(PetscOptionsGetInt(NULL, NULL, "-meshtype", &meshtype, NULL));
249: PetscCheck(meshtype <= 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "-meshtype <value> must be 0 or 1");
251: PetscCall(pic_advect(ppcell, meshtype));
253: PetscCall(PetscFinalize());
254: return 0;
255: }