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: PetscClassId DMSWARMCELLDM_CLASSID;
13: /*@
14: DMSwarmCellDMDestroy - destroy a `DMSwarmCellDM`
16: Collective
18: Input Parameter:
19: . celldm - address of `DMSwarmCellDM`
21: Level: advanced
23: .seealso: `DMSwarmCellDM`, `DMSwarmCellDMCreate()`
24: @*/
25: PetscErrorCode DMSwarmCellDMDestroy(DMSwarmCellDM *celldm)
26: {
27: PetscFunctionBegin;
28: if (!*celldm) PetscFunctionReturn(PETSC_SUCCESS);
30: if (--((PetscObject)*celldm)->refct > 0) {
31: *celldm = NULL;
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
34: PetscTryTypeMethod(*celldm, destroy);
35: for (PetscInt f = 0; f < (*celldm)->Nf; ++f) PetscCall(PetscFree((*celldm)->dmFields[f]));
36: PetscCall(PetscFree((*celldm)->dmFields));
37: for (PetscInt f = 0; f < (*celldm)->Nfc; ++f) PetscCall(PetscFree((*celldm)->coordFields[f]));
38: PetscCall(PetscFree((*celldm)->coordFields));
39: PetscCall(PetscFree((*celldm)->cellid));
40: PetscCall(DMSwarmSortDestroy(&(*celldm)->sort));
41: PetscCall(DMDestroy(&(*celldm)->dm));
42: PetscCall(PetscHeaderDestroy(celldm));
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*@
47: DMSwarmCellDMView - view a `DMSwarmCellDM`
49: Collective
51: Input Parameters:
52: + celldm - `DMSwarmCellDM`
53: - viewer - viewer to display field, for example `PETSC_VIEWER_STDOUT_WORLD`
55: Level: advanced
57: .seealso: `DMSwarmCellDM`, `DMSwarmCellDMCreate()`
58: @*/
59: PetscErrorCode DMSwarmCellDMView(DMSwarmCellDM celldm, PetscViewer viewer)
60: {
61: PetscBool iascii;
63: PetscFunctionBegin;
65: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)celldm), &viewer));
67: PetscCheckSameComm(celldm, 1, viewer, 2);
68: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
69: if (iascii) {
70: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)celldm, viewer));
71: PetscCall(PetscViewerASCIIPushTab(viewer));
72: PetscCall(PetscViewerASCIIPrintf(viewer, "solution field%s:", celldm->Nf > 1 ? "s" : ""));
73: for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", celldm->dmFields[f]));
74: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
75: PetscCall(PetscViewerASCIIPrintf(viewer, "coordinate field%s:", celldm->Nfc > 1 ? "s" : ""));
76: for (PetscInt f = 0; f < celldm->Nfc; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", celldm->coordFields[f]));
77: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
78: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT));
79: PetscCall(DMView(celldm->dm, viewer));
80: PetscCall(PetscViewerPopFormat(viewer));
81: }
82: PetscTryTypeMethod(celldm, view, viewer);
83: if (iascii) PetscCall(PetscViewerASCIIPopTab(viewer));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: /*@
88: DMSwarmCellDMGetDM - Returns the background `DM` for the `DMSwarm`
90: Not Collective
92: Input Parameter:
93: . celldm - The `DMSwarmCellDM` object
95: Output Parameter:
96: . dm - The `DM` object
98: Level: intermediate
100: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
101: @*/
102: PetscErrorCode DMSwarmCellDMGetDM(DMSwarmCellDM celldm, DM *dm)
103: {
104: PetscFunctionBegin;
106: PetscAssertPointer(dm, 2);
107: *dm = celldm->dm;
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: /*@C
112: DMSwarmCellDMGetFields - Returns the `DM` fields for the `DMSwarm`
114: Not Collective
116: Input Parameter:
117: . celldm - The `DMSwarmCellDM` object
119: Output Parameters:
120: + Nf - The number of fields
121: - names - The array of field names in the `DMSWARM`
123: Level: intermediate
125: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
126: @*/
127: PetscErrorCode DMSwarmCellDMGetFields(DMSwarmCellDM celldm, PetscInt *Nf, const char **names[])
128: {
129: PetscFunctionBegin;
131: if (Nf) {
132: PetscAssertPointer(Nf, 2);
133: *Nf = celldm->Nf;
134: }
135: if (names) {
136: PetscAssertPointer(names, 3);
137: *names = (const char **)celldm->dmFields;
138: }
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /*@C
143: DMSwarmCellDMGetCoordinateFields - Returns the `DM` coordinate fields for the `DMSwarm`
145: Not Collective
147: Input Parameter:
148: . celldm - The `DMSwarmCellDM` object
150: Output Parameters:
151: + Nfc - The number of coordinate fields
152: - names - The array of coordinate field names in the `DMSWARM`
154: Level: intermediate
156: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
157: @*/
158: PetscErrorCode DMSwarmCellDMGetCoordinateFields(DMSwarmCellDM celldm, PetscInt *Nfc, const char **names[])
159: {
160: PetscFunctionBegin;
162: if (Nfc) {
163: PetscAssertPointer(Nfc, 2);
164: *Nfc = celldm->Nfc;
165: }
166: if (names) {
167: PetscAssertPointer(names, 3);
168: *names = (const char **)celldm->coordFields;
169: }
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /*@C
174: DMSwarmCellDMGetCellID - Returns the cell id field name for the `DMSwarm`
176: Not Collective
178: Input Parameter:
179: . celldm - The `DMSwarmCellDM` object
181: Output Parameters:
182: . cellid - The cell id field name in the `DMSWARM`
184: Level: intermediate
186: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
187: @*/
188: PetscErrorCode DMSwarmCellDMGetCellID(DMSwarmCellDM celldm, const char *cellid[])
189: {
190: PetscFunctionBegin;
192: PetscAssertPointer(cellid, 2);
193: *cellid = celldm->cellid;
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*@C
198: DMSwarmCellDMGetSort - Returns the sort context over the active `DMSwarmCellDM` for the `DMSwarm`
200: Not Collective
202: Input Parameter:
203: . celldm - The `DMSwarmCellDM` object
205: Output Parameter:
206: . sort - The `DMSwarmSort` object
208: Level: intermediate
210: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
211: @*/
212: PetscErrorCode DMSwarmCellDMGetSort(DMSwarmCellDM celldm, DMSwarmSort *sort)
213: {
214: PetscFunctionBegin;
216: PetscAssertPointer(sort, 2);
217: *sort = celldm->sort;
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: /*@C
222: DMSwarmCellDMSetSort - Sets the sort context over the active `DMSwarmCellDM` for the `DMSwarm`
224: Not Collective
226: Input Parameters:
227: + celldm - The `DMSwarmCellDM` object
228: - sort - The `DMSwarmSort` object
230: Level: intermediate
232: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
233: @*/
234: PetscErrorCode DMSwarmCellDMSetSort(DMSwarmCellDM celldm, DMSwarmSort sort)
235: {
236: PetscFunctionBegin;
238: PetscAssertPointer(sort, 2);
239: celldm->sort = sort;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@
244: DMSwarmCellDMGetBlockSize - Returns the total blocksize for the `DM` fields
246: Not Collective
248: Input Parameters:
249: + celldm - The `DMSwarmCellDM` object
250: - sw - The `DMSwarm` object
252: Output Parameter:
253: . bs - The total block size
255: Level: intermediate
257: .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
258: @*/
259: PetscErrorCode DMSwarmCellDMGetBlockSize(DMSwarmCellDM celldm, DM sw, PetscInt *bs)
260: {
261: PetscFunctionBegin;
264: PetscAssertPointer(bs, 3);
265: *bs = 0;
266: for (PetscInt f = 0; f < celldm->Nf; ++f) {
267: PetscInt fbs;
269: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
270: *bs += fbs;
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*@
276: DMSwarmCellDMCreate - create a `DMSwarmCellDM`
278: Collective
280: Input Parameter:
281: + dm - The background `DM` for the `DMSwarm`
282: . Nf - The number of swarm fields defined over `dm`
283: . dmFields - The swarm field names for the `dm` fields
284: . Nfc - The number of swarm fields to use for coordinates over `dm`
285: - coordFields - The swarm field names for the `dm` coordinate fields
287: Output Parameter:
288: . celldm - The new `DMSwarmCellDM`
290: Level: advanced
292: .seealso: `DMSwarmCellDM`, `DMSWARM`, `DMSetType()`
293: @*/
294: PetscErrorCode DMSwarmCellDMCreate(DM dm, PetscInt Nf, const char *dmFields[], PetscInt Nfc, const char *coordFields[], DMSwarmCellDM *celldm)
295: {
296: DMSwarmCellDM b;
297: const char *name;
298: char cellid[PETSC_MAX_PATH_LEN];
300: PetscFunctionBegin;
302: if (Nf) PetscAssertPointer(dmFields, 3);
303: if (Nfc) PetscAssertPointer(coordFields, 5);
304: PetscCall(DMInitializePackage());
306: PetscCall(PetscHeaderCreate(b, DMSWARMCELLDM_CLASSID, "DMSwarmCellDM", "Background DM for a Swarm", "DM", PetscObjectComm((PetscObject)dm), DMSwarmCellDMDestroy, DMSwarmCellDMView));
307: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
308: PetscCall(PetscObjectSetName((PetscObject)b, name));
309: PetscCall(PetscObjectReference((PetscObject)dm));
310: b->dm = dm;
311: b->Nf = Nf;
312: b->Nfc = Nfc;
313: PetscCall(PetscMalloc1(b->Nf, &b->dmFields));
314: for (PetscInt f = 0; f < b->Nf; ++f) PetscCall(PetscStrallocpy(dmFields[f], &b->dmFields[f]));
315: PetscCall(PetscMalloc1(b->Nfc, &b->coordFields));
316: for (PetscInt f = 0; f < b->Nfc; ++f) PetscCall(PetscStrallocpy(coordFields[f], &b->coordFields[f]));
317: PetscCall(PetscSNPrintf(cellid, PETSC_MAX_PATH_LEN, "%s_cellid", name));
318: PetscCall(PetscStrallocpy(cellid, &b->cellid));
319: *celldm = b;
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /* Coordinate insertition/addition API */
324: /*@
325: DMSwarmSetPointsUniformCoordinates - Set point coordinates in a `DMSWARM` on a regular (ijk) grid
327: Collective
329: Input Parameters:
330: + sw - the `DMSWARM`
331: . min - minimum coordinate values in the x, y, z directions (array of length dim)
332: . max - maximum coordinate values in the x, y, z directions (array of length dim)
333: . npoints - number of points in each spatial direction (array of length dim)
334: - mode - indicates whether to append points to the swarm (`ADD_VALUES`), or over-ride existing points (`INSERT_VALUES`)
336: Level: beginner
338: Notes:
339: When using mode = `INSERT_VALUES`, this method will reset the number of particles in the `DMSWARM`
340: to be `npoints[0]` x `npoints[1]` (2D) or `npoints[0]` x `npoints[1]` x `npoints[2]` (3D). When using mode = `ADD_VALUES`,
341: new points will be appended to any already existing in the `DMSWARM`
343: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
344: @*/
345: PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM sw, PetscReal min[], PetscReal max[], PetscInt npoints[], InsertMode mode)
346: {
347: PetscReal lmin[] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
348: PetscReal lmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
349: PetscInt i, j, k, bs, b, n_estimate, n_curr, n_new_est, p, n_found, Nfc;
350: const PetscScalar *_coor;
351: DMSwarmCellDM celldm;
352: DM dm;
353: PetscReal dx[3];
354: PetscInt _npoints[] = {0, 0, 1};
355: Vec pos;
356: PetscScalar *_pos;
357: PetscReal *swarm_coor;
358: PetscInt *swarm_cellid;
359: PetscSF sfcell = NULL;
360: const PetscSFNode *LA_sfcell;
361: const char **coordFields, *cellid;
363: PetscFunctionBegin;
364: DMSWARMPICVALID(sw);
365: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
366: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
367: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
369: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
370: PetscCall(DMGetLocalBoundingBox(dm, lmin, lmax));
371: PetscCall(DMGetCoordinateDim(dm, &bs));
373: for (b = 0; b < bs; b++) {
374: if (npoints[b] > 1) {
375: dx[b] = (max[b] - min[b]) / ((PetscReal)(npoints[b] - 1));
376: } else {
377: dx[b] = 0.0;
378: }
379: _npoints[b] = npoints[b];
380: }
382: /* determine number of points living in the bounding box */
383: n_estimate = 0;
384: for (k = 0; k < _npoints[2]; k++) {
385: for (j = 0; j < _npoints[1]; j++) {
386: for (i = 0; i < _npoints[0]; i++) {
387: PetscReal xp[] = {0.0, 0.0, 0.0};
388: PetscInt ijk[3];
389: PetscBool point_inside = PETSC_TRUE;
391: ijk[0] = i;
392: ijk[1] = j;
393: ijk[2] = k;
394: for (b = 0; b < bs; b++) xp[b] = min[b] + ijk[b] * dx[b];
395: for (b = 0; b < bs; b++) {
396: if (xp[b] < lmin[b]) point_inside = PETSC_FALSE;
397: if (xp[b] > lmax[b]) point_inside = PETSC_FALSE;
398: }
399: if (point_inside) n_estimate++;
400: }
401: }
402: }
404: /* create candidate list */
405: PetscCall(VecCreate(PETSC_COMM_SELF, &pos));
406: PetscCall(VecSetSizes(pos, bs * n_estimate, PETSC_DECIDE));
407: PetscCall(VecSetBlockSize(pos, bs));
408: PetscCall(VecSetFromOptions(pos));
409: PetscCall(VecGetArray(pos, &_pos));
411: n_estimate = 0;
412: for (k = 0; k < _npoints[2]; k++) {
413: for (j = 0; j < _npoints[1]; j++) {
414: for (i = 0; i < _npoints[0]; i++) {
415: PetscReal xp[] = {0.0, 0.0, 0.0};
416: PetscInt ijk[3];
417: PetscBool point_inside = PETSC_TRUE;
419: ijk[0] = i;
420: ijk[1] = j;
421: ijk[2] = k;
422: for (b = 0; b < bs; b++) xp[b] = min[b] + ijk[b] * dx[b];
423: for (b = 0; b < bs; b++) {
424: if (xp[b] < lmin[b]) point_inside = PETSC_FALSE;
425: if (xp[b] > lmax[b]) point_inside = PETSC_FALSE;
426: }
427: if (point_inside) {
428: for (b = 0; b < bs; b++) _pos[bs * n_estimate + b] = xp[b];
429: n_estimate++;
430: }
431: }
432: }
433: }
434: PetscCall(VecRestoreArray(pos, &_pos));
436: /* locate points */
437: PetscCall(DMLocatePoints(dm, pos, DM_POINTLOCATION_NONE, &sfcell));
438: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
439: n_found = 0;
440: for (p = 0; p < n_estimate; p++) {
441: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) n_found++;
442: }
444: /* adjust size */
445: if (mode == ADD_VALUES) {
446: PetscCall(DMSwarmGetLocalSize(sw, &n_curr));
447: n_new_est = n_curr + n_found;
448: PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
449: }
450: if (mode == INSERT_VALUES) {
451: n_curr = 0;
452: n_new_est = n_found;
453: PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
454: }
456: /* initialize new coords, cell owners, pid */
457: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
458: PetscCall(VecGetArrayRead(pos, &_coor));
459: PetscCall(DMSwarmGetField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
460: PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
461: n_found = 0;
462: for (p = 0; p < n_estimate; p++) {
463: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
464: for (b = 0; b < bs; b++) swarm_coor[bs * (n_curr + n_found) + b] = PetscRealPart(_coor[bs * p + b]);
465: swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
466: n_found++;
467: }
468: }
469: PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
470: PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
471: PetscCall(VecRestoreArrayRead(pos, &_coor));
473: PetscCall(PetscSFDestroy(&sfcell));
474: PetscCall(VecDestroy(&pos));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: /*@
479: DMSwarmSetPointCoordinates - Set point coordinates in a `DMSWARM` from a user defined list
481: Collective
483: Input Parameters:
484: + sw - the `DMSWARM`
485: . npoints - the number of points to insert
486: . coor - the coordinate values
487: . 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
488: - mode - indicates whether to append points to the swarm (`ADD_VALUES`), or over-ride existing points (`INSERT_VALUES`)
490: Level: beginner
492: Notes:
493: If the user has specified `redundant` as `PETSC_FALSE`, the cell `DM` will attempt to locate the coordinates provided by `coor` within
494: its sub-domain. If they any values within `coor` are not located in the sub-domain, they will be ignored and will not get
495: added to the `DMSWARM`.
497: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`, `DMSwarmSetPointsUniformCoordinates()`
498: @*/
499: PetscErrorCode DMSwarmSetPointCoordinates(DM sw, PetscInt npoints, PetscReal coor[], PetscBool redundant, InsertMode mode)
500: {
501: PetscReal gmin[] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
502: PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
503: PetscInt i, N, bs, b, n_estimate, n_curr, n_new_est, p, n_found;
504: Vec coorlocal;
505: const PetscScalar *_coor;
506: DMSwarmCellDM celldm;
507: DM dm;
508: Vec pos;
509: PetscScalar *_pos;
510: PetscReal *swarm_coor;
511: PetscInt *swarm_cellid;
512: PetscSF sfcell = NULL;
513: const PetscSFNode *LA_sfcell;
514: PetscReal *my_coor;
515: PetscInt my_npoints, Nfc;
516: PetscMPIInt rank;
517: MPI_Comm comm;
518: const char **coordFields, *cellid;
520: PetscFunctionBegin;
521: DMSWARMPICVALID(sw);
522: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
523: PetscCallMPI(MPI_Comm_rank(comm, &rank));
525: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
526: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
527: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
529: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
530: PetscCall(DMGetCoordinatesLocal(dm, &coorlocal));
531: PetscCall(VecGetSize(coorlocal, &N));
532: PetscCall(VecGetBlockSize(coorlocal, &bs));
533: N = N / bs;
534: PetscCall(VecGetArrayRead(coorlocal, &_coor));
535: for (i = 0; i < N; i++) {
536: for (b = 0; b < bs; b++) {
537: gmin[b] = PetscMin(gmin[b], PetscRealPart(_coor[bs * i + b]));
538: gmax[b] = PetscMax(gmax[b], PetscRealPart(_coor[bs * i + b]));
539: }
540: }
541: PetscCall(VecRestoreArrayRead(coorlocal, &_coor));
543: /* broadcast points from rank 0 if requested */
544: if (redundant) {
545: PetscMPIInt imy;
547: my_npoints = npoints;
548: PetscCallMPI(MPI_Bcast(&my_npoints, 1, MPIU_INT, 0, comm));
550: if (rank > 0) { /* allocate space */
551: PetscCall(PetscMalloc1(bs * my_npoints, &my_coor));
552: } else {
553: my_coor = coor;
554: }
555: PetscCall(PetscMPIIntCast(bs * my_npoints, &imy));
556: PetscCallMPI(MPI_Bcast(my_coor, imy, MPIU_REAL, 0, comm));
557: } else {
558: my_npoints = npoints;
559: my_coor = coor;
560: }
562: /* determine the number of points living in the bounding box */
563: n_estimate = 0;
564: for (i = 0; i < my_npoints; i++) {
565: PetscBool point_inside = PETSC_TRUE;
567: for (b = 0; b < bs; b++) {
568: if (my_coor[bs * i + b] < gmin[b]) point_inside = PETSC_FALSE;
569: if (my_coor[bs * i + b] > gmax[b]) point_inside = PETSC_FALSE;
570: }
571: if (point_inside) n_estimate++;
572: }
574: /* create candidate list */
575: PetscCall(VecCreate(PETSC_COMM_SELF, &pos));
576: PetscCall(VecSetSizes(pos, bs * n_estimate, PETSC_DECIDE));
577: PetscCall(VecSetBlockSize(pos, bs));
578: PetscCall(VecSetFromOptions(pos));
579: PetscCall(VecGetArray(pos, &_pos));
581: n_estimate = 0;
582: for (i = 0; i < my_npoints; i++) {
583: PetscBool point_inside = PETSC_TRUE;
585: for (b = 0; b < bs; b++) {
586: if (my_coor[bs * i + b] < gmin[b]) point_inside = PETSC_FALSE;
587: if (my_coor[bs * i + b] > gmax[b]) point_inside = PETSC_FALSE;
588: }
589: if (point_inside) {
590: for (b = 0; b < bs; b++) _pos[bs * n_estimate + b] = my_coor[bs * i + b];
591: n_estimate++;
592: }
593: }
594: PetscCall(VecRestoreArray(pos, &_pos));
596: /* locate points */
597: PetscCall(DMLocatePoints(dm, pos, DM_POINTLOCATION_NONE, &sfcell));
599: PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
600: n_found = 0;
601: for (p = 0; p < n_estimate; p++) {
602: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) n_found++;
603: }
605: /* adjust size */
606: if (mode == ADD_VALUES) {
607: PetscCall(DMSwarmGetLocalSize(sw, &n_curr));
608: n_new_est = n_curr + n_found;
609: PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
610: }
611: if (mode == INSERT_VALUES) {
612: n_curr = 0;
613: n_new_est = n_found;
614: PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
615: }
617: /* initialize new coords, cell owners, pid */
618: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
619: PetscCall(VecGetArrayRead(pos, &_coor));
620: PetscCall(DMSwarmGetField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
621: PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
622: n_found = 0;
623: for (p = 0; p < n_estimate; p++) {
624: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
625: for (b = 0; b < bs; b++) swarm_coor[bs * (n_curr + n_found) + b] = PetscRealPart(_coor[bs * p + b]);
626: swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
627: n_found++;
628: }
629: }
630: PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
631: PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
632: PetscCall(VecRestoreArrayRead(pos, &_coor));
634: if (redundant) {
635: if (rank > 0) PetscCall(PetscFree(my_coor));
636: }
637: PetscCall(PetscSFDestroy(&sfcell));
638: PetscCall(VecDestroy(&pos));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM, DM, DMSwarmPICLayoutType, PetscInt);
643: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM, DM, DMSwarmPICLayoutType, PetscInt);
645: /*@
646: DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell
648: Not Collective
650: Input Parameters:
651: + dm - the `DMSWARM`
652: . layout_type - method used to fill each cell with the cell `DM`
653: - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type)
655: Level: beginner
657: Notes:
658: The insert method will reset any previous defined points within the `DMSWARM`.
660: When using a `DMDA` both 2D and 3D are supported for all layout types provided you are using `DMDA_ELEMENT_Q1`.
662: When using a `DMPLEX` the following case are supported\:
663: .vb
664: (i) DMSWARMPIC_LAYOUT_REGULAR: 2D (triangle),
665: (ii) DMSWARMPIC_LAYOUT_GAUSS: 2D and 3D provided the cell is a tri/tet or a quad/hex,
666: (iii) DMSWARMPIC_LAYOUT_SUBDIVISION: 2D and 3D for quad/hex and 2D tri.
667: .ve
669: .seealso: `DMSWARM`, `DMSwarmPICLayoutType`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
670: @*/
671: PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm, DMSwarmPICLayoutType layout_type, PetscInt fill_param)
672: {
673: DM celldm;
674: PetscBool isDA, isPLEX;
676: PetscFunctionBegin;
677: DMSWARMPICVALID(dm);
678: PetscCall(DMSwarmGetCellDM(dm, &celldm));
679: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isDA));
680: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isPLEX));
681: if (isDA) {
682: PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA(dm, celldm, layout_type, fill_param));
683: } else if (isPLEX) {
684: PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX(dm, celldm, layout_type, fill_param));
685: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }
689: extern PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *);
691: /*@C
692: DMSwarmSetPointCoordinatesCellwise - Insert point coordinates (defined over the reference cell) within each cell
694: Not Collective
696: Input Parameters:
697: + dm - the `DMSWARM`
698: . npoints - the number of points to insert in each cell
699: - xi - the coordinates (defined in the local coordinate system for each cell) to insert
701: Level: beginner
703: Notes:
704: The method will reset any previous defined points within the `DMSWARM`.
705: Only supported for `DMPLEX`. If you are using a `DMDA` it is recommended to either use
706: `DMSwarmInsertPointsUsingCellDM()`, or extract and set the coordinates yourself the following code
707: .vb
708: PetscReal *coor;
709: const char *coordname;
710: DMSwarmGetCoordinateField(dm, &coordname);
711: DMSwarmGetField(dm,coordname,NULL,NULL,(void**)&coor);
712: // user code to define the coordinates here
713: DMSwarmRestoreField(dm,coordname,NULL,NULL,(void**)&coor);
714: .ve
716: .seealso: `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmInsertPointsUsingCellDM()`
717: @*/
718: PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm, PetscInt npoints, PetscReal xi[])
719: {
720: DM celldm;
721: PetscBool isDA, isPLEX;
723: PetscFunctionBegin;
724: DMSWARMPICVALID(dm);
725: PetscCall(DMSwarmGetCellDM(dm, &celldm));
726: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isDA));
727: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isPLEX));
728: PetscCheck(!isDA, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMPLEX. Recommended you use DMSwarmInsertPointsUsingCellDM()");
729: if (isPLEX) {
730: PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, npoints, xi));
731: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
732: PetscFunctionReturn(PETSC_SUCCESS);
733: }
735: /*@
736: DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM
738: Not Collective
740: Input Parameter:
741: . sw - the `DMSWARM`
743: Output Parameters:
744: + ncells - the number of cells in the cell `DM` (optional argument, pass `NULL` to ignore)
745: - count - array of length ncells containing the number of points per cell
747: Level: beginner
749: Notes:
750: The array count is allocated internally and must be free'd by the user.
752: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
753: @*/
754: PetscErrorCode DMSwarmCreatePointPerCellCount(DM sw, PetscInt *ncells, PetscInt **count)
755: {
756: DMSwarmCellDM celldm;
757: PetscBool isvalid;
758: PetscInt nel;
759: PetscInt *sum;
760: const char *cellid;
762: PetscFunctionBegin;
763: PetscCall(DMSwarmSortGetIsValid(sw, &isvalid));
764: nel = 0;
765: if (isvalid) {
766: PetscInt e;
768: PetscCall(DMSwarmSortGetSizes(sw, &nel, NULL));
770: PetscCall(PetscMalloc1(nel, &sum));
771: for (e = 0; e < nel; e++) PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, e, &sum[e]));
772: } else {
773: DM dm;
774: PetscBool isda, isplex, isshell;
775: PetscInt p, npoints;
776: PetscInt *swarm_cellid;
778: /* get the number of cells */
779: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
780: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
781: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
782: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
783: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMSHELL, &isshell));
784: if (isda) {
785: PetscInt _nel, _npe;
786: const PetscInt *_element;
788: PetscCall(DMDAGetElements(dm, &_nel, &_npe, &_element));
789: nel = _nel;
790: PetscCall(DMDARestoreElements(dm, &_nel, &_npe, &_element));
791: } else if (isplex) {
792: PetscInt ps, pe;
794: PetscCall(DMPlexGetHeightStratum(dm, 0, &ps, &pe));
795: nel = pe - ps;
796: } else if (isshell) {
797: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
799: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
800: if (method_DMShellGetNumberOfCells) {
801: PetscCall(method_DMShellGetNumberOfCells(dm, &nel));
802: } else
803: 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);");
804: } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
806: PetscCall(PetscMalloc1(nel, &sum));
807: PetscCall(PetscArrayzero(sum, nel));
808: PetscCall(DMSwarmGetLocalSize(sw, &npoints));
809: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
810: PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
811: for (p = 0; p < npoints; p++) {
812: if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) sum[swarm_cellid[p]]++;
813: }
814: PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
815: }
816: if (ncells) *ncells = nel;
817: *count = sum;
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: /*@
822: DMSwarmGetNumSpecies - Get the number of particle species
824: Not Collective
826: Input Parameter:
827: . sw - the `DMSWARM`
829: Output Parameters:
830: . Ns - the number of species
832: Level: intermediate
834: .seealso: `DMSWARM`, `DMSwarmSetNumSpecies()`, `DMSwarmSetType()`, `DMSwarmType`
835: @*/
836: PetscErrorCode DMSwarmGetNumSpecies(DM sw, PetscInt *Ns)
837: {
838: DM_Swarm *swarm = (DM_Swarm *)sw->data;
840: PetscFunctionBegin;
841: *Ns = swarm->Ns;
842: PetscFunctionReturn(PETSC_SUCCESS);
843: }
845: /*@
846: DMSwarmSetNumSpecies - Set the number of particle species
848: Not Collective
850: Input Parameters:
851: + sw - the `DMSWARM`
852: - Ns - the number of species
854: Level: intermediate
856: .seealso: `DMSWARM`, `DMSwarmGetNumSpecies()`, `DMSwarmSetType()`, `DMSwarmType`
857: @*/
858: PetscErrorCode DMSwarmSetNumSpecies(DM sw, PetscInt Ns)
859: {
860: DM_Swarm *swarm = (DM_Swarm *)sw->data;
862: PetscFunctionBegin;
863: swarm->Ns = Ns;
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@C
868: DMSwarmGetCoordinateFunction - Get the function setting initial particle positions, if it exists
870: Not Collective
872: Input Parameter:
873: . sw - the `DMSWARM`
875: Output Parameter:
876: . coordFunc - the function setting initial particle positions, or `NULL`, see `PetscSimplePointFn` for the calling sequence
878: Level: intermediate
880: .seealso: `DMSWARM`, `DMSwarmSetCoordinateFunction()`, `DMSwarmGetVelocityFunction()`, `DMSwarmInitializeCoordinates()`, `PetscSimplePointFn`
881: @*/
882: PetscErrorCode DMSwarmGetCoordinateFunction(DM sw, PetscSimplePointFn **coordFunc)
883: {
884: DM_Swarm *swarm = (DM_Swarm *)sw->data;
886: PetscFunctionBegin;
888: PetscAssertPointer(coordFunc, 2);
889: *coordFunc = swarm->coordFunc;
890: PetscFunctionReturn(PETSC_SUCCESS);
891: }
893: /*@C
894: DMSwarmSetCoordinateFunction - Set the function setting initial particle positions
896: Not Collective
898: Input Parameters:
899: + sw - the `DMSWARM`
900: - coordFunc - the function setting initial particle positions, see `PetscSimplePointFn` for the calling sequence
902: Level: intermediate
904: .seealso: `DMSWARM`, `DMSwarmGetCoordinateFunction()`, `DMSwarmSetVelocityFunction()`, `DMSwarmInitializeCoordinates()`, `PetscSimplePointFn`
905: @*/
906: PetscErrorCode DMSwarmSetCoordinateFunction(DM sw, PetscSimplePointFn *coordFunc)
907: {
908: DM_Swarm *swarm = (DM_Swarm *)sw->data;
910: PetscFunctionBegin;
913: swarm->coordFunc = coordFunc;
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: /*@C
918: DMSwarmGetVelocityFunction - Get the function setting initial particle velocities, if it exists
920: Not Collective
922: Input Parameter:
923: . sw - the `DMSWARM`
925: Output Parameter:
926: . velFunc - the function setting initial particle velocities, or `NULL`, see `PetscSimplePointFn` for the calling sequence
928: Level: intermediate
930: .seealso: `DMSWARM`, `DMSwarmSetVelocityFunction()`, `DMSwarmGetCoordinateFunction()`, `DMSwarmInitializeVelocities()`, `PetscSimplePointFn`
931: @*/
932: PetscErrorCode DMSwarmGetVelocityFunction(DM sw, PetscSimplePointFn **velFunc)
933: {
934: DM_Swarm *swarm = (DM_Swarm *)sw->data;
936: PetscFunctionBegin;
938: PetscAssertPointer(velFunc, 2);
939: *velFunc = swarm->velFunc;
940: PetscFunctionReturn(PETSC_SUCCESS);
941: }
943: /*@C
944: DMSwarmSetVelocityFunction - Set the function setting initial particle velocities
946: Not Collective
948: Input Parameters:
949: + sw - the `DMSWARM`
950: - velFunc - the function setting initial particle velocities, see `PetscSimplePointFn` for the calling sequence
952: Level: intermediate
954: .seealso: `DMSWARM`, `DMSwarmGetVelocityFunction()`, `DMSwarmSetCoordinateFunction()`, `DMSwarmInitializeVelocities()`, `PetscSimplePointFn`
955: @*/
956: PetscErrorCode DMSwarmSetVelocityFunction(DM sw, PetscSimplePointFn *velFunc)
957: {
958: DM_Swarm *swarm = (DM_Swarm *)sw->data;
960: PetscFunctionBegin;
963: swarm->velFunc = velFunc;
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: /*@C
968: DMSwarmComputeLocalSize - Compute the local number and distribution of particles based upon a density function
970: Not Collective
972: Input Parameters:
973: + sw - The `DMSWARM`
974: . N - The target number of particles
975: - density - The density field for the particle layout, normalized to unity
977: Level: advanced
979: Note:
980: One particle will be created for each species.
982: .seealso: `DMSWARM`, `DMSwarmComputeLocalSizeFromOptions()`
983: @*/
984: PetscErrorCode DMSwarmComputeLocalSize(DM sw, PetscInt N, PetscProbFunc density)
985: {
986: DM dm;
987: DMSwarmCellDM celldm;
988: PetscQuadrature quad;
989: const PetscReal *xq, *wq;
990: PetscReal *n_int;
991: PetscInt *npc_s, *swarm_cellid, Ni;
992: PetscReal gmin[3], gmax[3], xi0[3];
993: PetscInt Ns, cStart, cEnd, c, dim, d, Nq, q, Np = 0, p, s;
994: PetscBool simplex;
995: const char *cellid;
997: PetscFunctionBegin;
998: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
999: PetscCall(DMSwarmGetCellDM(sw, &dm));
1000: PetscCall(DMGetDimension(dm, &dim));
1001: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
1002: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1003: PetscCall(DMPlexIsSimplex(dm, &simplex));
1004: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1005: if (simplex) PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 5, -1.0, 1.0, &quad));
1006: else PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 5, -1.0, 1.0, &quad));
1007: PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xq, &wq));
1008: PetscCall(PetscCalloc2(Ns, &n_int, (cEnd - cStart) * Ns, &npc_s));
1009: /* Integrate the density function to get the number of particles in each cell */
1010: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1011: for (c = 0; c < cEnd - cStart; ++c) {
1012: const PetscInt cell = c + cStart;
1013: PetscReal v0[3], J[9], invJ[9], detJ, detJp = 2. / (gmax[0] - gmin[0]), xr[3], den;
1015: /* Have to transform quadrature points/weights to cell domain */
1016: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
1017: PetscCall(PetscArrayzero(n_int, Ns));
1018: for (q = 0; q < Nq; ++q) {
1019: CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], xr);
1020: /* Have to transform mesh to domain of definition of PDF, [-1, 1], and weight PDF by |J|/2 */
1021: xr[0] = detJp * (xr[0] - gmin[0]) - 1.;
1023: for (s = 0; s < Ns; ++s) {
1024: PetscCall(density(xr, NULL, &den));
1025: n_int[s] += (detJp * den) * (detJ * wq[q]) / (PetscReal)Ns;
1026: }
1027: }
1028: for (s = 0; s < Ns; ++s) {
1029: Ni = N;
1030: npc_s[c * Ns + s] += (PetscInt)(Ni * n_int[s]);
1031: Np += npc_s[c * Ns + s];
1032: }
1033: }
1034: PetscCall(PetscQuadratureDestroy(&quad));
1035: PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
1036: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1037: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
1038: PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
1039: for (c = 0, p = 0; c < cEnd - cStart; ++c) {
1040: for (s = 0; s < Ns; ++s) {
1041: for (q = 0; q < npc_s[c * Ns + s]; ++q, ++p) swarm_cellid[p] = c;
1042: }
1043: }
1044: PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
1045: PetscCall(PetscFree2(n_int, npc_s));
1046: PetscFunctionReturn(PETSC_SUCCESS);
1047: }
1049: /*@
1050: DMSwarmComputeLocalSizeFromOptions - Compute the local number and distribution of particles based upon a density function determined by options
1052: Not Collective
1054: Input Parameter:
1055: . sw - The `DMSWARM`
1057: Level: advanced
1059: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`
1060: @*/
1061: PetscErrorCode DMSwarmComputeLocalSizeFromOptions(DM sw)
1062: {
1063: PetscProbFunc pdf;
1064: const char *prefix;
1065: char funcname[PETSC_MAX_PATH_LEN];
1066: PetscInt *N, Ns, dim, n;
1067: PetscBool flg;
1068: PetscMPIInt size, rank;
1070: PetscFunctionBegin;
1071: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sw), &size));
1072: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1073: PetscCall(PetscCalloc1(size, &N));
1074: PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
1075: n = size;
1076: PetscCall(PetscOptionsIntArray("-dm_swarm_num_particles", "The target number of particles", "", N, &n, NULL));
1077: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1078: PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
1079: if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
1080: PetscCall(PetscOptionsString("-dm_swarm_coordinate_function", "Function to determine particle coordinates", "DMSwarmSetCoordinateFunction", funcname, funcname, sizeof(funcname), &flg));
1081: PetscOptionsEnd();
1082: if (flg) {
1083: PetscSimplePointFn *coordFunc;
1085: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1086: PetscCall(PetscDLSym(NULL, funcname, (void **)&coordFunc));
1087: PetscCheck(coordFunc, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
1088: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1089: PetscCall(DMSwarmSetLocalSizes(sw, N[rank] * Ns, 0));
1090: PetscCall(DMSwarmSetCoordinateFunction(sw, coordFunc));
1091: } else {
1092: PetscCall(DMGetDimension(sw, &dim));
1093: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
1094: PetscCall(PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_coordinate_density", &pdf, NULL, NULL));
1095: PetscCall(DMSwarmComputeLocalSize(sw, N[rank], pdf));
1096: }
1097: PetscCall(PetscFree(N));
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: /*@
1102: DMSwarmInitializeCoordinates - Determine the initial coordinates of particles for a PIC method
1104: Not Collective
1106: Input Parameter:
1107: . sw - The `DMSWARM`
1109: Level: advanced
1111: Note:
1112: Currently, we randomly place particles in their assigned cell
1114: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeVelocities()`
1115: @*/
1116: PetscErrorCode DMSwarmInitializeCoordinates(DM sw)
1117: {
1118: DMSwarmCellDM celldm;
1119: PetscSimplePointFn *coordFunc;
1120: PetscScalar *weight;
1121: PetscReal *x;
1122: PetscInt *species;
1123: void *ctx;
1124: PetscBool removePoints = PETSC_TRUE;
1125: PetscDataType dtype;
1126: PetscInt Nfc, Np, p, Ns, dim, d, bs;
1127: const char **coordFields;
1129: PetscFunctionBeginUser;
1130: PetscCall(DMGetDimension(sw, &dim));
1131: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1132: PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1133: PetscCall(DMSwarmGetCoordinateFunction(sw, &coordFunc));
1135: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1136: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1137: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
1139: PetscCall(DMSwarmGetField(sw, coordFields[0], &bs, &dtype, (void **)&x));
1140: PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&weight));
1141: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1142: if (coordFunc) {
1143: PetscCall(DMGetApplicationContext(sw, &ctx));
1144: for (p = 0; p < Np; ++p) {
1145: PetscScalar X[3];
1147: PetscCall((*coordFunc)(dim, 0., NULL, p, X, ctx));
1148: for (d = 0; d < dim; ++d) x[p * dim + d] = PetscRealPart(X[d]);
1149: weight[p] = 1.0;
1150: species[p] = p % Ns;
1151: }
1152: } else {
1153: DM dm;
1154: PetscRandom rnd;
1155: PetscReal xi0[3];
1156: PetscInt cStart, cEnd, c;
1158: PetscCall(DMSwarmGetCellDM(sw, &dm));
1159: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1160: PetscCall(DMGetApplicationContext(sw, &ctx));
1162: /* Set particle position randomly in cell, set weights to 1 */
1163: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
1164: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
1165: PetscCall(PetscRandomSetFromOptions(rnd));
1166: PetscCall(DMSwarmSortGetAccess(sw));
1167: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1168: for (c = cStart; c < cEnd; ++c) {
1169: PetscReal v0[3], J[9], invJ[9], detJ;
1170: PetscInt *pidx, Npc, q;
1172: PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1173: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
1174: for (q = 0; q < Npc; ++q) {
1175: const PetscInt p = pidx[q];
1176: PetscReal xref[3];
1178: for (d = 0; d < dim; ++d) PetscCall(PetscRandomGetValueReal(rnd, &xref[d]));
1179: CoordinatesRefToReal(dim, dim, xi0, v0, J, xref, &x[p * dim]);
1181: weight[p] = 1.0 / Np;
1182: species[p] = p % Ns;
1183: }
1184: PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1185: }
1186: PetscCall(PetscRandomDestroy(&rnd));
1187: PetscCall(DMSwarmSortRestoreAccess(sw));
1188: }
1189: PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&x));
1190: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1191: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1193: PetscCall(DMSwarmMigrate(sw, removePoints));
1194: PetscCall(DMLocalizeCoordinates(sw));
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: /*@C
1199: DMSwarmInitializeVelocities - Set the initial velocities of particles using a distribution.
1201: Collective
1203: Input Parameters:
1204: + sw - The `DMSWARM` object
1205: . sampler - A function which uniformly samples the velocity PDF
1206: - v0 - The velocity scale for nondimensionalization for each species
1208: Level: advanced
1210: Note:
1211: 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.
1213: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeCoordinates()`, `DMSwarmInitializeVelocitiesFromOptions()`
1214: @*/
1215: PetscErrorCode DMSwarmInitializeVelocities(DM sw, PetscProbFunc sampler, const PetscReal v0[])
1216: {
1217: PetscSimplePointFn *velFunc;
1218: PetscReal *v;
1219: PetscInt *species;
1220: void *ctx;
1221: PetscInt dim, Np, p;
1223: PetscFunctionBegin;
1224: PetscCall(DMSwarmGetVelocityFunction(sw, &velFunc));
1226: PetscCall(DMGetDimension(sw, &dim));
1227: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1228: PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
1229: PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1230: if (v0[0] == 0.) {
1231: PetscCall(PetscArrayzero(v, Np * dim));
1232: } else if (velFunc) {
1233: PetscCall(DMGetApplicationContext(sw, &ctx));
1234: for (p = 0; p < Np; ++p) {
1235: PetscInt s = species[p], d;
1236: PetscScalar vel[3];
1238: PetscCall((*velFunc)(dim, 0., NULL, p, vel, ctx));
1239: for (d = 0; d < dim; ++d) v[p * dim + d] = (v0[s] / v0[0]) * PetscRealPart(vel[d]);
1240: }
1241: } else {
1242: PetscRandom rnd;
1244: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd));
1245: PetscCall(PetscRandomSetInterval(rnd, 0, 1.));
1246: PetscCall(PetscRandomSetFromOptions(rnd));
1248: for (p = 0; p < Np; ++p) {
1249: PetscInt s = species[p], d;
1250: PetscReal a[3], vel[3];
1252: for (d = 0; d < dim; ++d) PetscCall(PetscRandomGetValueReal(rnd, &a[d]));
1253: PetscCall(sampler(a, NULL, vel));
1254: for (d = 0; d < dim; ++d) v[p * dim + d] = (v0[s] / v0[0]) * vel[d];
1255: }
1256: PetscCall(PetscRandomDestroy(&rnd));
1257: }
1258: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1259: PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: /*@
1264: DMSwarmInitializeVelocitiesFromOptions - Set the initial velocities of particles using a distribution determined from options.
1266: Collective
1268: Input Parameters:
1269: + sw - The `DMSWARM` object
1270: - v0 - The velocity scale for nondimensionalization for each species
1272: Level: advanced
1274: .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeCoordinates()`, `DMSwarmInitializeVelocities()`
1275: @*/
1276: PetscErrorCode DMSwarmInitializeVelocitiesFromOptions(DM sw, const PetscReal v0[])
1277: {
1278: PetscProbFunc sampler;
1279: PetscInt dim;
1280: const char *prefix;
1281: char funcname[PETSC_MAX_PATH_LEN];
1282: PetscBool flg;
1284: PetscFunctionBegin;
1285: PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
1286: PetscCall(PetscOptionsString("-dm_swarm_velocity_function", "Function to determine particle velocities", "DMSwarmSetVelocityFunction", funcname, funcname, sizeof(funcname), &flg));
1287: PetscOptionsEnd();
1288: if (flg) {
1289: PetscSimplePointFn *velFunc;
1291: PetscCall(PetscDLSym(NULL, funcname, (void **)&velFunc));
1292: PetscCheck(velFunc, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
1293: PetscCall(DMSwarmSetVelocityFunction(sw, velFunc));
1294: }
1295: PetscCall(DMGetDimension(sw, &dim));
1296: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
1297: PetscCall(PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_velocity_density", NULL, NULL, &sampler));
1298: PetscCall(DMSwarmInitializeVelocities(sw, sampler, v0));
1299: PetscFunctionReturn(PETSC_SUCCESS);
1300: }