Actual source code: swarmpic.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petscsf.h>
4: #include <petscdmda.h>
5: #include <petscdmplex.h>
6: #include <petscdt.h>
7: #include "../src/dm/impls/swarm/data_bucket.h"
9: #include <petsc/private/petscfeimpl.h>
11: /* Coordinate insertition/addition API */
12: /*@
13: DMSwarmSetPointsUniformCoordinates - Set point coordinates in a `DMSWARM` on a regular (ijk) grid
15: Collective
17: Input Parameters:
18: + dm - the `DMSWARM`
19: . min - minimum coordinate values in the x, y, z directions (array of length dim)
20: . max - maximum coordinate values in the x, y, z directions (array of length dim)
21: . npoints - number of points in each spatial direction (array of length dim)
22: - mode - indicates whether to append points to the swarm (`ADD_VALUES`), or over-ride existing points (`INSERT_VALUES`)
24: Level: beginner
26: Notes:
27: When using mode = `INSERT_VALUES`, this method will reset the number of particles in the `DMSWARM`
28: to be `npoints[0]` x `npoints[1]` (2D) or `npoints[0]` x `npoints[1]` x `npoints[2]` (3D). When using mode = `ADD_VALUES`,
29: new points will be appended to any already existing in the `DMSWARM`
31: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
32: @*/
33: PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM dm, PetscReal min[], PetscReal max[], PetscInt npoints[], InsertMode mode)
34: {
35: PetscReal lmin[] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
36: PetscReal lmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
37: PetscInt i, j, k, bs, b, n_estimate, n_curr, n_new_est, p, n_found;
38: const PetscScalar *_coor;
39: DM celldm;
40: PetscReal dx[3];
41: PetscInt _npoints[] = {0, 0, 1};
42: Vec pos;
43: PetscScalar *_pos;
44: PetscReal *swarm_coor;
45: PetscInt *swarm_cellid;
46: PetscSF sfcell = NULL;
47: const PetscSFNode *LA_sfcell;
48: const char *coordname;
50: PetscFunctionBegin;
51: DMSWARMPICVALID(dm);
52: PetscCall(DMSwarmGetCellDM(dm, &celldm));
53: PetscCall(DMGetLocalBoundingBox(celldm, lmin, lmax));
54: PetscCall(DMGetCoordinateDim(celldm, &bs));
55: PetscCall(DMSwarmGetCoordinateField(dm, &coordname));
57: for (b = 0; b < bs; b++) {
58: if (npoints[b] > 1) {
59: dx[b] = (max[b] - min[b]) / ((PetscReal)(npoints[b] - 1));
60: } else {
61: dx[b] = 0.0;
62: }
63: _npoints[b] = npoints[b];
64: }
66: /* determine number of points living in the bounding box */
67: n_estimate = 0;
68: for (k = 0; k < _npoints[2]; k++) {
69: for (j = 0; j < _npoints[1]; j++) {
70: for (i = 0; i < _npoints[0]; i++) {
71: PetscReal xp[] = {0.0, 0.0, 0.0};
72: PetscInt ijk[3];
73: PetscBool point_inside = PETSC_TRUE;
75: ijk[0] = i;
76: ijk[1] = j;
77: ijk[2] = k;
78: for (b = 0; b < bs; b++) xp[b] = min[b] + ijk[b] * dx[b];
79: for (b = 0; b < bs; b++) {
80: if (xp[b] < lmin[b]) point_inside = PETSC_FALSE;
81: if (xp[b] > lmax[b]) point_inside = PETSC_FALSE;
82: }
83: if (point_inside) n_estimate++;
84: }
85: }
86: }
88: /* create candidate list */
89: PetscCall(VecCreate(PETSC_COMM_SELF, &pos));
90: PetscCall(VecSetSizes(pos, bs * n_estimate, PETSC_DECIDE));
91: PetscCall(VecSetBlockSize(pos, bs));
92: PetscCall(VecSetFromOptions(pos));
93: PetscCall(VecGetArray(pos, &_pos));
95: n_estimate = 0;
96: for (k = 0; k < _npoints[2]; k++) {
97: for (j = 0; j < _npoints[1]; j++) {
98: for (i = 0; i < _npoints[0]; i++) {
99: PetscReal xp[] = {0.0, 0.0, 0.0};
100: PetscInt ijk[3];
101: PetscBool point_inside = PETSC_TRUE;
103: ijk[0] = i;
104: ijk[1] = j;
105: ijk[2] = k;
106: for (b = 0; b < bs; b++) xp[b] = min[b] + ijk[b] * dx[b];
107: for (b = 0; b < bs; b++) {
108: if (xp[b] < lmin[b]) point_inside = PETSC_FALSE;
109: if (xp[b] > lmax[b]) point_inside = PETSC_FALSE;
110: }
111: if (point_inside) {
112: for (b = 0; b < bs; b++) _pos[bs * n_estimate + b] = xp[b];
113: n_estimate++;
114: }
115: }
116: }
117: }
118: PetscCall(VecRestoreArray(pos, &_pos));
120: /* locate points */
121: PetscCall(DMLocatePoints(celldm, pos, DM_POINTLOCATION_NONE, &sfcell));
122: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
123: n_found = 0;
124: for (p = 0; p < n_estimate; p++) {
125: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) n_found++;
126: }
128: /* adjust size */
129: if (mode == ADD_VALUES) {
130: PetscCall(DMSwarmGetLocalSize(dm, &n_curr));
131: n_new_est = n_curr + n_found;
132: PetscCall(DMSwarmSetLocalSizes(dm, n_new_est, -1));
133: }
134: if (mode == INSERT_VALUES) {
135: n_curr = 0;
136: n_new_est = n_found;
137: PetscCall(DMSwarmSetLocalSizes(dm, n_new_est, -1));
138: }
140: /* initialize new coords, cell owners, pid */
141: PetscCall(VecGetArrayRead(pos, &_coor));
142: PetscCall(DMSwarmGetField(dm, coordname, NULL, NULL, (void **)&swarm_coor));
143: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
144: n_found = 0;
145: for (p = 0; p < n_estimate; p++) {
146: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
147: for (b = 0; b < bs; b++) swarm_coor[bs * (n_curr + n_found) + b] = PetscRealPart(_coor[bs * p + b]);
148: swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
149: n_found++;
150: }
151: }
152: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
153: PetscCall(DMSwarmRestoreField(dm, coordname, NULL, NULL, (void **)&swarm_coor));
154: PetscCall(VecRestoreArrayRead(pos, &_coor));
156: PetscCall(PetscSFDestroy(&sfcell));
157: PetscCall(VecDestroy(&pos));
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: /*@
162: DMSwarmSetPointCoordinates - Set point coordinates in a `DMSWARM` from a user defined list
164: Collective
166: Input Parameters:
167: + dm - the `DMSWARM`
168: . npoints - the number of points to insert
169: . coor - the coordinate values
170: . redundant - if set to `PETSC_TRUE`, it is assumed that `npoints` and `coor` are only valid on rank 0 and should be broadcast to other ranks
171: - mode - indicates whether to append points to the swarm (`ADD_VALUES`), or over-ride existing points (`INSERT_VALUES`)
173: Level: beginner
175: Notes:
176: If the user has specified `redundant` as `PETSC_FALSE`, the cell `DM` will attempt to locate the coordinates provided by `coor` within
177: its sub-domain. If they any values within `coor` are not located in the sub-domain, they will be ignored and will not get
178: added to the `DMSWARM`.
180: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`, `DMSwarmSetPointsUniformCoordinates()`
181: @*/
182: PetscErrorCode DMSwarmSetPointCoordinates(DM dm, PetscInt npoints, PetscReal coor[], PetscBool redundant, InsertMode mode)
183: {
184: PetscReal gmin[] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
185: PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
186: PetscInt i, N, bs, b, n_estimate, n_curr, n_new_est, p, n_found;
187: Vec coorlocal;
188: const PetscScalar *_coor;
189: DM celldm;
190: Vec pos;
191: PetscScalar *_pos;
192: PetscReal *swarm_coor;
193: PetscInt *swarm_cellid;
194: PetscSF sfcell = NULL;
195: const PetscSFNode *LA_sfcell;
196: PetscReal *my_coor;
197: PetscInt my_npoints;
198: PetscMPIInt rank;
199: MPI_Comm comm;
200: const char *coordname;
202: PetscFunctionBegin;
203: DMSWARMPICVALID(dm);
204: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
205: PetscCallMPI(MPI_Comm_rank(comm, &rank));
207: PetscCall(DMSwarmGetCellDM(dm, &celldm));
208: PetscCall(DMSwarmGetCoordinateField(dm, &coordname));
209: PetscCall(DMGetCoordinatesLocal(celldm, &coorlocal));
210: PetscCall(VecGetSize(coorlocal, &N));
211: PetscCall(VecGetBlockSize(coorlocal, &bs));
212: N = N / bs;
213: PetscCall(VecGetArrayRead(coorlocal, &_coor));
214: for (i = 0; i < N; i++) {
215: for (b = 0; b < bs; b++) {
216: gmin[b] = PetscMin(gmin[b], PetscRealPart(_coor[bs * i + b]));
217: gmax[b] = PetscMax(gmax[b], PetscRealPart(_coor[bs * i + b]));
218: }
219: }
220: PetscCall(VecRestoreArrayRead(coorlocal, &_coor));
222: /* broadcast points from rank 0 if requested */
223: if (redundant) {
224: PetscMPIInt imy;
226: my_npoints = npoints;
227: PetscCallMPI(MPI_Bcast(&my_npoints, 1, MPIU_INT, 0, comm));
229: if (rank > 0) { /* allocate space */
230: PetscCall(PetscMalloc1(bs * my_npoints, &my_coor));
231: } else {
232: my_coor = coor;
233: }
234: PetscCall(PetscMPIIntCast(bs * my_npoints, &imy));
235: PetscCallMPI(MPI_Bcast(my_coor, imy, MPIU_REAL, 0, comm));
236: } else {
237: my_npoints = npoints;
238: my_coor = coor;
239: }
241: /* determine the number of points living in the bounding box */
242: n_estimate = 0;
243: for (i = 0; i < my_npoints; i++) {
244: PetscBool point_inside = PETSC_TRUE;
246: for (b = 0; b < bs; b++) {
247: if (my_coor[bs * i + b] < gmin[b]) point_inside = PETSC_FALSE;
248: if (my_coor[bs * i + b] > gmax[b]) point_inside = PETSC_FALSE;
249: }
250: if (point_inside) n_estimate++;
251: }
253: /* create candidate list */
254: PetscCall(VecCreate(PETSC_COMM_SELF, &pos));
255: PetscCall(VecSetSizes(pos, bs * n_estimate, PETSC_DECIDE));
256: PetscCall(VecSetBlockSize(pos, bs));
257: PetscCall(VecSetFromOptions(pos));
258: PetscCall(VecGetArray(pos, &_pos));
260: n_estimate = 0;
261: for (i = 0; i < my_npoints; i++) {
262: PetscBool point_inside = PETSC_TRUE;
264: for (b = 0; b < bs; b++) {
265: if (my_coor[bs * i + b] < gmin[b]) point_inside = PETSC_FALSE;
266: if (my_coor[bs * i + b] > gmax[b]) point_inside = PETSC_FALSE;
267: }
268: if (point_inside) {
269: for (b = 0; b < bs; b++) _pos[bs * n_estimate + b] = my_coor[bs * i + b];
270: n_estimate++;
271: }
272: }
273: PetscCall(VecRestoreArray(pos, &_pos));
275: /* locate points */
276: PetscCall(DMLocatePoints(celldm, pos, DM_POINTLOCATION_NONE, &sfcell));
278: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
279: n_found = 0;
280: for (p = 0; p < n_estimate; p++) {
281: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) n_found++;
282: }
284: /* adjust size */
285: if (mode == ADD_VALUES) {
286: PetscCall(DMSwarmGetLocalSize(dm, &n_curr));
287: n_new_est = n_curr + n_found;
288: PetscCall(DMSwarmSetLocalSizes(dm, n_new_est, -1));
289: }
290: if (mode == INSERT_VALUES) {
291: n_curr = 0;
292: n_new_est = n_found;
293: PetscCall(DMSwarmSetLocalSizes(dm, n_new_est, -1));
294: }
296: /* initialize new coords, cell owners, pid */
297: PetscCall(VecGetArrayRead(pos, &_coor));
298: PetscCall(DMSwarmGetField(dm, coordname, NULL, NULL, (void **)&swarm_coor));
299: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
300: n_found = 0;
301: for (p = 0; p < n_estimate; p++) {
302: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
303: for (b = 0; b < bs; b++) swarm_coor[bs * (n_curr + n_found) + b] = PetscRealPart(_coor[bs * p + b]);
304: swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
305: n_found++;
306: }
307: }
308: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
309: PetscCall(DMSwarmRestoreField(dm, coordname, NULL, NULL, (void **)&swarm_coor));
310: PetscCall(VecRestoreArrayRead(pos, &_coor));
312: if (redundant) {
313: if (rank > 0) PetscCall(PetscFree(my_coor));
314: }
315: PetscCall(PetscSFDestroy(&sfcell));
316: PetscCall(VecDestroy(&pos));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM, DM, DMSwarmPICLayoutType, PetscInt);
321: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM, DM, DMSwarmPICLayoutType, PetscInt);
323: /*@
324: DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell
326: Not Collective
328: Input Parameters:
329: + dm - the `DMSWARM`
330: . layout_type - method used to fill each cell with the cell `DM`
331: - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type)
333: Level: beginner
335: Notes:
336: The insert method will reset any previous defined points within the `DMSWARM`.
338: When using a `DMDA` both 2D and 3D are supported for all layout types provided you are using `DMDA_ELEMENT_Q1`.
340: When using a `DMPLEX` the following case are supported\:
341: .vb
342: (i) DMSWARMPIC_LAYOUT_REGULAR: 2D (triangle),
343: (ii) DMSWARMPIC_LAYOUT_GAUSS: 2D and 3D provided the cell is a tri/tet or a quad/hex,
344: (iii) DMSWARMPIC_LAYOUT_SUBDIVISION: 2D and 3D for quad/hex and 2D tri.
345: .ve
347: .seealso: `DMSWARM`, `DMSwarmPICLayoutType`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
348: @*/
349: PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm, DMSwarmPICLayoutType layout_type, PetscInt fill_param)
350: {
351: DM celldm;
352: PetscBool isDA, isPLEX;
354: PetscFunctionBegin;
355: DMSWARMPICVALID(dm);
356: PetscCall(DMSwarmGetCellDM(dm, &celldm));
357: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isDA));
358: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isPLEX));
359: if (isDA) {
360: PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA(dm, celldm, layout_type, fill_param));
361: } else if (isPLEX) {
362: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX(dm, celldm, layout_type, fill_param));
363: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: extern PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *);
369: /*@C
370: DMSwarmSetPointCoordinatesCellwise - Insert point coordinates (defined over the reference cell) within each cell
372: Not Collective
374: Input Parameters:
375: + dm - the `DMSWARM`
376: . npoints - the number of points to insert in each cell
377: - xi - the coordinates (defined in the local coordinate system for each cell) to insert
379: Level: beginner
381: Notes:
382: The method will reset any previous defined points within the `DMSWARM`.
383: Only supported for `DMPLEX`. If you are using a `DMDA` it is recommended to either use
384: `DMSwarmInsertPointsUsingCellDM()`, or extract and set the coordinates yourself the following code
385: .vb
386: PetscReal *coor;
387: const char *coordname;
388: DMSwarmGetCoordinateField(dm, &coordname);
389: DMSwarmGetField(dm,coordname,NULL,NULL,(void**)&coor);
390: // user code to define the coordinates here
391: DMSwarmRestoreField(dm,coordname,NULL,NULL,(void**)&coor);
392: .ve
394: .seealso: `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmInsertPointsUsingCellDM()`
395: @*/
396: PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm, PetscInt npoints, PetscReal xi[])
397: {
398: DM celldm;
399: PetscBool isDA, isPLEX;
401: PetscFunctionBegin;
402: DMSWARMPICVALID(dm);
403: PetscCall(DMSwarmGetCellDM(dm, &celldm));
404: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isDA));
405: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isPLEX));
406: PetscCheck(!isDA, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMPLEX. Recommended you use DMSwarmInsertPointsUsingCellDM()");
407: if (isPLEX) {
408: PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, npoints, xi));
409: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*@
414: DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM
416: Not Collective
418: Input Parameter:
419: . dm - the `DMSWARM`
421: Output Parameters:
422: + ncells - the number of cells in the cell `DM` (optional argument, pass `NULL` to ignore)
423: - count - array of length ncells containing the number of points per cell
425: Level: beginner
427: Notes:
428: The array count is allocated internally and must be free'd by the user.
430: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
431: @*/
432: PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm, PetscInt *ncells, PetscInt **count)
433: {
434: PetscBool isvalid;
435: PetscInt nel;
436: PetscInt *sum;
438: PetscFunctionBegin;
439: PetscCall(DMSwarmSortGetIsValid(dm, &isvalid));
440: nel = 0;
441: if (isvalid) {
442: PetscInt e;
444: PetscCall(DMSwarmSortGetSizes(dm, &nel, NULL));
446: PetscCall(PetscMalloc1(nel, &sum));
447: for (e = 0; e < nel; e++) PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm, e, &sum[e]));
448: } else {
449: DM celldm;
450: PetscBool isda, isplex, isshell;
451: PetscInt p, npoints;
452: PetscInt *swarm_cellid;
454: /* get the number of cells */
455: PetscCall(DMSwarmGetCellDM(dm, &celldm));
456: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isda));
457: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isplex));
458: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMSHELL, &isshell));
459: if (isda) {
460: PetscInt _nel, _npe;
461: const PetscInt *_element;
463: PetscCall(DMDAGetElements(celldm, &_nel, &_npe, &_element));
464: nel = _nel;
465: PetscCall(DMDARestoreElements(celldm, &_nel, &_npe, &_element));
466: } else if (isplex) {
467: PetscInt ps, pe;
469: PetscCall(DMPlexGetHeightStratum(celldm, 0, &ps, &pe));
470: nel = pe - ps;
471: } else if (isshell) {
472: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
474: PetscCall(PetscObjectQueryFunction((PetscObject)celldm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
475: if (method_DMShellGetNumberOfCells) {
476: PetscCall(method_DMShellGetNumberOfCells(celldm, &nel));
477: } else
478: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine the number of cells for the DMSHELL object. User must provide a method via PetscObjectComposeFunction( (PetscObject)shelldm, \"DMGetNumberOfCells_C\", your_function_to_compute_number_of_cells);");
479: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
481: PetscCall(PetscMalloc1(nel, &sum));
482: PetscCall(PetscArrayzero(sum, nel));
483: PetscCall(DMSwarmGetLocalSize(dm, &npoints));
484: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
485: for (p = 0; p < npoints; p++) {
486: if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) sum[swarm_cellid[p]]++;
487: }
488: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
489: }
490: if (ncells) *ncells = nel;
491: *count = sum;
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@
496: DMSwarmGetNumSpecies - Get the number of particle species
498: Not Collective
500: Input Parameter:
501: . sw - the `DMSWARM`
503: Output Parameters:
504: . Ns - the number of species
506: Level: intermediate
508: .seealso: `DMSWARM`, `DMSwarmSetNumSpecies()`, `DMSwarmSetType()`, `DMSwarmType`
509: @*/
510: PetscErrorCode DMSwarmGetNumSpecies(DM sw, PetscInt *Ns)
511: {
512: DM_Swarm *swarm = (DM_Swarm *)sw->data;
514: PetscFunctionBegin;
515: *Ns = swarm->Ns;
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: /*@
520: DMSwarmSetNumSpecies - Set the number of particle species
522: Not Collective
524: Input Parameters:
525: + sw - the `DMSWARM`
526: - Ns - the number of species
528: Level: intermediate
530: .seealso: `DMSWARM`, `DMSwarmGetNumSpecies()`, `DMSwarmSetType()`, `DMSwarmType`
531: @*/
532: PetscErrorCode DMSwarmSetNumSpecies(DM sw, PetscInt Ns)
533: {
534: DM_Swarm *swarm = (DM_Swarm *)sw->data;
536: PetscFunctionBegin;
537: swarm->Ns = Ns;
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: /*@C
542: DMSwarmGetCoordinateFunction - Get the function setting initial particle positions, if it exists
544: Not Collective
546: Input Parameter:
547: . sw - the `DMSWARM`
549: Output Parameter:
550: . coordFunc - the function setting initial particle positions, or `NULL`, see `PetscSimplePointFn` for the calling sequence
552: Level: intermediate
554: .seealso: `DMSWARM`, `DMSwarmSetCoordinateFunction()`, `DMSwarmGetVelocityFunction()`, `DMSwarmInitializeCoordinates()`, `PetscSimplePointFn`
555: @*/
556: PetscErrorCode DMSwarmGetCoordinateFunction(DM sw, PetscSimplePointFn **coordFunc)
557: {
558: DM_Swarm *swarm = (DM_Swarm *)sw->data;
560: PetscFunctionBegin;
562: PetscAssertPointer(coordFunc, 2);
563: *coordFunc = swarm->coordFunc;
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: /*@C
568: DMSwarmSetCoordinateFunction - Set the function setting initial particle positions
570: Not Collective
572: Input Parameters:
573: + sw - the `DMSWARM`
574: - coordFunc - the function setting initial particle positions, see `PetscSimplePointFn` for the calling sequence
576: Level: intermediate
578: .seealso: `DMSWARM`, `DMSwarmGetCoordinateFunction()`, `DMSwarmSetVelocityFunction()`, `DMSwarmInitializeCoordinates()`, `PetscSimplePointFn`
579: @*/
580: PetscErrorCode DMSwarmSetCoordinateFunction(DM sw, PetscSimplePointFn *coordFunc)
581: {
582: DM_Swarm *swarm = (DM_Swarm *)sw->data;
584: PetscFunctionBegin;
587: swarm->coordFunc = coordFunc;
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: /*@C
592: DMSwarmGetVelocityFunction - Get the function setting initial particle velocities, if it exists
594: Not Collective
596: Input Parameter:
597: . sw - the `DMSWARM`
599: Output Parameter:
600: . velFunc - the function setting initial particle velocities, or `NULL`, see `PetscSimplePointFn` for the calling sequence
602: Level: intermediate
604: .seealso: `DMSWARM`, `DMSwarmSetVelocityFunction()`, `DMSwarmGetCoordinateFunction()`, `DMSwarmInitializeVelocities()`, `PetscSimplePointFn`
605: @*/
606: PetscErrorCode DMSwarmGetVelocityFunction(DM sw, PetscSimplePointFn **velFunc)
607: {
608: DM_Swarm *swarm = (DM_Swarm *)sw->data;
610: PetscFunctionBegin;
612: PetscAssertPointer(velFunc, 2);
613: *velFunc = swarm->velFunc;
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: /*@C
618: DMSwarmSetVelocityFunction - Set the function setting initial particle velocities
620: Not Collective
622: Input Parameters:
623: + sw - the `DMSWARM`
624: - velFunc - the function setting initial particle velocities, see `PetscSimplePointFn` for the calling sequence
626: Level: intermediate
628: .seealso: `DMSWARM`, `DMSwarmGetVelocityFunction()`, `DMSwarmSetCoordinateFunction()`, `DMSwarmInitializeVelocities()`, `PetscSimplePointFn`
629: @*/
630: PetscErrorCode DMSwarmSetVelocityFunction(DM sw, PetscSimplePointFn *velFunc)
631: {
632: DM_Swarm *swarm = (DM_Swarm *)sw->data;
634: PetscFunctionBegin;
637: swarm->velFunc = velFunc;
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: /*@C
642: DMSwarmComputeLocalSize - Compute the local number and distribution of particles based upon a density function
644: Not Collective
646: Input Parameters:
647: + sw - The `DMSWARM`
648: . N - The target number of particles
649: - density - The density field for the particle layout, normalized to unity
651: Level: advanced
653: Note:
654: One particle will be created for each species.
656: .seealso: `DMSWARM`, `DMSwarmComputeLocalSizeFromOptions()`
657: @*/
658: PetscErrorCode DMSwarmComputeLocalSize(DM sw, PetscInt N, PetscProbFunc density)
659: {
660: DM dm;
661: PetscQuadrature quad;
662: const PetscReal *xq, *wq;
663: PetscReal *n_int;
664: PetscInt *npc_s, *cellid, Ni;
665: PetscReal gmin[3], gmax[3], xi0[3];
666: PetscInt Ns, cStart, cEnd, c, dim, d, Nq, q, Np = 0, p, s;
667: PetscBool simplex;
669: PetscFunctionBegin;
670: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
671: PetscCall(DMSwarmGetCellDM(sw, &dm));
672: PetscCall(DMGetDimension(dm, &dim));
673: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
674: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
675: PetscCall(DMPlexIsSimplex(dm, &simplex));
676: PetscCall(DMGetCoordinatesLocalSetUp(dm));
677: if (simplex) PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 5, -1.0, 1.0, &quad));
678: else PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 5, -1.0, 1.0, &quad));
679: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xq, &wq));
680: PetscCall(PetscCalloc2(Ns, &n_int, (cEnd - cStart) * Ns, &npc_s));
681: /* Integrate the density function to get the number of particles in each cell */
682: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
683: for (c = 0; c < cEnd - cStart; ++c) {
684: const PetscInt cell = c + cStart;
685: PetscReal v0[3], J[9], invJ[9], detJ, detJp = 2. / (gmax[0] - gmin[0]), xr[3], den;
687: /* Have to transform quadrature points/weights to cell domain */
688: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
689: PetscCall(PetscArrayzero(n_int, Ns));
690: for (q = 0; q < Nq; ++q) {
691: CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], xr);
692: /* Have to transform mesh to domain of definition of PDF, [-1, 1], and weight PDF by |J|/2 */
693: xr[0] = detJp * (xr[0] - gmin[0]) - 1.;
695: for (s = 0; s < Ns; ++s) {
696: PetscCall(density(xr, NULL, &den));
697: n_int[s] += (detJp * den) * (detJ * wq[q]) / (PetscReal)Ns;
698: }
699: }
700: for (s = 0; s < Ns; ++s) {
701: Ni = N;
702: npc_s[c * Ns + s] += (PetscInt)(Ni * n_int[s]);
703: Np += npc_s[c * Ns + s];
704: }
705: }
706: PetscCall(PetscQuadratureDestroy(&quad));
707: PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
708: PetscCall(DMSwarmGetField(sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
709: for (c = 0, p = 0; c < cEnd - cStart; ++c) {
710: for (s = 0; s < Ns; ++s) {
711: for (q = 0; q < npc_s[c * Ns + s]; ++q, ++p) cellid[p] = c;
712: }
713: }
714: PetscCall(DMSwarmRestoreField(sw, DMSwarmPICField_cellid, NULL, NULL, (void **)&cellid));
715: PetscCall(PetscFree2(n_int, npc_s));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: /*@
720: DMSwarmComputeLocalSizeFromOptions - Compute the local number and distribution of particles based upon a density function determined by options
722: Not Collective
724: Input Parameter:
725: . sw - The `DMSWARM`
727: Level: advanced
729: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`
730: @*/
731: PetscErrorCode DMSwarmComputeLocalSizeFromOptions(DM sw)
732: {
733: PetscProbFunc pdf;
734: const char *prefix;
735: char funcname[PETSC_MAX_PATH_LEN];
736: PetscInt *N, Ns, dim, n;
737: PetscBool flg;
738: PetscMPIInt size, rank;
740: PetscFunctionBegin;
741: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sw), &size));
742: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
743: PetscCall(PetscCalloc1(size, &N));
744: PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
745: n = size;
746: PetscCall(PetscOptionsIntArray("-dm_swarm_num_particles", "The target number of particles", "", N, &n, NULL));
747: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
748: PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
749: if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
750: PetscCall(PetscOptionsString("-dm_swarm_coordinate_function", "Function to determine particle coordinates", "DMSwarmSetCoordinateFunction", funcname, funcname, sizeof(funcname), &flg));
751: PetscOptionsEnd();
752: if (flg) {
753: PetscSimplePointFn *coordFunc;
755: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
756: PetscCall(PetscDLSym(NULL, funcname, (void **)&coordFunc));
757: PetscCheck(coordFunc, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
758: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
759: PetscCall(DMSwarmSetLocalSizes(sw, N[rank] * Ns, 0));
760: PetscCall(DMSwarmSetCoordinateFunction(sw, coordFunc));
761: } else {
762: PetscCall(DMGetDimension(sw, &dim));
763: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
764: PetscCall(PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_coordinate_density", &pdf, NULL, NULL));
765: PetscCall(DMSwarmComputeLocalSize(sw, N[rank], pdf));
766: }
767: PetscCall(PetscFree(N));
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: /*@
772: DMSwarmInitializeCoordinates - Determine the initial coordinates of particles for a PIC method
774: Not Collective
776: Input Parameter:
777: . sw - The `DMSWARM`
779: Level: advanced
781: Note:
782: Currently, we randomly place particles in their assigned cell
784: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeVelocities()`
785: @*/
786: PetscErrorCode DMSwarmInitializeCoordinates(DM sw)
787: {
788: PetscSimplePointFn *coordFunc;
789: PetscScalar *weight;
790: PetscReal *x;
791: PetscInt *species;
792: void *ctx;
793: PetscBool removePoints = PETSC_TRUE;
794: PetscDataType dtype;
795: PetscInt Np, p, Ns, dim, d, bs;
796: const char *coordname;
798: PetscFunctionBeginUser;
799: PetscCall(DMGetDimension(sw, &dim));
800: PetscCall(DMSwarmGetLocalSize(sw, &Np));
801: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
802: PetscCall(DMSwarmGetCoordinateField(sw, &coordname));
803: PetscCall(DMSwarmGetCoordinateFunction(sw, &coordFunc));
805: PetscCall(DMSwarmGetField(sw, coordname, &bs, &dtype, (void **)&x));
806: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&weight));
807: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
808: if (coordFunc) {
809: PetscCall(DMGetApplicationContext(sw, &ctx));
810: for (p = 0; p < Np; ++p) {
811: PetscScalar X[3];
813: PetscCall((*coordFunc)(dim, 0., NULL, p, X, ctx));
814: for (d = 0; d < dim; ++d) x[p * dim + d] = PetscRealPart(X[d]);
815: weight[p] = 1.0;
816: species[p] = p % Ns;
817: }
818: } else {
819: DM dm;
820: PetscRandom rnd;
821: PetscReal xi0[3];
822: PetscInt cStart, cEnd, c;
824: PetscCall(DMSwarmGetCellDM(sw, &dm));
825: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
826: PetscCall(DMGetApplicationContext(sw, &ctx));
828: /* Set particle position randomly in cell, set weights to 1 */
829: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
830: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
831: PetscCall(PetscRandomSetFromOptions(rnd));
832: PetscCall(DMSwarmSortGetAccess(sw));
833: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
834: for (c = cStart; c < cEnd; ++c) {
835: PetscReal v0[3], J[9], invJ[9], detJ;
836: PetscInt *pidx, Npc, q;
838: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
839: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
840: for (q = 0; q < Npc; ++q) {
841: const PetscInt p = pidx[q];
842: PetscReal xref[3];
844: for (d = 0; d < dim; ++d) PetscCall(PetscRandomGetValueReal(rnd, &xref[d]));
845: CoordinatesRefToReal(dim, dim, xi0, v0, J, xref, &x[p * dim]);
847: weight[p] = 1.0 / Np;
848: species[p] = p % Ns;
849: }
850: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
851: }
852: PetscCall(PetscRandomDestroy(&rnd));
853: PetscCall(DMSwarmSortRestoreAccess(sw));
854: }
855: PetscCall(DMSwarmRestoreField(sw, coordname, NULL, NULL, (void **)&x));
856: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
857: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
859: PetscCall(DMSwarmMigrate(sw, removePoints));
860: PetscCall(DMLocalizeCoordinates(sw));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: /*@C
865: DMSwarmInitializeVelocities - Set the initial velocities of particles using a distribution.
867: Collective
869: Input Parameters:
870: + sw - The `DMSWARM` object
871: . sampler - A function which uniformly samples the velocity PDF
872: - v0 - The velocity scale for nondimensionalization for each species
874: Level: advanced
876: Note:
877: If `v0` is zero for the first species, all velocities are set to zero. If it is zero for any other species, the effect will be to give that species zero velocity.
879: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeCoordinates()`, `DMSwarmInitializeVelocitiesFromOptions()`
880: @*/
881: PetscErrorCode DMSwarmInitializeVelocities(DM sw, PetscProbFunc sampler, const PetscReal v0[])
882: {
883: PetscSimplePointFn *velFunc;
884: PetscReal *v;
885: PetscInt *species;
886: void *ctx;
887: PetscInt dim, Np, p;
889: PetscFunctionBegin;
890: PetscCall(DMSwarmGetVelocityFunction(sw, &velFunc));
892: PetscCall(DMGetDimension(sw, &dim));
893: PetscCall(DMSwarmGetLocalSize(sw, &Np));
894: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
895: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
896: if (v0[0] == 0.) {
897: PetscCall(PetscArrayzero(v, Np * dim));
898: } else if (velFunc) {
899: PetscCall(DMGetApplicationContext(sw, &ctx));
900: for (p = 0; p < Np; ++p) {
901: PetscInt s = species[p], d;
902: PetscScalar vel[3];
904: PetscCall((*velFunc)(dim, 0., NULL, p, vel, ctx));
905: for (d = 0; d < dim; ++d) v[p * dim + d] = (v0[s] / v0[0]) * PetscRealPart(vel[d]);
906: }
907: } else {
908: PetscRandom rnd;
910: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd));
911: PetscCall(PetscRandomSetInterval(rnd, 0, 1.));
912: PetscCall(PetscRandomSetFromOptions(rnd));
914: for (p = 0; p < Np; ++p) {
915: PetscInt s = species[p], d;
916: PetscReal a[3], vel[3];
918: for (d = 0; d < dim; ++d) PetscCall(PetscRandomGetValueReal(rnd, &a[d]));
919: PetscCall(sampler(a, NULL, vel));
920: for (d = 0; d < dim; ++d) v[p * dim + d] = (v0[s] / v0[0]) * vel[d];
921: }
922: PetscCall(PetscRandomDestroy(&rnd));
923: }
924: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
925: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: /*@
930: DMSwarmInitializeVelocitiesFromOptions - Set the initial velocities of particles using a distribution determined from options.
932: Collective
934: Input Parameters:
935: + sw - The `DMSWARM` object
936: - v0 - The velocity scale for nondimensionalization for each species
938: Level: advanced
940: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeCoordinates()`, `DMSwarmInitializeVelocities()`
941: @*/
942: PetscErrorCode DMSwarmInitializeVelocitiesFromOptions(DM sw, const PetscReal v0[])
943: {
944: PetscProbFunc sampler;
945: PetscInt dim;
946: const char *prefix;
947: char funcname[PETSC_MAX_PATH_LEN];
948: PetscBool flg;
950: PetscFunctionBegin;
951: PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
952: PetscCall(PetscOptionsString("-dm_swarm_velocity_function", "Function to determine particle velocities", "DMSwarmSetVelocityFunction", funcname, funcname, sizeof(funcname), &flg));
953: PetscOptionsEnd();
954: if (flg) {
955: PetscSimplePointFn *velFunc;
957: PetscCall(PetscDLSym(NULL, funcname, (void **)&velFunc));
958: PetscCheck(velFunc, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
959: PetscCall(DMSwarmSetVelocityFunction(sw, velFunc));
960: }
961: PetscCall(DMGetDimension(sw, &dim));
962: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
963: PetscCall(PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_velocity_density", NULL, NULL, &sampler));
964: PetscCall(DMSwarmInitializeVelocities(sw, sampler, v0));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }