Actual source code: swarmpic.c
1: #include <petsc/private/dmswarmimpl.h>
2: #include <petscsf.h>
3: #include <petscdmda.h>
4: #include <petscdmplex.h>
5: #include <petscdt.h>
6: #include "../src/dm/impls/swarm/data_bucket.h"
8: #include <petsc/private/petscfeimpl.h>
10: PetscClassId DMSWARMCELLDM_CLASSID;
12: /*@
13: DMSwarmCellDMDestroy - destroy a `DMSwarmCellDM`
15: Collective
17: Input Parameter:
18: . celldm - address of `DMSwarmCellDM`
20: Level: advanced
22: .seealso: `DMSwarmCellDM`, `DMSwarmCellDMCreate()`
23: @*/
24: PetscErrorCode DMSwarmCellDMDestroy(DMSwarmCellDM *celldm)
25: {
26: PetscFunctionBegin;
27: if (!*celldm) PetscFunctionReturn(PETSC_SUCCESS);
29: if (--((PetscObject)*celldm)->refct > 0) {
30: *celldm = NULL;
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
33: PetscTryTypeMethod(*celldm, destroy);
34: for (PetscInt f = 0; f < (*celldm)->Nf; ++f) PetscCall(PetscFree((*celldm)->dmFields[f]));
35: PetscCall(PetscFree((*celldm)->dmFields));
36: for (PetscInt f = 0; f < (*celldm)->Nfc; ++f) PetscCall(PetscFree((*celldm)->coordFields[f]));
37: PetscCall(PetscFree((*celldm)->coordFields));
38: PetscCall(PetscFree((*celldm)->cellid));
39: PetscCall(DMSwarmSortDestroy(&(*celldm)->sort));
40: PetscCall(DMDestroy(&(*celldm)->dm));
41: PetscCall(PetscHeaderDestroy(celldm));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: /*@
46: DMSwarmCellDMView - view a `DMSwarmCellDM`
48: Collective
50: Input Parameters:
51: + celldm - `DMSwarmCellDM`
52: - viewer - viewer to display field, for example `PETSC_VIEWER_STDOUT_WORLD`
54: Level: advanced
56: .seealso: `DMSwarmCellDM`, `DMSwarmCellDMCreate()`
57: @*/
58: PetscErrorCode DMSwarmCellDMView(DMSwarmCellDM celldm, PetscViewer viewer)
59: {
60: PetscBool isascii;
62: PetscFunctionBegin;
64: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)celldm), &viewer));
66: PetscCheckSameComm(celldm, 1, viewer, 2);
67: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
68: if (isascii) {
69: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)celldm, viewer));
70: PetscCall(PetscViewerASCIIPushTab(viewer));
71: PetscCall(PetscViewerASCIIPrintf(viewer, "solution field%s:", celldm->Nf > 1 ? "s" : ""));
72: for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", celldm->dmFields[f]));
73: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
74: PetscCall(PetscViewerASCIIPrintf(viewer, "coordinate field%s:", celldm->Nfc > 1 ? "s" : ""));
75: for (PetscInt f = 0; f < celldm->Nfc; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", celldm->coordFields[f]));
76: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
77: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT));
78: PetscCall(DMView(celldm->dm, viewer));
79: PetscCall(PetscViewerPopFormat(viewer));
80: }
81: PetscTryTypeMethod(celldm, view, viewer);
82: if (isascii) PetscCall(PetscViewerASCIIPopTab(viewer));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: /*@
87: DMSwarmCellDMGetDM - Returns the background `DM` for the `DMSwarm`
89: Not Collective
91: Input Parameter:
92: . celldm - The `DMSwarmCellDM` object
94: Output Parameter:
95: . dm - The `DM` object
97: Level: intermediate
99: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
100: @*/
101: PetscErrorCode DMSwarmCellDMGetDM(DMSwarmCellDM celldm, DM *dm)
102: {
103: PetscFunctionBegin;
105: PetscAssertPointer(dm, 2);
106: *dm = celldm->dm;
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: /*@C
111: DMSwarmCellDMGetFields - Returns the `DM` fields for the `DMSwarm`
113: Not Collective
115: Input Parameter:
116: . celldm - The `DMSwarmCellDM` object
118: Output Parameters:
119: + Nf - The number of fields
120: - names - The array of field names in the `DMSWARM`
122: Level: intermediate
124: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
125: @*/
126: PetscErrorCode DMSwarmCellDMGetFields(DMSwarmCellDM celldm, PetscInt *Nf, const char **names[])
127: {
128: PetscFunctionBegin;
130: if (Nf) {
131: PetscAssertPointer(Nf, 2);
132: *Nf = celldm->Nf;
133: }
134: if (names) {
135: PetscAssertPointer(names, 3);
136: *names = (const char **)celldm->dmFields;
137: }
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*@C
142: DMSwarmCellDMGetCoordinateFields - Returns the `DM` coordinate fields for the `DMSwarm`
144: Not Collective
146: Input Parameter:
147: . celldm - The `DMSwarmCellDM` object
149: Output Parameters:
150: + Nfc - The number of coordinate fields
151: - names - The array of coordinate field names in the `DMSWARM`
153: Level: intermediate
155: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
156: @*/
157: PetscErrorCode DMSwarmCellDMGetCoordinateFields(DMSwarmCellDM celldm, PetscInt *Nfc, const char **names[])
158: {
159: PetscFunctionBegin;
161: if (Nfc) {
162: PetscAssertPointer(Nfc, 2);
163: *Nfc = celldm->Nfc;
164: }
165: if (names) {
166: PetscAssertPointer(names, 3);
167: *names = (const char **)celldm->coordFields;
168: }
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: /*@C
173: DMSwarmCellDMGetCellID - Returns the cell id field name for the `DMSwarm`
175: Not Collective
177: Input Parameter:
178: . celldm - The `DMSwarmCellDM` object
180: Output Parameters:
181: . cellid - The cell id field name in the `DMSWARM`
183: Level: intermediate
185: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
186: @*/
187: PetscErrorCode DMSwarmCellDMGetCellID(DMSwarmCellDM celldm, const char *cellid[])
188: {
189: PetscFunctionBegin;
191: PetscAssertPointer(cellid, 2);
192: *cellid = celldm->cellid;
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: /*@C
197: DMSwarmCellDMGetSort - Returns the sort context over the active `DMSwarmCellDM` for the `DMSwarm`
199: Not Collective
201: Input Parameter:
202: . celldm - The `DMSwarmCellDM` object
204: Output Parameter:
205: . sort - The `DMSwarmSort` object
207: Level: intermediate
209: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`, `DMSwarmSortGetAccess()`, `DMSwarmSortDestroy()`
210: @*/
211: PetscErrorCode DMSwarmCellDMGetSort(DMSwarmCellDM celldm, DMSwarmSort *sort)
212: {
213: PetscFunctionBegin;
215: PetscAssertPointer(sort, 2);
216: *sort = celldm->sort;
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /*@C
221: DMSwarmCellDMSetSort - Sets the sort context over the active `DMSwarmCellDM` for the `DMSwarm`
223: Not Collective
225: Input Parameters:
226: + celldm - The `DMSwarmCellDM` object
227: - sort - The `DMSwarmSort` object, or `NULL` to clear the context
229: Level: intermediate
231: Note:
232: User should destroy `sort` afterward with `DMSwarmSortDestroy()`, as the `DMSwarmCellDM` will hold a reference.
234: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`, `DMSwarmSortGetAccess()`, `DMSwarmSortDestroy()`
235: @*/
236: PetscErrorCode DMSwarmCellDMSetSort(DMSwarmCellDM celldm, DMSwarmSort sort)
237: {
238: PetscFunctionBegin;
240: if (sort) PetscAssertPointer(sort, 2);
241: PetscCall(PetscObjectReference((PetscObject)sort));
242: PetscCall(DMSwarmSortDestroy(&celldm->sort));
243: celldm->sort = sort;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*@
248: DMSwarmCellDMGetBlockSize - Returns the total blocksize for the `DM` fields
250: Not Collective
252: Input Parameters:
253: + celldm - The `DMSwarmCellDM` object
254: - sw - The `DMSwarm` object
256: Output Parameter:
257: . bs - The total block size
259: Level: intermediate
261: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
262: @*/
263: PetscErrorCode DMSwarmCellDMGetBlockSize(DMSwarmCellDM celldm, DM sw, PetscInt *bs)
264: {
265: PetscFunctionBegin;
268: PetscAssertPointer(bs, 3);
269: *bs = 0;
270: for (PetscInt f = 0; f < celldm->Nf; ++f) {
271: PetscInt fbs;
273: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
274: *bs += fbs;
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: DMSwarmCellDMCreate - create a `DMSwarmCellDM`
282: Collective
284: Input Parameters:
285: + dm - The background `DM` for the `DMSwarm`
286: . Nf - The number of swarm fields defined over `dm`
287: . dmFields - The swarm field names for the `dm` fields
288: . Nfc - The number of swarm fields to use for coordinates over `dm`
289: - coordFields - The swarm field names for the `dm` coordinate fields
291: Output Parameter:
292: . celldm - The new `DMSwarmCellDM`
294: Level: advanced
296: .seealso: `DMSwarmCellDM`, `DMSWARM`, `DMSetType()`
297: @*/
298: PetscErrorCode DMSwarmCellDMCreate(DM dm, PetscInt Nf, const char *dmFields[], PetscInt Nfc, const char *coordFields[], DMSwarmCellDM *celldm)
299: {
300: DMSwarmCellDM b;
301: const char *name;
302: char cellid[PETSC_MAX_PATH_LEN];
304: PetscFunctionBegin;
306: if (Nf) PetscAssertPointer(dmFields, 3);
307: if (Nfc) PetscAssertPointer(coordFields, 5);
308: PetscCall(DMInitializePackage());
310: PetscCall(PetscHeaderCreate(b, DMSWARMCELLDM_CLASSID, "DMSwarmCellDM", "Background DM for a Swarm", "DM", PetscObjectComm((PetscObject)dm), DMSwarmCellDMDestroy, DMSwarmCellDMView));
311: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
312: PetscCall(PetscObjectSetName((PetscObject)b, name));
313: PetscCall(PetscObjectReference((PetscObject)dm));
314: b->dm = dm;
315: b->Nf = Nf;
316: b->Nfc = Nfc;
317: PetscCall(PetscMalloc1(b->Nf, &b->dmFields));
318: for (PetscInt f = 0; f < b->Nf; ++f) PetscCall(PetscStrallocpy(dmFields[f], &b->dmFields[f]));
319: PetscCall(PetscMalloc1(b->Nfc, &b->coordFields));
320: for (PetscInt f = 0; f < b->Nfc; ++f) PetscCall(PetscStrallocpy(coordFields[f], &b->coordFields[f]));
321: PetscCall(PetscSNPrintf(cellid, PETSC_MAX_PATH_LEN, "%s_cellid", name));
322: PetscCall(PetscStrallocpy(cellid, &b->cellid));
323: *celldm = b;
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: /* Coordinate insertition/addition API */
328: /*@
329: DMSwarmSetPointsUniformCoordinates - Set point coordinates in a `DMSWARM` on a regular (ijk) grid
331: Collective
333: Input Parameters:
334: + sw - the `DMSWARM`
335: . min - minimum coordinate values in the x, y, z directions (array of length dim)
336: . max - maximum coordinate values in the x, y, z directions (array of length dim)
337: . npoints - number of points in each spatial direction (array of length dim)
338: - mode - indicates whether to append points to the swarm (`ADD_VALUES`), or over-ride existing points (`INSERT_VALUES`)
340: Level: beginner
342: Notes:
343: When using mode = `INSERT_VALUES`, this method will reset the number of particles in the `DMSWARM`
344: to be `npoints[0]` x `npoints[1]` (2D) or `npoints[0]` x `npoints[1]` x `npoints[2]` (3D). When using mode = `ADD_VALUES`,
345: new points will be appended to any already existing in the `DMSWARM`
347: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
348: @*/
349: PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM sw, PetscReal min[], PetscReal max[], PetscInt npoints[], InsertMode mode)
350: {
351: PetscReal lmin[] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
352: PetscReal lmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
353: PetscInt i, j, k, bs, b, n_estimate, n_curr, n_new_est, p, n_found, Nfc;
354: const PetscScalar *_coor;
355: DMSwarmCellDM celldm;
356: DM dm;
357: PetscReal dx[3];
358: PetscInt _npoints[] = {0, 0, 1};
359: Vec pos;
360: PetscScalar *_pos;
361: PetscReal *swarm_coor;
362: PetscInt *swarm_cellid;
363: PetscSF sfcell = NULL;
364: const PetscSFNode *LA_sfcell;
365: const char **coordFields, *cellid;
367: PetscFunctionBegin;
368: DMSWARMPICVALID(sw);
369: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
370: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
371: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
373: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
374: PetscCall(DMGetLocalBoundingBox(dm, lmin, lmax));
375: PetscCall(DMGetCoordinateDim(dm, &bs));
377: for (b = 0; b < bs; b++) {
378: if (npoints[b] > 1) {
379: dx[b] = (max[b] - min[b]) / ((PetscReal)(npoints[b] - 1));
380: } else {
381: dx[b] = 0.0;
382: }
383: _npoints[b] = npoints[b];
384: }
386: /* determine number of points living in the bounding box */
387: n_estimate = 0;
388: for (k = 0; k < _npoints[2]; k++) {
389: for (j = 0; j < _npoints[1]; j++) {
390: for (i = 0; i < _npoints[0]; i++) {
391: PetscReal xp[] = {0.0, 0.0, 0.0};
392: PetscInt ijk[3];
393: PetscBool point_inside = PETSC_TRUE;
395: ijk[0] = i;
396: ijk[1] = j;
397: ijk[2] = k;
398: for (b = 0; b < bs; b++) xp[b] = min[b] + ijk[b] * dx[b];
399: for (b = 0; b < bs; b++) {
400: if (xp[b] < lmin[b]) point_inside = PETSC_FALSE;
401: if (xp[b] > lmax[b]) point_inside = PETSC_FALSE;
402: }
403: if (point_inside) n_estimate++;
404: }
405: }
406: }
408: /* create candidate list */
409: PetscCall(VecCreate(PETSC_COMM_SELF, &pos));
410: PetscCall(VecSetSizes(pos, bs * n_estimate, PETSC_DECIDE));
411: PetscCall(VecSetBlockSize(pos, bs));
412: PetscCall(VecSetFromOptions(pos));
413: PetscCall(VecGetArray(pos, &_pos));
415: n_estimate = 0;
416: for (k = 0; k < _npoints[2]; k++) {
417: for (j = 0; j < _npoints[1]; j++) {
418: for (i = 0; i < _npoints[0]; i++) {
419: PetscReal xp[] = {0.0, 0.0, 0.0};
420: PetscInt ijk[3];
421: PetscBool point_inside = PETSC_TRUE;
423: ijk[0] = i;
424: ijk[1] = j;
425: ijk[2] = k;
426: for (b = 0; b < bs; b++) xp[b] = min[b] + ijk[b] * dx[b];
427: for (b = 0; b < bs; b++) {
428: if (xp[b] < lmin[b]) point_inside = PETSC_FALSE;
429: if (xp[b] > lmax[b]) point_inside = PETSC_FALSE;
430: }
431: if (point_inside) {
432: for (b = 0; b < bs; b++) _pos[bs * n_estimate + b] = xp[b];
433: n_estimate++;
434: }
435: }
436: }
437: }
438: PetscCall(VecRestoreArray(pos, &_pos));
440: /* locate points */
441: PetscCall(DMLocatePoints(dm, pos, DM_POINTLOCATION_NONE, &sfcell));
442: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
443: n_found = 0;
444: for (p = 0; p < n_estimate; p++) {
445: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) n_found++;
446: }
448: /* adjust size */
449: if (mode == ADD_VALUES) {
450: PetscCall(DMSwarmGetLocalSize(sw, &n_curr));
451: n_new_est = n_curr + n_found;
452: PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
453: }
454: if (mode == INSERT_VALUES) {
455: n_curr = 0;
456: n_new_est = n_found;
457: PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
458: }
460: /* initialize new coords, cell owners, pid */
461: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
462: PetscCall(VecGetArrayRead(pos, &_coor));
463: PetscCall(DMSwarmGetField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
464: PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
465: n_found = 0;
466: for (p = 0; p < n_estimate; p++) {
467: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
468: for (b = 0; b < bs; b++) swarm_coor[bs * (n_curr + n_found) + b] = PetscRealPart(_coor[bs * p + b]);
469: swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
470: n_found++;
471: }
472: }
473: PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
474: PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
475: PetscCall(VecRestoreArrayRead(pos, &_coor));
477: PetscCall(PetscSFDestroy(&sfcell));
478: PetscCall(VecDestroy(&pos));
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: /*@
483: DMSwarmSetPointCoordinates - Set point coordinates in a `DMSWARM` from a user defined list
485: Collective
487: Input Parameters:
488: + sw - the `DMSWARM`
489: . npoints - the number of points to insert
490: . coor - the coordinate values
491: . 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
492: - mode - indicates whether to append points to the swarm (`ADD_VALUES`), or over-ride existing points (`INSERT_VALUES`)
494: Level: beginner
496: Notes:
497: If the user has specified `redundant` as `PETSC_FALSE`, the cell `DM` will attempt to locate the coordinates provided by `coor` within
498: its sub-domain. If they any values within `coor` are not located in the sub-domain, they will be ignored and will not get
499: added to the `DMSWARM`.
501: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`, `DMSwarmSetPointsUniformCoordinates()`
502: @*/
503: PetscErrorCode DMSwarmSetPointCoordinates(DM sw, PetscInt npoints, PetscReal coor[], PetscBool redundant, InsertMode mode)
504: {
505: PetscReal gmin[] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
506: PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
507: PetscInt i, N, bs, b, n_estimate, n_curr, n_new_est, p, n_found;
508: Vec coorlocal;
509: const PetscScalar *_coor;
510: DMSwarmCellDM celldm;
511: DM dm;
512: Vec pos;
513: PetscScalar *_pos;
514: PetscReal *swarm_coor;
515: PetscInt *swarm_cellid;
516: PetscSF sfcell = NULL;
517: const PetscSFNode *LA_sfcell;
518: PetscReal *my_coor;
519: PetscInt my_npoints, Nfc;
520: PetscMPIInt rank;
521: MPI_Comm comm;
522: const char **coordFields, *cellid;
524: PetscFunctionBegin;
525: DMSWARMPICVALID(sw);
526: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
527: PetscCallMPI(MPI_Comm_rank(comm, &rank));
529: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
530: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
531: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
533: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
534: PetscCall(DMGetCoordinatesLocal(dm, &coorlocal));
535: PetscCall(VecGetSize(coorlocal, &N));
536: PetscCall(VecGetBlockSize(coorlocal, &bs));
537: N = N / bs;
538: PetscCall(VecGetArrayRead(coorlocal, &_coor));
539: for (i = 0; i < N; i++) {
540: for (b = 0; b < bs; b++) {
541: gmin[b] = PetscMin(gmin[b], PetscRealPart(_coor[bs * i + b]));
542: gmax[b] = PetscMax(gmax[b], PetscRealPart(_coor[bs * i + b]));
543: }
544: }
545: PetscCall(VecRestoreArrayRead(coorlocal, &_coor));
547: /* broadcast points from rank 0 if requested */
548: if (redundant) {
549: PetscMPIInt imy;
551: my_npoints = npoints;
552: PetscCallMPI(MPI_Bcast(&my_npoints, 1, MPIU_INT, 0, comm));
554: if (rank > 0) { /* allocate space */
555: PetscCall(PetscMalloc1(bs * my_npoints, &my_coor));
556: } else {
557: my_coor = coor;
558: }
559: PetscCall(PetscMPIIntCast(bs * my_npoints, &imy));
560: PetscCallMPI(MPI_Bcast(my_coor, imy, MPIU_REAL, 0, comm));
561: } else {
562: my_npoints = npoints;
563: my_coor = coor;
564: }
566: /* determine the number of points living in the bounding box */
567: n_estimate = 0;
568: for (i = 0; i < my_npoints; i++) {
569: PetscBool point_inside = PETSC_TRUE;
571: for (b = 0; b < bs; b++) {
572: if (my_coor[bs * i + b] < gmin[b]) point_inside = PETSC_FALSE;
573: if (my_coor[bs * i + b] > gmax[b]) point_inside = PETSC_FALSE;
574: }
575: if (point_inside) n_estimate++;
576: }
578: /* create candidate list */
579: PetscCall(VecCreate(PETSC_COMM_SELF, &pos));
580: PetscCall(VecSetSizes(pos, bs * n_estimate, PETSC_DECIDE));
581: PetscCall(VecSetBlockSize(pos, bs));
582: PetscCall(VecSetFromOptions(pos));
583: PetscCall(VecGetArray(pos, &_pos));
585: n_estimate = 0;
586: for (i = 0; i < my_npoints; i++) {
587: PetscBool point_inside = PETSC_TRUE;
589: for (b = 0; b < bs; b++) {
590: if (my_coor[bs * i + b] < gmin[b]) point_inside = PETSC_FALSE;
591: if (my_coor[bs * i + b] > gmax[b]) point_inside = PETSC_FALSE;
592: }
593: if (point_inside) {
594: for (b = 0; b < bs; b++) _pos[bs * n_estimate + b] = my_coor[bs * i + b];
595: n_estimate++;
596: }
597: }
598: PetscCall(VecRestoreArray(pos, &_pos));
600: /* locate points */
601: PetscCall(DMLocatePoints(dm, pos, DM_POINTLOCATION_NONE, &sfcell));
603: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
604: n_found = 0;
605: for (p = 0; p < n_estimate; p++) {
606: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) n_found++;
607: }
609: /* adjust size */
610: if (mode == ADD_VALUES) {
611: PetscCall(DMSwarmGetLocalSize(sw, &n_curr));
612: n_new_est = n_curr + n_found;
613: PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
614: }
615: if (mode == INSERT_VALUES) {
616: n_curr = 0;
617: n_new_est = n_found;
618: PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
619: }
621: /* initialize new coords, cell owners, pid */
622: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
623: PetscCall(VecGetArrayRead(pos, &_coor));
624: PetscCall(DMSwarmGetField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
625: PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
626: n_found = 0;
627: for (p = 0; p < n_estimate; p++) {
628: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
629: for (b = 0; b < bs; b++) swarm_coor[bs * (n_curr + n_found) + b] = PetscRealPart(_coor[bs * p + b]);
630: swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
631: n_found++;
632: }
633: }
634: PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
635: PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
636: PetscCall(VecRestoreArrayRead(pos, &_coor));
638: if (redundant) {
639: if (rank > 0) PetscCall(PetscFree(my_coor));
640: }
641: PetscCall(PetscSFDestroy(&sfcell));
642: PetscCall(VecDestroy(&pos));
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM, DM, DMSwarmPICLayoutType, PetscInt);
647: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM, DM, DMSwarmPICLayoutType, PetscInt);
649: /*@
650: DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell
652: Not Collective
654: Input Parameters:
655: + dm - the `DMSWARM`
656: . layout_type - method used to fill each cell with the cell `DM`
657: - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type)
659: Level: beginner
661: Notes:
662: The insert method will reset any previous defined points within the `DMSWARM`.
664: When using a `DMDA` both 2D and 3D are supported for all layout types provided you are using `DMDA_ELEMENT_Q1`.
666: When using a `DMPLEX` the following case are supported\:
667: .vb
668: (i) DMSWARMPIC_LAYOUT_REGULAR: 2D (triangle),
669: (ii) DMSWARMPIC_LAYOUT_GAUSS: 2D and 3D provided the cell is a tri/tet or a quad/hex,
670: (iii) DMSWARMPIC_LAYOUT_SUBDIVISION: 2D and 3D for quad/hex and 2D tri.
671: .ve
673: .seealso: `DMSWARM`, `DMSwarmPICLayoutType`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
674: @*/
675: PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm, DMSwarmPICLayoutType layout_type, PetscInt fill_param)
676: {
677: DM celldm;
678: PetscBool isDA, isPLEX;
680: PetscFunctionBegin;
681: DMSWARMPICVALID(dm);
682: PetscCall(DMSwarmGetCellDM(dm, &celldm));
683: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isDA));
684: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isPLEX));
685: if (isDA) {
686: PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA(dm, celldm, layout_type, fill_param));
687: } else if (isPLEX) {
688: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX(dm, celldm, layout_type, fill_param));
689: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: extern PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *);
695: /*@C
696: DMSwarmSetPointCoordinatesCellwise - Insert point coordinates (defined over the reference cell) within each cell
698: Not Collective
700: Input Parameters:
701: + dm - the `DMSWARM`
702: . npoints - the number of points to insert in each cell
703: - xi - the coordinates (defined in the local coordinate system for each cell) to insert
705: Level: beginner
707: Notes:
708: The method will reset any previous defined points within the `DMSWARM`.
709: Only supported for `DMPLEX`. If you are using a `DMDA` it is recommended to either use
710: `DMSwarmInsertPointsUsingCellDM()`, or extract and set the coordinates yourself the following code
711: .vb
712: PetscReal *coor;
713: const char *coordname;
714: DMSwarmGetCoordinateField(dm, &coordname);
715: DMSwarmGetField(dm,coordname,NULL,NULL,(void**)&coor);
716: // user code to define the coordinates here
717: DMSwarmRestoreField(dm,coordname,NULL,NULL,(void**)&coor);
718: .ve
720: .seealso: `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmInsertPointsUsingCellDM()`
721: @*/
722: PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm, PetscInt npoints, PetscReal xi[])
723: {
724: DM celldm;
725: PetscBool isDA, isPLEX;
727: PetscFunctionBegin;
728: DMSWARMPICVALID(dm);
729: PetscCall(DMSwarmGetCellDM(dm, &celldm));
730: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isDA));
731: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isPLEX));
732: PetscCheck(!isDA, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMPLEX. Recommended you use DMSwarmInsertPointsUsingCellDM()");
733: PetscCheck(isPLEX, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
734: PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, npoints, xi));
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: /*@
739: DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM
741: Not Collective
743: Input Parameter:
744: . sw - the `DMSWARM`
746: Output Parameters:
747: + ncells - the number of cells in the cell `DM` (optional argument, pass `NULL` to ignore)
748: - count - array of length ncells containing the number of points per cell
750: Level: beginner
752: Notes:
753: The array count is allocated internally and must be free'd by the user.
755: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
756: @*/
757: PetscErrorCode DMSwarmCreatePointPerCellCount(DM sw, PetscInt *ncells, PetscInt **count)
758: {
759: DMSwarmCellDM celldm;
760: PetscBool isvalid;
761: PetscInt nel;
762: PetscInt *sum;
763: const char *cellid;
765: PetscFunctionBegin;
766: PetscCall(DMSwarmSortGetIsValid(sw, &isvalid));
767: nel = 0;
768: if (isvalid) {
769: PetscCall(DMSwarmSortGetSizes(sw, &nel, NULL));
771: PetscCall(PetscMalloc1(nel, &sum));
772: for (PetscInt e = 0; e < nel; e++) PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, e, &sum[e]));
773: } else {
774: DM dm;
775: PetscBool isda, isplex, isshell;
776: PetscInt p, npoints;
777: PetscInt *swarm_cellid;
779: /* get the number of cells */
780: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
781: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
782: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
783: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
784: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMSHELL, &isshell));
785: if (isda) {
786: PetscInt _nel, _npe;
787: const PetscInt *_element;
789: PetscCall(DMDAGetElements(dm, &_nel, &_npe, &_element));
790: nel = _nel;
791: PetscCall(DMDARestoreElements(dm, &_nel, &_npe, &_element));
792: } else if (isplex) {
793: PetscInt ps, pe;
795: PetscCall(DMPlexGetHeightStratum(dm, 0, &ps, &pe));
796: nel = pe - ps;
797: } else if (isshell) {
798: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
800: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
801: if (method_DMShellGetNumberOfCells) {
802: PetscCall(method_DMShellGetNumberOfCells(dm, &nel));
803: } else
804: SETERRQ(PetscObjectComm((PetscObject)sw), 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);");
805: } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
807: PetscCall(PetscMalloc1(nel, &sum));
808: PetscCall(PetscArrayzero(sum, nel));
809: PetscCall(DMSwarmGetLocalSize(sw, &npoints));
810: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
811: PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
812: for (p = 0; p < npoints; p++) {
813: if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) sum[swarm_cellid[p]]++;
814: }
815: PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
816: }
817: if (ncells) *ncells = nel;
818: *count = sum;
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: /*@
823: DMSwarmGetNumSpecies - Get the number of particle species
825: Not Collective
827: Input Parameter:
828: . sw - the `DMSWARM`
830: Output Parameters:
831: . Ns - the number of species
833: Level: intermediate
835: .seealso: `DMSWARM`, `DMSwarmSetNumSpecies()`, `DMSwarmSetType()`, `DMSwarmType`
836: @*/
837: PetscErrorCode DMSwarmGetNumSpecies(DM sw, PetscInt *Ns)
838: {
839: DM_Swarm *swarm = (DM_Swarm *)sw->data;
841: PetscFunctionBegin;
842: *Ns = swarm->Ns;
843: PetscFunctionReturn(PETSC_SUCCESS);
844: }
846: /*@
847: DMSwarmSetNumSpecies - Set the number of particle species
849: Not Collective
851: Input Parameters:
852: + sw - the `DMSWARM`
853: - Ns - the number of species
855: Level: intermediate
857: .seealso: `DMSWARM`, `DMSwarmGetNumSpecies()`, `DMSwarmSetType()`, `DMSwarmType`
858: @*/
859: PetscErrorCode DMSwarmSetNumSpecies(DM sw, PetscInt Ns)
860: {
861: DM_Swarm *swarm = (DM_Swarm *)sw->data;
863: PetscFunctionBegin;
864: swarm->Ns = Ns;
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*@C
869: DMSwarmGetCoordinateFunction - Get the function setting initial particle positions, if it exists
871: Not Collective
873: Input Parameter:
874: . sw - the `DMSWARM`
876: Output Parameter:
877: . coordFunc - the function setting initial particle positions, or `NULL`, see `PetscSimplePointFn` for the calling sequence
879: Level: intermediate
881: .seealso: `DMSWARM`, `DMSwarmSetCoordinateFunction()`, `DMSwarmGetVelocityFunction()`, `DMSwarmInitializeCoordinates()`, `PetscSimplePointFn`
882: @*/
883: PetscErrorCode DMSwarmGetCoordinateFunction(DM sw, PetscSimplePointFn **coordFunc)
884: {
885: DM_Swarm *swarm = (DM_Swarm *)sw->data;
887: PetscFunctionBegin;
889: PetscAssertPointer(coordFunc, 2);
890: *coordFunc = swarm->coordFunc;
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: /*@C
895: DMSwarmSetCoordinateFunction - Set the function setting initial particle positions
897: Not Collective
899: Input Parameters:
900: + sw - the `DMSWARM`
901: - coordFunc - the function setting initial particle positions, see `PetscSimplePointFn` for the calling sequence
903: Level: intermediate
905: .seealso: `DMSWARM`, `DMSwarmGetCoordinateFunction()`, `DMSwarmSetVelocityFunction()`, `DMSwarmInitializeCoordinates()`, `PetscSimplePointFn`
906: @*/
907: PetscErrorCode DMSwarmSetCoordinateFunction(DM sw, PetscSimplePointFn *coordFunc)
908: {
909: DM_Swarm *swarm = (DM_Swarm *)sw->data;
911: PetscFunctionBegin;
914: swarm->coordFunc = coordFunc;
915: PetscFunctionReturn(PETSC_SUCCESS);
916: }
918: /*@C
919: DMSwarmGetVelocityFunction - Get the function setting initial particle velocities, if it exists
921: Not Collective
923: Input Parameter:
924: . sw - the `DMSWARM`
926: Output Parameter:
927: . velFunc - the function setting initial particle velocities, or `NULL`, see `PetscSimplePointFn` for the calling sequence
929: Level: intermediate
931: .seealso: `DMSWARM`, `DMSwarmSetVelocityFunction()`, `DMSwarmGetCoordinateFunction()`, `DMSwarmInitializeVelocities()`, `PetscSimplePointFn`
932: @*/
933: PetscErrorCode DMSwarmGetVelocityFunction(DM sw, PetscSimplePointFn **velFunc)
934: {
935: DM_Swarm *swarm = (DM_Swarm *)sw->data;
937: PetscFunctionBegin;
939: PetscAssertPointer(velFunc, 2);
940: *velFunc = swarm->velFunc;
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: /*@C
945: DMSwarmSetVelocityFunction - Set the function setting initial particle velocities
947: Not Collective
949: Input Parameters:
950: + sw - the `DMSWARM`
951: - velFunc - the function setting initial particle velocities, see `PetscSimplePointFn` for the calling sequence
953: Level: intermediate
955: .seealso: `DMSWARM`, `DMSwarmGetVelocityFunction()`, `DMSwarmSetCoordinateFunction()`, `DMSwarmInitializeVelocities()`, `PetscSimplePointFn`
956: @*/
957: PetscErrorCode DMSwarmSetVelocityFunction(DM sw, PetscSimplePointFn *velFunc)
958: {
959: DM_Swarm *swarm = (DM_Swarm *)sw->data;
961: PetscFunctionBegin;
964: swarm->velFunc = velFunc;
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: /*@C
969: DMSwarmComputeLocalSize - Compute the local number and distribution of particles based upon a density function
971: Not Collective
973: Input Parameters:
974: + sw - The `DMSWARM`
975: . N - The target number of particles
976: - density - The density field for the particle layout, normalized to unity
978: Level: advanced
980: Note:
981: One particle will be created for each species.
983: .seealso: `DMSWARM`, `DMSwarmComputeLocalSizeFromOptions()`
984: @*/
985: PetscErrorCode DMSwarmComputeLocalSize(DM sw, PetscInt N, PetscProbFn *density)
986: {
987: DM dm;
988: DMSwarmCellDM celldm;
989: PetscQuadrature quad;
990: const PetscReal *xq, *wq;
991: PetscReal *n_int;
992: PetscInt *npc_s, *swarm_cellid, Ni;
993: PetscReal gmin[3], gmax[3], xi0[3];
994: PetscInt Ns, cStart, cEnd, c, dim, d, Nq, q, Np = 0, p, s;
995: PetscBool simplex;
996: const char *cellid;
998: PetscFunctionBegin;
999: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1000: PetscCall(DMSwarmGetCellDM(sw, &dm));
1001: PetscCall(DMGetDimension(dm, &dim));
1002: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
1003: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1004: PetscCall(DMPlexIsSimplex(dm, &simplex));
1005: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1006: if (simplex) PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 5, -1.0, 1.0, &quad));
1007: else PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 5, -1.0, 1.0, &quad));
1008: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xq, &wq));
1009: PetscCall(PetscCalloc2(Ns, &n_int, (cEnd - cStart) * Ns, &npc_s));
1010: /* Integrate the density function to get the number of particles in each cell */
1011: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1012: for (c = 0; c < cEnd - cStart; ++c) {
1013: const PetscInt cell = c + cStart;
1014: PetscReal v0[3], J[9], invJ[9], detJ, detJp = 2. / (gmax[0] - gmin[0]), xr[3], den;
1016: /* Have to transform quadrature points/weights to cell domain */
1017: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
1018: PetscCall(PetscArrayzero(n_int, Ns));
1019: for (q = 0; q < Nq; ++q) {
1020: CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], xr);
1021: /* Have to transform mesh to domain of definition of PDF, [-1, 1], and weight PDF by |J|/2 */
1022: xr[0] = detJp * (xr[0] - gmin[0]) - 1.;
1024: for (s = 0; s < Ns; ++s) {
1025: PetscCall(density(xr, NULL, &den));
1026: n_int[s] += (detJp * den) * (detJ * wq[q]) / (PetscReal)Ns;
1027: }
1028: }
1029: for (s = 0; s < Ns; ++s) {
1030: Ni = N;
1031: npc_s[c * Ns + s] += (PetscInt)(Ni * n_int[s] + 0.5); // TODO Wish we wrapped round()
1032: Np += npc_s[c * Ns + s];
1033: }
1034: }
1035: PetscCall(PetscQuadratureDestroy(&quad));
1036: PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
1037: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1038: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
1039: PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
1040: for (c = 0, p = 0; c < cEnd - cStart; ++c) {
1041: for (s = 0; s < Ns; ++s) {
1042: for (q = 0; q < npc_s[c * Ns + s]; ++q, ++p) swarm_cellid[p] = c;
1043: }
1044: }
1045: PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
1046: PetscCall(PetscFree2(n_int, npc_s));
1047: PetscFunctionReturn(PETSC_SUCCESS);
1048: }
1050: /*@
1051: DMSwarmComputeLocalSizeFromOptions - Compute the local number and distribution of particles based upon a density function determined by options
1053: Not Collective
1055: Input Parameter:
1056: . sw - The `DMSWARM`
1058: Level: advanced
1060: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`
1061: @*/
1062: PetscErrorCode DMSwarmComputeLocalSizeFromOptions(DM sw)
1063: {
1064: PetscProbFn *pdf;
1065: const char *prefix;
1066: char funcname[PETSC_MAX_PATH_LEN];
1067: PetscInt *N, Ns, dim, n;
1068: PetscBool flg;
1069: PetscMPIInt size, rank;
1071: PetscFunctionBegin;
1072: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sw), &size));
1073: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1074: PetscCall(PetscCalloc1(size, &N));
1075: PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
1076: n = size;
1077: PetscCall(PetscOptionsIntArray("-dm_swarm_num_particles", "The target number of particles", "", N, &n, NULL));
1078: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1079: PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
1080: if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
1081: PetscCall(PetscOptionsString("-dm_swarm_coordinate_function", "Function to determine particle coordinates", "DMSwarmSetCoordinateFunction", funcname, funcname, sizeof(funcname), &flg));
1082: PetscOptionsEnd();
1083: if (flg) {
1084: PetscSimplePointFn *coordFunc;
1086: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1087: PetscCall(PetscDLSym(NULL, funcname, (void **)&coordFunc));
1088: PetscCheck(coordFunc, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
1089: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1090: PetscCall(DMSwarmSetLocalSizes(sw, N[rank] * Ns, 0));
1091: PetscCall(DMSwarmSetCoordinateFunction(sw, coordFunc));
1092: } else {
1093: PetscCall(DMGetDimension(sw, &dim));
1094: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
1095: PetscCall(PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_coordinate_density", &pdf, NULL, NULL));
1096: PetscCall(DMSwarmComputeLocalSize(sw, N[rank], pdf));
1097: }
1098: PetscCall(PetscFree(N));
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: /*@
1103: DMSwarmInitializeCoordinates - Determine the initial coordinates of particles for a PIC method
1105: Not Collective
1107: Input Parameter:
1108: . sw - The `DMSWARM`
1110: Level: advanced
1112: Note:
1113: Currently, we randomly place particles in their assigned cell
1115: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeVelocities()`
1116: @*/
1117: PetscErrorCode DMSwarmInitializeCoordinates(DM sw)
1118: {
1119: DMSwarmCellDM celldm;
1120: PetscSimplePointFn *coordFunc;
1121: PetscScalar *weight;
1122: PetscReal *x;
1123: PetscInt *species;
1124: void *ctx;
1125: PetscBool removePoints = PETSC_TRUE;
1126: PetscDataType dtype;
1127: PetscInt Nfc, Np, p, Ns, dim, d, bs;
1128: const char **coordFields;
1130: PetscFunctionBeginUser;
1131: PetscCall(DMGetDimension(sw, &dim));
1132: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1133: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1134: PetscCall(DMSwarmGetCoordinateFunction(sw, &coordFunc));
1136: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1137: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1138: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
1140: PetscCall(DMSwarmGetField(sw, coordFields[0], &bs, &dtype, (void **)&x));
1141: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&weight));
1142: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1143: if (coordFunc) {
1144: PetscCall(DMGetApplicationContext(sw, &ctx));
1145: for (p = 0; p < Np; ++p) {
1146: PetscScalar X[3];
1148: PetscCall((*coordFunc)(dim, 0., NULL, p, X, ctx));
1149: for (d = 0; d < dim; ++d) x[p * dim + d] = PetscRealPart(X[d]);
1150: weight[p] = 1.0;
1151: species[p] = p % Ns;
1152: }
1153: } else {
1154: DM dm;
1155: PetscRandom rnd;
1156: PetscReal xi0[3];
1157: PetscInt cStart, cEnd, c;
1159: PetscCall(DMSwarmGetCellDM(sw, &dm));
1160: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1161: PetscCall(DMGetApplicationContext(sw, &ctx));
1163: /* Set particle position randomly in cell, set weights to 1 */
1164: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
1165: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
1166: PetscCall(PetscRandomSetFromOptions(rnd));
1167: PetscCall(DMSwarmSortGetAccess(sw));
1168: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1169: for (c = cStart; c < cEnd; ++c) {
1170: PetscReal v0[3], J[9], invJ[9], detJ;
1171: PetscInt *pidx, Npc, q;
1173: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1174: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
1175: for (q = 0; q < Npc; ++q) {
1176: const PetscInt p = pidx[q];
1177: PetscReal xref[3];
1179: for (d = 0; d < dim; ++d) PetscCall(PetscRandomGetValueReal(rnd, &xref[d]));
1180: CoordinatesRefToReal(dim, dim, xi0, v0, J, xref, &x[p * dim]);
1182: weight[p] = 1.0 / Np;
1183: species[p] = p % Ns;
1184: }
1185: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1186: }
1187: PetscCall(PetscRandomDestroy(&rnd));
1188: PetscCall(DMSwarmSortRestoreAccess(sw));
1189: }
1190: PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&x));
1191: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1192: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1194: PetscCall(DMSwarmMigrate(sw, removePoints));
1195: PetscCall(DMLocalizeCoordinates(sw));
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: /*@C
1200: DMSwarmInitializeVelocities - Set the initial velocities of particles using a distribution.
1202: Collective
1204: Input Parameters:
1205: + sw - The `DMSWARM` object
1206: . sampler - A function which uniformly samples the velocity PDF
1207: - v0 - The velocity scale for nondimensionalization for each species
1209: Level: advanced
1211: Note:
1212: 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.
1214: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeCoordinates()`, `DMSwarmInitializeVelocitiesFromOptions()`
1215: @*/
1216: PetscErrorCode DMSwarmInitializeVelocities(DM sw, PetscProbFn *sampler, const PetscReal v0[])
1217: {
1218: PetscSimplePointFn *velFunc;
1219: PetscReal *v;
1220: PetscInt *species;
1221: void *ctx;
1222: PetscInt dim, Np, p;
1224: PetscFunctionBegin;
1225: PetscCall(DMSwarmGetVelocityFunction(sw, &velFunc));
1227: PetscCall(DMGetDimension(sw, &dim));
1228: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1229: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
1230: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1231: if (v0[0] == 0.) {
1232: PetscCall(PetscArrayzero(v, Np * dim));
1233: } else if (velFunc) {
1234: PetscCall(DMGetApplicationContext(sw, &ctx));
1235: for (p = 0; p < Np; ++p) {
1236: PetscInt s = species[p], d;
1237: PetscScalar vel[3];
1239: PetscCall((*velFunc)(dim, 0., NULL, p, vel, ctx));
1240: for (d = 0; d < dim; ++d) v[p * dim + d] = (v0[s] / v0[0]) * PetscRealPart(vel[d]);
1241: }
1242: } else {
1243: PetscRandom rnd;
1245: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd));
1246: PetscCall(PetscRandomSetInterval(rnd, 0, 1.));
1247: PetscCall(PetscRandomSetFromOptions(rnd));
1249: for (p = 0; p < Np; ++p) {
1250: PetscInt s = species[p], d;
1251: PetscReal a[3], vel[3];
1253: for (d = 0; d < dim; ++d) PetscCall(PetscRandomGetValueReal(rnd, &a[d]));
1254: PetscCall(sampler(a, NULL, vel));
1255: for (d = 0; d < dim; ++d) v[p * dim + d] = (v0[s] / v0[0]) * vel[d];
1256: }
1257: PetscCall(PetscRandomDestroy(&rnd));
1258: }
1259: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1260: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1261: PetscFunctionReturn(PETSC_SUCCESS);
1262: }
1264: /*@
1265: DMSwarmInitializeVelocitiesFromOptions - Set the initial velocities of particles using a distribution determined from options.
1267: Collective
1269: Input Parameters:
1270: + sw - The `DMSWARM` object
1271: - v0 - The velocity scale for nondimensionalization for each species
1273: Level: advanced
1275: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeCoordinates()`, `DMSwarmInitializeVelocities()`
1276: @*/
1277: PetscErrorCode DMSwarmInitializeVelocitiesFromOptions(DM sw, const PetscReal v0[])
1278: {
1279: PetscProbFn *sampler;
1280: PetscInt dim;
1281: const char *prefix;
1282: char funcname[PETSC_MAX_PATH_LEN];
1283: PetscBool flg;
1285: PetscFunctionBegin;
1286: PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
1287: PetscCall(PetscOptionsString("-dm_swarm_velocity_function", "Function to determine particle velocities", "DMSwarmSetVelocityFunction", funcname, funcname, sizeof(funcname), &flg));
1288: PetscOptionsEnd();
1289: if (flg) {
1290: PetscSimplePointFn *velFunc;
1292: PetscCall(PetscDLSym(NULL, funcname, (void **)&velFunc));
1293: PetscCheck(velFunc, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
1294: PetscCall(DMSwarmSetVelocityFunction(sw, velFunc));
1295: }
1296: PetscCall(DMGetDimension(sw, &dim));
1297: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
1298: PetscCall(PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_velocity_density", NULL, NULL, &sampler));
1299: PetscCall(DMSwarmInitializeVelocities(sw, sampler, v0));
1300: PetscFunctionReturn(PETSC_SUCCESS);
1301: }
1303: // The input vector U is assumed to be from a PetscFE. The Swarm fields are input as auxiliary values.
1304: PetscErrorCode DMProjectFieldLocal_Swarm(DM dm, PetscReal time, Vec U, PetscPointFn **funcs, InsertMode mode, Vec X)
1305: {
1306: MPI_Comm comm;
1307: DM dmIn;
1308: PetscDS ds;
1309: PetscTabulation *T;
1310: DMSwarmCellDM celldm;
1311: PetscScalar *a, *val, *u, *u_x;
1312: PetscFEGeom fegeom;
1313: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
1314: PetscInt dim, dE, Np, n, Nf, Nfc, Nfu, cStart, cEnd, maxC = 0, totbs = 0;
1315: const char **coordFields, **fields;
1316: PetscReal **coordVals, **vals;
1317: PetscInt *cbs, *bs, *uOff, *uOff_x;
1319: PetscFunctionBegin;
1320: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1321: PetscCall(VecGetDM(U, &dmIn));
1322: PetscCall(DMGetDimension(dmIn, &dim));
1323: PetscCall(DMGetCoordinateDim(dmIn, &dE));
1324: PetscCall(DMGetDS(dmIn, &ds));
1325: PetscCall(PetscDSGetNumFields(ds, &Nfu));
1326: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
1327: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
1328: PetscCall(PetscDSGetTabulation(ds, &T));
1329: PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x));
1330: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
1331: PetscCall(DMPlexGetHeightStratum(dmIn, 0, &cStart, &cEnd));
1333: fegeom.dim = dim;
1334: fegeom.dimEmbed = dE;
1335: fegeom.v = v0;
1336: fegeom.xi = v0ref;
1337: fegeom.J = J;
1338: fegeom.invJ = invJ;
1339: fegeom.detJ = &detJ;
1341: PetscCall(DMSwarmGetLocalSize(dm, &Np));
1342: PetscCall(VecGetLocalSize(X, &n));
1343: PetscCheck(n == Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Output vector local size %" PetscInt_FMT " != %" PetscInt_FMT " number of local particles", n, Np);
1344: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
1345: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1346: PetscCall(DMSwarmCellDMGetFields(celldm, &Nf, &fields));
1348: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &cbs));
1349: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dm, coordFields[i], &cbs[i], NULL, (void **)&coordVals[i]));
1350: PetscCall(PetscMalloc2(Nf, &vals, Nfc, &bs));
1351: for (PetscInt i = 0; i < Nf; ++i) {
1352: PetscCall(DMSwarmGetField(dm, fields[i], &bs[i], NULL, (void **)&vals[i]));
1353: totbs += bs[i];
1354: }
1356: PetscCall(DMSwarmSortGetAccess(dm));
1357: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1358: PetscInt *pindices, Npc;
1360: PetscCall(DMSwarmSortGetPointsPerCell(dm, cell, &Npc, &pindices));
1361: maxC = PetscMax(maxC, Npc);
1362: PetscCall(DMSwarmSortRestorePointsPerCell(dm, cell, &Npc, &pindices));
1363: }
1364: PetscCall(PetscMalloc3(maxC * dim, &xi, maxC * totbs, &val, Nfu, &T));
1365: PetscCall(VecGetArray(X, &a));
1366: {
1367: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1368: PetscInt *pindices, Npc;
1370: // TODO: Use DMField instead of assuming affine
1371: PetscCall(DMPlexComputeCellGeometryFEM(dmIn, cell, NULL, v0, J, invJ, &detJ));
1372: PetscCall(DMSwarmSortGetPointsPerCell(dm, cell, &Npc, &pindices));
1374: PetscScalar *closure = NULL;
1375: PetscInt Ncl;
1377: // Get fields from input vector and auxiliary fields from swarm
1378: for (PetscInt p = 0; p < Npc; ++p) {
1379: PetscReal xr[8];
1380: PetscInt off;
1382: off = 0;
1383: for (PetscInt i = 0; i < Nfc; ++i) {
1384: for (PetscInt b = 0; b < cbs[i]; ++b, ++off) xr[off] = coordVals[i][pindices[p] * cbs[i] + b];
1385: }
1386: PetscCheck(off == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, dim);
1387: CoordinatesRealToRef(dE, dim, fegeom.xi, fegeom.v, fegeom.invJ, xr, &xi[p * dim]);
1388: off = 0;
1389: for (PetscInt i = 0; i < Nf; ++i) {
1390: for (PetscInt b = 0; b < bs[i]; ++b, ++off) val[p * totbs + off] = vals[i][pindices[p] * bs[i] + b];
1391: }
1392: PetscCheck(off == totbs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of swarm fields is %" PetscInt_FMT " != %" PetscInt_FMT " the computed total block size", off, totbs);
1393: }
1394: PetscCall(DMPlexVecGetClosure(dmIn, NULL, U, cell, &Ncl, &closure));
1395: for (PetscInt field = 0; field < Nfu; ++field) {
1396: PetscFE fe;
1398: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1399: PetscCall(PetscFECreateTabulation(fe, 1, Npc, xi, 1, &T[field]));
1400: }
1401: for (PetscInt p = 0; p < Npc; ++p) {
1402: // Get fields from input vector
1403: PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nfu, 0, p, T, &fegeom, closure, NULL, u, u_x, NULL));
1404: (*funcs[0])(dim, 1, 1, uOff, uOff_x, u, NULL, u_x, bs, NULL, &val[p * totbs], NULL, NULL, time, &xi[p * dim], 0, NULL, &a[pindices[p]]);
1405: }
1406: PetscCall(DMSwarmSortRestorePointsPerCell(dm, cell, &Npc, &pindices));
1407: PetscCall(DMPlexVecRestoreClosure(dmIn, NULL, U, cell, &Ncl, &closure));
1408: for (PetscInt field = 0; field < Nfu; ++field) PetscCall(PetscTabulationDestroy(&T[field]));
1409: }
1410: }
1411: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dm, coordFields[i], &cbs[i], NULL, (void **)&coordVals[i]));
1412: for (PetscInt i = 0; i < Nf; ++i) PetscCall(DMSwarmRestoreField(dm, fields[i], &bs[i], NULL, (void **)&vals[i]));
1413: PetscCall(VecRestoreArray(X, &a));
1414: PetscCall(DMSwarmSortRestoreAccess(dm));
1415: PetscCall(PetscFree3(xi, val, T));
1416: PetscCall(PetscFree3(v0, J, invJ));
1417: PetscCall(PetscFree2(coordVals, cbs));
1418: PetscCall(PetscFree2(vals, bs));
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }