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";
7
8 #include <petsc.h>
9 #include <petscdm.h>
10 #include <petscdmda.h>
11 #include <petscdmswarm.h>
12
pic_advect(PetscInt ppcell,PetscInt meshtype)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;
25
26 PetscFunctionBegin;
27 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL));
28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL));
29
30 /* Create the background cell DM */
31 if (meshtype == 0) { /* DA */
32 PetscInt nxy;
33 PetscInt dof = 1;
34 PetscInt stencil_width = 1;
35
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));
39
40 PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1));
41
42 PetscCall(DMSetFromOptions(celldm));
43
44 PetscCall(DMSetUp(celldm));
45
46 PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.5));
47
48 minradius = 1.0 / ((PetscReal)(nxy - 1));
49
50 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "DA(minradius) %1.4e\n", (double)minradius));
51 }
52
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;
62
63 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mesh type: DMPLEX\n"));
64 PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, faces, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &celldm));
65
66 /* Distribute mesh over processes */
67 PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh));
68 if (distributedMesh) {
69 PetscCall(DMDestroy(&celldm));
70 celldm = distributedMesh;
71 }
72
73 PetscCall(DMSetFromOptions(celldm));
74
75 PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion));
76 PetscCall(DMSetLocalSection(celldm, section));
77
78 PetscCall(DMSetUp(celldm));
79
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 }
88
89 /* Create the DMSwarm */
90 PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
91 PetscCall(DMSetType(swarm, DMSWARM));
92 PetscCall(DMSetDimension(swarm, dim));
93
94 /* Configure swarm to be of type PIC */
95 PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
96 PetscCall(DMSwarmSetCellDM(swarm, celldm));
97
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));
102
103 /* Set initial local sizes of the DMSwarm with a buffer length of zero */
104 PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));
105
106 /* Insert swarm coordinates cell-wise */
107 /*PetscCall(DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,ppcell));*/
108 PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell));
109
110 /* Define initial conditions for th swarm fields "phi" and "region" */
111 {
112 PetscReal *s_coor, *s_phi, *s_region;
113 PetscInt npoints, p;
114
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];
123
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 }
131
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));
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;
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];
196
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 }
203
204 /* Project swarm field "phi" onto the cell DM */
205 PetscCall(DMSwarmProjectFields(swarm, NULL, 1, fieldnames, pfields, SCATTER_FORWARD));
206
207 if (view) {
208 PetscViewer viewer;
209 char fname[PETSC_MAX_PATH_LEN];
210
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));
214
215 /* View projected field */
216 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
217 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
218 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
219
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 }
239
main(int argc,char ** args)240 int main(int argc, char **args)
241 {
242 PetscInt ppcell = 1;
243 PetscInt meshtype = 0;
244
245 PetscFunctionBeginUser;
246 PetscCall(PetscInitialize(&argc, &args, NULL, 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");
250
251 PetscCall(pic_advect(ppcell, meshtype));
252
253 PetscCall(PetscFinalize());
254 return 0;
255 }
256