Actual source code: swarm.c
1: #include "petscdmswarm.h"
2: #define PETSCDM_DLL
3: #include <petsc/private/dmswarmimpl.h>
4: #include <petsc/private/hashsetij.h>
5: #include <petsc/private/petscfeimpl.h>
6: #include <petscviewer.h>
7: #include <petscdraw.h>
8: #include <petscdmplex.h>
9: #include <petscblaslapack.h>
10: #include "../src/dm/impls/swarm/data_bucket.h"
11: #include <petscdmlabel.h>
12: #include <petscsection.h>
14: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
15: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
16: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
18: const char *DMSwarmTypeNames[] = {"basic", "pic", NULL};
19: const char *DMSwarmMigrateTypeNames[] = {"basic", "dmcellnscatter", "dmcellexact", "user", NULL};
20: const char *DMSwarmCollectTypeNames[] = {"basic", "boundingbox", "general", "user", NULL};
21: const char *DMSwarmRemapTypeNames[] = {"none", "pfak", "colella", "DMSwarmRemapType", "DMSWARM_REMAP_", NULL};
22: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
24: const char DMSwarmField_pid[] = "DMSwarm_pid";
25: const char DMSwarmField_rank[] = "DMSwarm_rank";
26: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
28: PetscInt SwarmDataFieldId = -1;
30: #if defined(PETSC_HAVE_HDF5)
31: #include <petscviewerhdf5.h>
33: static PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
34: {
35: DM dm;
36: PetscReal seqval;
37: PetscInt seqnum, bs;
38: PetscBool isseq, ists;
40: PetscFunctionBegin;
41: PetscCall(VecGetDM(v, &dm));
42: PetscCall(VecGetBlockSize(v, &bs));
43: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particle_fields"));
44: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
45: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
46: PetscCall(PetscViewerHDF5IsTimestepping(viewer, &ists));
47: if (ists) PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
48: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
49: PetscCall(VecViewNative(v, viewer));
50: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
51: PetscCall(PetscViewerHDF5PopGroup(viewer));
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
56: {
57: DMSwarmCellDM celldm;
58: Vec coordinates;
59: PetscInt Np, Nfc;
60: PetscBool isseq;
61: const char **coordFields;
63: PetscFunctionBegin;
64: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
65: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
66: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
67: PetscCall(DMSwarmGetSize(dm, &Np));
68: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, coordFields[0], &coordinates));
69: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
70: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
71: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
72: PetscCall(VecViewNative(coordinates, viewer));
73: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
74: PetscCall(PetscViewerHDF5PopGroup(viewer));
75: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, coordFields[0], &coordinates));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
78: #endif
80: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
81: {
82: DM dm;
83: #if defined(PETSC_HAVE_HDF5)
84: PetscBool ishdf5;
85: #endif
87: PetscFunctionBegin;
88: PetscCall(VecGetDM(v, &dm));
89: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
90: #if defined(PETSC_HAVE_HDF5)
91: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
92: if (ishdf5) {
93: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
96: #endif
97: PetscCall(VecViewNative(v, viewer));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /*@C
102: DMSwarmVectorGetField - Gets the fields from which to define a `Vec` object
103: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
105: Not collective
107: Input Parameter:
108: . sw - a `DMSWARM`
110: Output Parameters:
111: + Nf - the number of fields
112: - fieldnames - the textual name given to each registered field, or NULL if it has not been set
114: Level: beginner
116: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
117: @*/
118: PetscErrorCode DMSwarmVectorGetField(DM sw, PetscInt *Nf, const char **fieldnames[])
119: {
120: DMSwarmCellDM celldm;
122: PetscFunctionBegin;
124: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
125: PetscCall(DMSwarmCellDMGetFields(celldm, Nf, fieldnames));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*@
130: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
131: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
133: Collective
135: Input Parameters:
136: + dm - a `DMSWARM`
137: - fieldname - the textual name given to each registered field
139: Level: beginner
141: Notes:
142: The field with name `fieldname` must be defined as having a data type of `PetscScalar`.
144: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
145: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
147: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineFields()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
148: @*/
149: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
150: {
151: PetscFunctionBegin;
152: PetscCall(DMSwarmVectorDefineFields(dm, 1, &fieldname));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: /*@C
157: DMSwarmVectorDefineFields - Sets the fields from which to define a `Vec` object
158: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
160: Collective, No Fortran support
162: Input Parameters:
163: + sw - a `DMSWARM`
164: . Nf - the number of fields
165: - fieldnames - the textual name given to each registered field
167: Level: beginner
169: Notes:
170: Each field with name in `fieldnames` must be defined as having a data type of `PetscScalar`.
172: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
173: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
175: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
176: @*/
177: PetscErrorCode DMSwarmVectorDefineFields(DM sw, PetscInt Nf, const char *fieldnames[])
178: {
179: DM_Swarm *swarm = (DM_Swarm *)sw->data;
180: DMSwarmCellDM celldm;
182: PetscFunctionBegin;
184: if (fieldnames) PetscAssertPointer(fieldnames, 3);
185: if (!swarm->issetup) PetscCall(DMSetUp(sw));
186: PetscCheck(Nf >= 0, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields must be non-negative, not %" PetscInt_FMT, Nf);
187: // Create a dummy cell DM if none has been specified (I think we should not support this mode)
188: if (!swarm->activeCellDM) {
189: DM dm;
190: DMSwarmCellDM celldm;
192: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), &dm));
193: PetscCall(DMSetType(dm, DMSHELL));
194: PetscCall(PetscObjectSetName((PetscObject)dm, "dummy"));
195: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 0, NULL, &celldm));
196: PetscCall(DMDestroy(&dm));
197: PetscCall(DMSwarmAddCellDM(sw, celldm));
198: PetscCall(DMSwarmCellDMDestroy(&celldm));
199: PetscCall(DMSwarmSetCellDMActive(sw, "dummy"));
200: }
201: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
202: for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscFree(celldm->dmFields[f]));
203: PetscCall(PetscFree(celldm->dmFields));
205: celldm->Nf = Nf;
206: PetscCall(PetscMalloc1(Nf, &celldm->dmFields));
207: for (PetscInt f = 0; f < Nf; ++f) {
208: PetscDataType type;
210: // Check all fields are of type PETSC_REAL or PETSC_SCALAR
211: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], NULL, &type));
212: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
213: PetscCall(PetscStrallocpy(fieldnames[f], (char **)&celldm->dmFields[f]));
214: }
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /* requires DMSwarmDefineFieldVector has been called */
219: static PetscErrorCode DMCreateGlobalVector_Swarm(DM sw, Vec *vec)
220: {
221: DM_Swarm *swarm = (DM_Swarm *)sw->data;
222: DMSwarmCellDM celldm;
223: Vec x;
224: char name[PETSC_MAX_PATH_LEN];
225: PetscInt bs = 0, n;
227: PetscFunctionBegin;
228: if (!swarm->issetup) PetscCall(DMSetUp(sw));
229: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
230: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
231: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
233: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
234: for (PetscInt f = 0; f < celldm->Nf; ++f) {
235: PetscInt fbs;
236: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
237: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
238: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
239: bs += fbs;
240: }
241: PetscCall(VecCreate(PetscObjectComm((PetscObject)sw), &x));
242: PetscCall(PetscObjectSetName((PetscObject)x, name));
243: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
244: PetscCall(VecSetBlockSize(x, bs));
245: PetscCall(VecSetDM(x, sw));
246: PetscCall(VecSetFromOptions(x));
247: PetscCall(VecSetOperation(x, VECOP_VIEW, (void (*)(void))VecView_Swarm));
248: *vec = x;
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /* requires DMSwarmDefineFieldVector has been called */
253: static PetscErrorCode DMCreateLocalVector_Swarm(DM sw, Vec *vec)
254: {
255: DM_Swarm *swarm = (DM_Swarm *)sw->data;
256: DMSwarmCellDM celldm;
257: Vec x;
258: char name[PETSC_MAX_PATH_LEN];
259: PetscInt bs = 0, n;
261: PetscFunctionBegin;
262: if (!swarm->issetup) PetscCall(DMSetUp(sw));
263: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
264: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
265: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
267: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
268: for (PetscInt f = 0; f < celldm->Nf; ++f) {
269: PetscInt fbs;
270: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
271: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
272: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
273: bs += fbs;
274: }
275: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
276: PetscCall(PetscObjectSetName((PetscObject)x, name));
277: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
278: PetscCall(VecSetBlockSize(x, bs));
279: PetscCall(VecSetDM(x, sw));
280: PetscCall(VecSetFromOptions(x));
281: *vec = x;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
286: {
287: DM_Swarm *swarm = (DM_Swarm *)dm->data;
288: DMSwarmDataField gfield;
289: PetscInt bs, nlocal, fid = -1, cfid = -2;
290: PetscBool flg;
292: PetscFunctionBegin;
293: /* check vector is an inplace array */
294: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
295: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
296: (void)flg; /* avoid compiler warning */
297: PetscCheck(cfid == fid, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)! %" PetscInt_FMT " != %" PetscInt_FMT, fieldname, cfid, fid);
298: PetscCall(VecGetLocalSize(*vec, &nlocal));
299: PetscCall(VecGetBlockSize(*vec, &bs));
300: PetscCheck(nlocal / bs == swarm->db->L, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
301: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
302: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
303: PetscCall(VecResetArray(*vec));
304: PetscCall(VecDestroy(vec));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
309: {
310: DM_Swarm *swarm = (DM_Swarm *)dm->data;
311: PetscDataType type;
312: PetscScalar *array;
313: PetscInt bs, n, fid;
314: char name[PETSC_MAX_PATH_LEN];
315: PetscMPIInt size;
316: PetscBool iscuda, iskokkos, iship;
318: PetscFunctionBegin;
319: if (!swarm->issetup) PetscCall(DMSetUp(dm));
320: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
321: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
322: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
324: PetscCallMPI(MPI_Comm_size(comm, &size));
325: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
326: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
327: PetscCall(PetscStrcmp(dm->vectype, VECHIP, &iship));
328: PetscCall(VecCreate(comm, vec));
329: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
330: PetscCall(VecSetBlockSize(*vec, bs));
331: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
332: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
333: else if (iship) PetscCall(VecSetType(*vec, VECHIP));
334: else PetscCall(VecSetType(*vec, VECSTANDARD));
335: PetscCall(VecPlaceArray(*vec, array));
337: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
338: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
340: /* Set guard */
341: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
342: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
344: PetscCall(VecSetDM(*vec, dm));
345: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: static PetscErrorCode DMSwarmDestroyVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], Vec *vec)
350: {
351: DM_Swarm *swarm = (DM_Swarm *)sw->data;
352: const PetscScalar *array;
353: PetscInt bs, n, id = 0, cid = -2;
354: PetscBool flg;
356: PetscFunctionBegin;
357: for (PetscInt f = 0; f < Nf; ++f) {
358: PetscInt fid;
360: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
361: id += fid;
362: }
363: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cid, flg));
364: (void)flg; /* avoid compiler warning */
365: PetscCheck(cid == id, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)! %" PetscInt_FMT " != %" PetscInt_FMT, fieldnames[0], cid, id);
366: PetscCall(VecGetLocalSize(*vec, &n));
367: PetscCall(VecGetBlockSize(*vec, &bs));
368: n /= bs;
369: PetscCheck(n == swarm->db->L, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
370: PetscCall(VecGetArrayRead(*vec, &array));
371: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
372: PetscScalar *farray;
373: PetscDataType ftype;
374: PetscInt fbs;
376: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
377: PetscCheck(off + fbs <= bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid blocksize %" PetscInt_FMT " < %" PetscInt_FMT, bs, off + fbs);
378: for (PetscInt i = 0; i < n; ++i) {
379: for (PetscInt b = 0; b < fbs; ++b) farray[i * fbs + b] = array[i * bs + off + b];
380: }
381: off += fbs;
382: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
383: }
384: PetscCall(VecRestoreArrayRead(*vec, &array));
385: PetscCall(VecDestroy(vec));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: static PetscErrorCode DMSwarmCreateVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], MPI_Comm comm, Vec *vec)
390: {
391: DM_Swarm *swarm = (DM_Swarm *)sw->data;
392: PetscScalar *array;
393: PetscInt n, bs = 0, id = 0;
394: char name[PETSC_MAX_PATH_LEN];
396: PetscFunctionBegin;
397: if (!swarm->issetup) PetscCall(DMSetUp(sw));
398: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
399: for (PetscInt f = 0; f < Nf; ++f) {
400: PetscDataType ftype;
401: PetscInt fbs;
403: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
404: PetscCheck(ftype == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
405: bs += fbs;
406: }
408: PetscCall(VecCreate(comm, vec));
409: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
410: PetscCall(VecSetBlockSize(*vec, bs));
411: PetscCall(VecSetType(*vec, sw->vectype));
413: PetscCall(VecGetArrayWrite(*vec, &array));
414: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
415: PetscScalar *farray;
416: PetscDataType ftype;
417: PetscInt fbs;
419: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
420: for (PetscInt i = 0; i < n; ++i) {
421: for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = farray[i * fbs + b];
422: }
423: off += fbs;
424: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
425: }
426: PetscCall(VecRestoreArrayWrite(*vec, &array));
428: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
429: for (PetscInt f = 0; f < Nf; ++f) {
430: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
431: PetscCall(PetscStrlcat(name, fieldnames[f], PETSC_MAX_PATH_LEN));
432: }
433: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
435: for (PetscInt f = 0; f < Nf; ++f) {
436: PetscInt fid;
438: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
439: id += fid;
440: }
441: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, id));
443: PetscCall(VecSetDM(*vec, sw));
444: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
450: \hat f = \sum_i f_i \phi_i
452: and a particle function is given by
454: f = \sum_i w_i \delta(x - x_i)
456: then we want to require that
458: M \hat f = M_p f
460: where the particle mass matrix is given by
462: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
464: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
465: his integral. We allow this with the boolean flag.
466: */
467: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
468: {
469: const char *name = "Mass Matrix";
470: MPI_Comm comm;
471: DMSwarmCellDM celldm;
472: PetscDS prob;
473: PetscSection fsection, globalFSection;
474: PetscHSetIJ ht;
475: PetscLayout rLayout, colLayout;
476: PetscInt *dnz, *onz;
477: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
478: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
479: PetscScalar *elemMat;
480: PetscInt dim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0, totNc = 0;
481: const char **coordFields;
482: PetscReal **coordVals;
483: PetscInt *bs;
485: PetscFunctionBegin;
486: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
487: PetscCall(DMGetCoordinateDim(dmf, &dim));
488: PetscCall(DMGetDS(dmf, &prob));
489: PetscCall(PetscDSGetNumFields(prob, &Nf));
490: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
491: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
492: PetscCall(DMGetLocalSection(dmf, &fsection));
493: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
494: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
495: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
497: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
498: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
499: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
501: PetscCall(PetscLayoutCreate(comm, &colLayout));
502: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
503: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
504: PetscCall(PetscLayoutSetUp(colLayout));
505: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
506: PetscCall(PetscLayoutDestroy(&colLayout));
508: PetscCall(PetscLayoutCreate(comm, &rLayout));
509: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
510: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
511: PetscCall(PetscLayoutSetUp(rLayout));
512: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
513: PetscCall(PetscLayoutDestroy(&rLayout));
515: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
516: PetscCall(PetscHSetIJCreate(&ht));
518: PetscCall(PetscSynchronizedFlush(comm, NULL));
519: for (PetscInt field = 0; field < Nf; ++field) {
520: PetscObject obj;
521: PetscClassId id;
522: PetscInt Nc;
524: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
525: PetscCall(PetscObjectGetClassId(obj, &id));
526: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
527: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
528: totNc += Nc;
529: }
530: /* count non-zeros */
531: PetscCall(DMSwarmSortGetAccess(dmc));
532: for (PetscInt field = 0; field < Nf; ++field) {
533: PetscObject obj;
534: PetscClassId id;
535: PetscInt Nc;
537: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
538: PetscCall(PetscObjectGetClassId(obj, &id));
539: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
540: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
542: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
543: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
544: PetscInt numFIndices, numCIndices;
546: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
547: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
548: maxC = PetscMax(maxC, numCIndices);
549: {
550: PetscHashIJKey key;
551: PetscBool missing;
552: for (PetscInt i = 0; i < numFIndices; ++i) {
553: key.j = findices[i]; /* global column (from Plex) */
554: if (key.j >= 0) {
555: /* Get indices for coarse elements */
556: for (PetscInt j = 0; j < numCIndices; ++j) {
557: for (PetscInt c = 0; c < Nc; ++c) {
558: // TODO Need field offset on particle here
559: key.i = cindices[j] * totNc + c + rStart; /* global cols (from Swarm) */
560: if (key.i < 0) continue;
561: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
562: if (missing) {
563: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
564: else ++onz[key.i - rStart];
565: } else SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
566: }
567: }
568: }
569: }
570: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
571: }
572: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
573: }
574: }
575: PetscCall(PetscHSetIJDestroy(&ht));
576: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
577: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
578: PetscCall(PetscFree2(dnz, onz));
579: PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
580: for (PetscInt field = 0; field < Nf; ++field) {
581: PetscTabulation Tcoarse;
582: PetscObject obj;
583: PetscClassId id;
584: PetscInt Nc;
586: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
587: PetscCall(PetscObjectGetClassId(obj, &id));
588: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
589: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
591: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
592: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
593: PetscInt *findices, *cindices;
594: PetscInt numFIndices, numCIndices;
596: /* TODO: Use DMField instead of assuming affine */
597: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
598: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
599: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
600: for (PetscInt j = 0; j < numCIndices; ++j) {
601: PetscReal xr[8];
602: PetscInt off = 0;
604: for (PetscInt i = 0; i < Nfc; ++i) {
605: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[j] * bs[i] + b];
606: }
607: 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);
608: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, xr, &xi[j * dim]);
609: }
610: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
611: else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
612: /* Get elemMat entries by multiplying by weight */
613: PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
614: for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
615: for (PetscInt j = 0; j < numCIndices; ++j) {
616: for (PetscInt c = 0; c < Nc; ++c) {
617: // TODO Need field offset on particle and field here
618: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
619: elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
620: }
621: }
622: }
623: for (PetscInt j = 0; j < numCIndices; ++j)
624: // TODO Need field offset on particle here
625: for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
626: if (0) PetscCall(DMPrintCellMatrix(cell, name, numCIndices * Nc, numFIndices, elemMat));
627: PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
628: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
629: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
630: PetscCall(PetscTabulationDestroy(&Tcoarse));
631: }
632: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
633: }
634: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
635: PetscCall(DMSwarmSortRestoreAccess(dmc));
636: PetscCall(PetscFree3(v0, J, invJ));
637: PetscCall(PetscFree2(coordVals, bs));
638: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
639: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: /* Returns empty matrix for use with SNES FD */
644: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
645: {
646: Vec field;
647: PetscInt size;
649: PetscFunctionBegin;
650: PetscCall(DMGetGlobalVector(sw, &field));
651: PetscCall(VecGetLocalSize(field, &size));
652: PetscCall(DMRestoreGlobalVector(sw, &field));
653: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
654: PetscCall(MatSetFromOptions(*m));
655: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
656: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
657: PetscCall(MatZeroEntries(*m));
658: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
659: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
660: PetscCall(MatShift(*m, 1.0));
661: PetscCall(MatSetDM(*m, sw));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: /* FEM cols, Particle rows */
666: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
667: {
668: DMSwarmCellDM celldm;
669: PetscSection gsf;
670: PetscInt m, n, Np, bs;
671: void *ctx;
673: PetscFunctionBegin;
674: PetscCall(DMSwarmGetCellDMActive(dmCoarse, &celldm));
675: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)dmCoarse), PETSC_ERR_USER, "Active cell DM does not define any fields");
676: PetscCall(DMGetGlobalSection(dmFine, &gsf));
677: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
678: PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
679: PetscCall(DMSwarmCellDMGetBlockSize(celldm, dmCoarse, &bs));
680: n = Np * bs;
681: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
682: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
683: PetscCall(MatSetType(*mass, dmCoarse->mattype));
684: PetscCall(DMGetApplicationContext(dmFine, &ctx));
686: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
687: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
692: {
693: const char *name = "Mass Matrix Square";
694: MPI_Comm comm;
695: DMSwarmCellDM celldm;
696: PetscDS prob;
697: PetscSection fsection, globalFSection;
698: PetscHSetIJ ht;
699: PetscLayout rLayout, colLayout;
700: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
701: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
702: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
703: PetscScalar *elemMat, *elemMatSq;
704: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0;
705: const char **coordFields;
706: PetscReal **coordVals;
707: PetscInt *bs;
709: PetscFunctionBegin;
710: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
711: PetscCall(DMGetCoordinateDim(dmf, &cdim));
712: PetscCall(DMGetDS(dmf, &prob));
713: PetscCall(PetscDSGetNumFields(prob, &Nf));
714: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
715: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
716: PetscCall(DMGetLocalSection(dmf, &fsection));
717: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
718: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
719: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
721: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
722: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
723: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
725: PetscCall(PetscLayoutCreate(comm, &colLayout));
726: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
727: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
728: PetscCall(PetscLayoutSetUp(colLayout));
729: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
730: PetscCall(PetscLayoutDestroy(&colLayout));
732: PetscCall(PetscLayoutCreate(comm, &rLayout));
733: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
734: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
735: PetscCall(PetscLayoutSetUp(rLayout));
736: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
737: PetscCall(PetscLayoutDestroy(&rLayout));
739: PetscCall(DMPlexGetDepth(dmf, &depth));
740: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
741: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
742: PetscCall(PetscMalloc1(maxAdjSize, &adj));
744: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
745: PetscCall(PetscHSetIJCreate(&ht));
746: /* Count nonzeros
747: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
748: */
749: PetscCall(DMSwarmSortGetAccess(dmc));
750: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
751: PetscInt *cindices;
752: PetscInt numCIndices;
753: #if 0
754: PetscInt adjSize = maxAdjSize, a, j;
755: #endif
757: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
758: maxC = PetscMax(maxC, numCIndices);
759: /* Diagonal block */
760: for (PetscInt i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
761: #if 0
762: /* Off-diagonal blocks */
763: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
764: for (a = 0; a < adjSize; ++a) {
765: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
766: const PetscInt ncell = adj[a];
767: PetscInt *ncindices;
768: PetscInt numNCIndices;
770: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
771: {
772: PetscHashIJKey key;
773: PetscBool missing;
775: for (i = 0; i < numCIndices; ++i) {
776: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
777: if (key.i < 0) continue;
778: for (j = 0; j < numNCIndices; ++j) {
779: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
780: if (key.j < 0) continue;
781: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
782: if (missing) {
783: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
784: else ++onz[key.i - rStart];
785: }
786: }
787: }
788: }
789: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
790: }
791: }
792: #endif
793: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
794: }
795: PetscCall(PetscHSetIJDestroy(&ht));
796: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
797: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
798: PetscCall(PetscFree2(dnz, onz));
799: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
800: /* Fill in values
801: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
802: Start just by producing block diagonal
803: Could loop over adjacent cells
804: Produce neighboring element matrix
805: TODO Determine which columns and rows correspond to shared dual vector
806: Do MatMatMult with rectangular matrices
807: Insert block
808: */
809: for (PetscInt field = 0; field < Nf; ++field) {
810: PetscTabulation Tcoarse;
811: PetscObject obj;
812: PetscInt Nc;
814: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
815: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
816: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
817: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
818: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
819: PetscInt *findices, *cindices;
820: PetscInt numFIndices, numCIndices;
822: /* TODO: Use DMField instead of assuming affine */
823: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
824: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
825: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
826: for (PetscInt p = 0; p < numCIndices; ++p) {
827: PetscReal xr[8];
828: PetscInt off = 0;
830: for (PetscInt i = 0; i < Nfc; ++i) {
831: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[p] * bs[i] + b];
832: }
833: PetscCheck(off == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, cdim);
834: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[p * cdim]);
835: }
836: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
837: /* Get elemMat entries by multiplying by weight */
838: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
839: for (PetscInt i = 0; i < numFIndices; ++i) {
840: for (PetscInt p = 0; p < numCIndices; ++p) {
841: for (PetscInt c = 0; c < Nc; ++c) {
842: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
843: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
844: }
845: }
846: }
847: PetscCall(PetscTabulationDestroy(&Tcoarse));
848: for (PetscInt p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
849: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
850: /* Block diagonal */
851: if (numCIndices) {
852: PetscBLASInt blasn, blask;
853: PetscScalar one = 1.0, zero = 0.0;
855: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
856: PetscCall(PetscBLASIntCast(numFIndices, &blask));
857: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
858: }
859: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
860: /* TODO off-diagonal */
861: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
862: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
863: }
864: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
865: }
866: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
867: PetscCall(PetscFree(adj));
868: PetscCall(DMSwarmSortRestoreAccess(dmc));
869: PetscCall(PetscFree3(v0, J, invJ));
870: PetscCall(PetscFree2(coordVals, bs));
871: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
872: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: /*@
877: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
879: Collective
881: Input Parameters:
882: + dmCoarse - a `DMSWARM`
883: - dmFine - a `DMPLEX`
885: Output Parameter:
886: . mass - the square of the particle mass matrix
888: Level: advanced
890: Note:
891: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
892: future to compute the full normal equations.
894: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
895: @*/
896: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
897: {
898: PetscInt n;
899: void *ctx;
901: PetscFunctionBegin;
902: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
903: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
904: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
905: PetscCall(MatSetType(*mass, dmCoarse->mattype));
906: PetscCall(DMGetApplicationContext(dmFine, &ctx));
908: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
909: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: /*@
914: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
916: Collective
918: Input Parameters:
919: + dm - a `DMSWARM`
920: - fieldname - the textual name given to a registered field
922: Output Parameter:
923: . vec - the vector
925: Level: beginner
927: Note:
928: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
930: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
931: @*/
932: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
933: {
934: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
936: PetscFunctionBegin;
938: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }
942: /*@
943: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
945: Collective
947: Input Parameters:
948: + dm - a `DMSWARM`
949: - fieldname - the textual name given to a registered field
951: Output Parameter:
952: . vec - the vector
954: Level: beginner
956: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
957: @*/
958: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
959: {
960: PetscFunctionBegin;
962: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: /*@
967: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
969: Collective
971: Input Parameters:
972: + dm - a `DMSWARM`
973: - fieldname - the textual name given to a registered field
975: Output Parameter:
976: . vec - the vector
978: Level: beginner
980: Note:
981: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
983: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
984: @*/
985: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
986: {
987: MPI_Comm comm = PETSC_COMM_SELF;
989: PetscFunctionBegin;
990: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: /*@
995: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
997: Collective
999: Input Parameters:
1000: + dm - a `DMSWARM`
1001: - fieldname - the textual name given to a registered field
1003: Output Parameter:
1004: . vec - the vector
1006: Level: beginner
1008: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
1009: @*/
1010: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1011: {
1012: PetscFunctionBegin;
1014: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: /*@
1019: DMSwarmCreateGlobalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1021: Collective
1023: Input Parameters:
1024: + dm - a `DMSWARM`
1025: . Nf - the number of fields
1026: - fieldnames - the textual names given to the registered fields
1028: Output Parameter:
1029: . vec - the vector
1031: Level: beginner
1033: Notes:
1034: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromFields()`.
1036: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateGlobalVectorFromField()`
1038: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromFields()`
1039: @*/
1040: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1041: {
1042: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1044: PetscFunctionBegin;
1046: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1047: PetscFunctionReturn(PETSC_SUCCESS);
1048: }
1050: /*@
1051: DMSwarmDestroyGlobalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1053: Collective
1055: Input Parameters:
1056: + dm - a `DMSWARM`
1057: . Nf - the number of fields
1058: - fieldnames - the textual names given to the registered fields
1060: Output Parameter:
1061: . vec - the vector
1063: Level: beginner
1065: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1066: @*/
1067: PetscErrorCode DMSwarmDestroyGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1068: {
1069: PetscFunctionBegin;
1071: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1072: PetscFunctionReturn(PETSC_SUCCESS);
1073: }
1075: /*@
1076: DMSwarmCreateLocalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1078: Collective
1080: Input Parameters:
1081: + dm - a `DMSWARM`
1082: . Nf - the number of fields
1083: - fieldnames - the textual names given to the registered fields
1085: Output Parameter:
1086: . vec - the vector
1088: Level: beginner
1090: Notes:
1091: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1093: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateLocalVectorFromField()`
1095: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1096: @*/
1097: PetscErrorCode DMSwarmCreateLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1098: {
1099: MPI_Comm comm = PETSC_COMM_SELF;
1101: PetscFunctionBegin;
1102: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1103: PetscFunctionReturn(PETSC_SUCCESS);
1104: }
1106: /*@
1107: DMSwarmDestroyLocalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1109: Collective
1111: Input Parameters:
1112: + dm - a `DMSWARM`
1113: . Nf - the number of fields
1114: - fieldnames - the textual names given to the registered fields
1116: Output Parameter:
1117: . vec - the vector
1119: Level: beginner
1121: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromFields()`
1122: @*/
1123: PetscErrorCode DMSwarmDestroyLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1124: {
1125: PetscFunctionBegin;
1127: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1128: PetscFunctionReturn(PETSC_SUCCESS);
1129: }
1131: /*@
1132: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
1134: Collective
1136: Input Parameter:
1137: . dm - a `DMSWARM`
1139: Level: beginner
1141: Note:
1142: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
1144: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1145: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1146: @*/
1147: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
1148: {
1149: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1151: PetscFunctionBegin;
1152: if (!swarm->field_registration_initialized) {
1153: swarm->field_registration_initialized = PETSC_TRUE;
1154: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
1155: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
1156: }
1157: PetscFunctionReturn(PETSC_SUCCESS);
1158: }
1160: /*@
1161: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
1163: Collective
1165: Input Parameter:
1166: . dm - a `DMSWARM`
1168: Level: beginner
1170: Note:
1171: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
1173: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1174: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1175: @*/
1176: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
1177: {
1178: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1180: PetscFunctionBegin;
1181: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
1182: swarm->field_registration_finalized = PETSC_TRUE;
1183: PetscFunctionReturn(PETSC_SUCCESS);
1184: }
1186: /*@
1187: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
1189: Not Collective
1191: Input Parameters:
1192: + sw - a `DMSWARM`
1193: . nlocal - the length of each registered field
1194: - buffer - the length of the buffer used to efficient dynamic re-sizing
1196: Level: beginner
1198: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
1199: @*/
1200: PetscErrorCode DMSwarmSetLocalSizes(DM sw, PetscInt nlocal, PetscInt buffer)
1201: {
1202: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1203: PetscMPIInt rank;
1204: PetscInt *rankval;
1206: PetscFunctionBegin;
1207: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
1208: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
1209: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
1211: // Initialize values in pid and rank placeholders
1212: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1213: PetscCall(DMSwarmGetField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1214: for (PetscInt p = 0; p < nlocal; p++) rankval[p] = rank;
1215: PetscCall(DMSwarmRestoreField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1216: /* TODO: [pid - use MPI_Scan] */
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: /*@
1221: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
1223: Collective
1225: Input Parameters:
1226: + sw - a `DMSWARM`
1227: - dm - the `DM` to attach to the `DMSWARM`
1229: Level: beginner
1231: Note:
1232: The attached `DM` (dm) will be queried for point location and
1233: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1235: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
1236: @*/
1237: PetscErrorCode DMSwarmSetCellDM(DM sw, DM dm)
1238: {
1239: DMSwarmCellDM celldm;
1240: const char *name;
1241: char *coordName;
1243: PetscFunctionBegin;
1246: PetscCall(PetscStrallocpy(DMSwarmPICField_coor, &coordName));
1247: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 1, (const char **)&coordName, &celldm));
1248: PetscCall(PetscFree(coordName));
1249: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1250: PetscCall(DMSwarmAddCellDM(sw, celldm));
1251: PetscCall(DMSwarmCellDMDestroy(&celldm));
1252: PetscCall(DMSwarmSetCellDMActive(sw, name));
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: /*@
1257: DMSwarmGetCellDM - Fetches the active cell `DM`
1259: Collective
1261: Input Parameter:
1262: . sw - a `DMSWARM`
1264: Output Parameter:
1265: . dm - the active `DM` for the `DMSWARM`
1267: Level: beginner
1269: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1270: @*/
1271: PetscErrorCode DMSwarmGetCellDM(DM sw, DM *dm)
1272: {
1273: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1274: DMSwarmCellDM celldm;
1276: PetscFunctionBegin;
1278: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)&celldm));
1279: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "There is no cell DM named %s in this Swarm", swarm->activeCellDM);
1280: PetscCall(DMSwarmCellDMGetDM(celldm, dm));
1281: PetscFunctionReturn(PETSC_SUCCESS);
1282: }
1284: /*@C
1285: DMSwarmGetCellDMNames - Get the list of cell `DM` names
1287: Not collective
1289: Input Parameter:
1290: . sw - a `DMSWARM`
1292: Output Parameters:
1293: + Ndm - the number of `DMSwarmCellDM` in the `DMSWARM`
1294: - celldms - the name of each `DMSwarmCellDM`
1296: Level: beginner
1298: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMByName()`
1299: @*/
1300: PetscErrorCode DMSwarmGetCellDMNames(DM sw, PetscInt *Ndm, const char **celldms[])
1301: {
1302: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1303: PetscObjectList next = swarm->cellDMs;
1304: PetscInt n = 0;
1306: PetscFunctionBegin;
1308: PetscAssertPointer(Ndm, 2);
1309: PetscAssertPointer(celldms, 3);
1310: while (next) {
1311: next = next->next;
1312: ++n;
1313: }
1314: PetscCall(PetscMalloc1(n, celldms));
1315: next = swarm->cellDMs;
1316: n = 0;
1317: while (next) {
1318: (*celldms)[n] = (const char *)next->obj->name;
1319: next = next->next;
1320: ++n;
1321: }
1322: *Ndm = n;
1323: PetscFunctionReturn(PETSC_SUCCESS);
1324: }
1326: /*@
1327: DMSwarmSetCellDMActive - Activates a cell `DM` for a `DMSWARM`
1329: Collective
1331: Input Parameters:
1332: + sw - a `DMSWARM`
1333: - name - name of the cell `DM` to active for the `DMSWARM`
1335: Level: beginner
1337: Note:
1338: The attached `DM` (dmcell) will be queried for point location and
1339: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1341: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1342: @*/
1343: PetscErrorCode DMSwarmSetCellDMActive(DM sw, const char name[])
1344: {
1345: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1346: DMSwarmCellDM celldm;
1348: PetscFunctionBegin;
1350: PetscCall(PetscInfo(sw, "Setting cell DM to %s\n", name));
1351: PetscCall(PetscFree(swarm->activeCellDM));
1352: PetscCall(PetscStrallocpy(name, (char **)&swarm->activeCellDM));
1353: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: /*@
1358: DMSwarmGetCellDMActive - Returns the active cell `DM` for a `DMSWARM`
1360: Collective
1362: Input Parameter:
1363: . sw - a `DMSWARM`
1365: Output Parameter:
1366: . celldm - the active `DMSwarmCellDM`
1368: Level: beginner
1370: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1371: @*/
1372: PetscErrorCode DMSwarmGetCellDMActive(DM sw, DMSwarmCellDM *celldm)
1373: {
1374: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1376: PetscFunctionBegin;
1378: PetscAssertPointer(celldm, 2);
1379: PetscCheck(swarm->activeCellDM, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no active cell DM");
1380: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)celldm));
1381: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", swarm->activeCellDM);
1382: PetscFunctionReturn(PETSC_SUCCESS);
1383: }
1385: /*@C
1386: DMSwarmGetCellDMByName - Get a `DMSwarmCellDM` from its name
1388: Not collective
1390: Input Parameters:
1391: + sw - a `DMSWARM`
1392: - name - the name
1394: Output Parameter:
1395: . celldm - the `DMSwarmCellDM`
1397: Level: beginner
1399: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMNames()`
1400: @*/
1401: PetscErrorCode DMSwarmGetCellDMByName(DM sw, const char name[], DMSwarmCellDM *celldm)
1402: {
1403: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1405: PetscFunctionBegin;
1407: PetscAssertPointer(name, 2);
1408: PetscAssertPointer(celldm, 3);
1409: PetscCall(PetscObjectListFind(swarm->cellDMs, name, (PetscObject *)celldm));
1410: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", name);
1411: PetscFunctionReturn(PETSC_SUCCESS);
1412: }
1414: /*@
1415: DMSwarmAddCellDM - Adds a cell `DM` to the `DMSWARM`
1417: Collective
1419: Input Parameters:
1420: + sw - a `DMSWARM`
1421: - celldm - the `DMSwarmCellDM`
1423: Level: beginner
1425: Note:
1426: Cell DMs with the same name will share the cellid field
1428: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPushCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1429: @*/
1430: PetscErrorCode DMSwarmAddCellDM(DM sw, DMSwarmCellDM celldm)
1431: {
1432: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1433: const char *name;
1434: PetscInt dim;
1435: PetscBool flg;
1436: MPI_Comm comm;
1438: PetscFunctionBegin;
1440: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1442: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1443: PetscCall(PetscObjectListAdd(&swarm->cellDMs, name, (PetscObject)celldm));
1444: PetscCall(DMGetDimension(sw, &dim));
1445: for (PetscInt f = 0; f < celldm->Nfc; ++f) {
1446: PetscCall(DMSwarmDataFieldStringInList(celldm->coordFields[f], swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1447: if (!flg) {
1448: PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->coordFields[f], dim, PETSC_DOUBLE));
1449: } else {
1450: PetscDataType dt;
1451: PetscInt bs;
1453: PetscCall(DMSwarmGetFieldInfo(sw, celldm->coordFields[f], &bs, &dt));
1454: PetscCheck(bs == dim, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has blocksize %" PetscInt_FMT " != %" PetscInt_FMT " spatial dimension", celldm->coordFields[f], bs, dim);
1455: PetscCheck(dt == PETSC_DOUBLE, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has datatype %s != PETSC_DOUBLE", celldm->coordFields[f], PetscDataTypes[dt]);
1456: }
1457: }
1458: // Assume that DMs with the same name share the cellid field
1459: PetscCall(DMSwarmDataFieldStringInList(celldm->cellid, swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1460: if (!flg) {
1461: PetscBool isShell, isDummy;
1462: const char *name;
1464: // Allow dummy DMSHELL (I don't think we should support this mode)
1465: PetscCall(PetscObjectTypeCompare((PetscObject)celldm->dm, DMSHELL, &isShell));
1466: PetscCall(PetscObjectGetName((PetscObject)celldm->dm, &name));
1467: PetscCall(PetscStrcmp(name, "dummy", &isDummy));
1468: if (!isShell || !isDummy) PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->cellid, 1, PETSC_INT));
1469: }
1470: PetscCall(DMSwarmSetCellDMActive(sw, name));
1471: PetscFunctionReturn(PETSC_SUCCESS);
1472: }
1474: /*@
1475: DMSwarmGetLocalSize - Retrieves the local length of fields registered
1477: Not Collective
1479: Input Parameter:
1480: . dm - a `DMSWARM`
1482: Output Parameter:
1483: . nlocal - the length of each registered field
1485: Level: beginner
1487: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
1488: @*/
1489: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
1490: {
1491: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1493: PetscFunctionBegin;
1494: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
1495: PetscFunctionReturn(PETSC_SUCCESS);
1496: }
1498: /*@
1499: DMSwarmGetSize - Retrieves the total length of fields registered
1501: Collective
1503: Input Parameter:
1504: . dm - a `DMSWARM`
1506: Output Parameter:
1507: . n - the total length of each registered field
1509: Level: beginner
1511: Note:
1512: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
1514: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
1515: @*/
1516: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
1517: {
1518: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1519: PetscInt nlocal;
1521: PetscFunctionBegin;
1522: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1523: PetscCallMPI(MPIU_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1524: PetscFunctionReturn(PETSC_SUCCESS);
1525: }
1527: /*@
1528: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
1530: Collective
1532: Input Parameters:
1533: + dm - a `DMSWARM`
1534: . fieldname - the textual name to identify this field
1535: . blocksize - the number of each data type
1536: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
1538: Level: beginner
1540: Notes:
1541: The textual name for each registered field must be unique.
1543: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1544: @*/
1545: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1546: {
1547: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1548: size_t size;
1550: PetscFunctionBegin;
1551: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1552: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1554: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1555: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1556: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1557: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1558: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1560: PetscCall(PetscDataTypeGetSize(type, &size));
1561: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1562: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1563: {
1564: DMSwarmDataField gfield;
1566: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1567: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1568: }
1569: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1570: PetscFunctionReturn(PETSC_SUCCESS);
1571: }
1573: /*@C
1574: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
1576: Collective
1578: Input Parameters:
1579: + dm - a `DMSWARM`
1580: . fieldname - the textual name to identify this field
1581: - size - the size in bytes of the user struct of each data type
1583: Level: beginner
1585: Note:
1586: The textual name for each registered field must be unique.
1588: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1589: @*/
1590: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1591: {
1592: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1594: PetscFunctionBegin;
1595: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1596: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1597: PetscFunctionReturn(PETSC_SUCCESS);
1598: }
1600: /*@C
1601: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1603: Collective
1605: Input Parameters:
1606: + dm - a `DMSWARM`
1607: . fieldname - the textual name to identify this field
1608: . size - the size in bytes of the user data type
1609: - blocksize - the number of each data type
1611: Level: beginner
1613: Note:
1614: The textual name for each registered field must be unique.
1616: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1617: @*/
1618: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1619: {
1620: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1622: PetscFunctionBegin;
1623: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1624: {
1625: DMSwarmDataField gfield;
1627: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1628: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1629: }
1630: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1631: PetscFunctionReturn(PETSC_SUCCESS);
1632: }
1634: /*@C
1635: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1637: Not Collective, No Fortran Support
1639: Input Parameters:
1640: + dm - a `DMSWARM`
1641: - fieldname - the textual name to identify this field
1643: Output Parameters:
1644: + blocksize - the number of each data type
1645: . type - the data type
1646: - data - pointer to raw array
1648: Level: beginner
1650: Notes:
1651: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1653: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1654: @*/
1655: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1656: {
1657: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1658: DMSwarmDataField gfield;
1660: PetscFunctionBegin;
1662: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1663: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1664: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1665: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1666: if (blocksize) *blocksize = gfield->bs;
1667: if (type) *type = gfield->petsc_type;
1668: PetscFunctionReturn(PETSC_SUCCESS);
1669: }
1671: /*@C
1672: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1674: Not Collective, No Fortran Support
1676: Input Parameters:
1677: + dm - a `DMSWARM`
1678: - fieldname - the textual name to identify this field
1680: Output Parameters:
1681: + blocksize - the number of each data type
1682: . type - the data type
1683: - data - pointer to raw array
1685: Level: beginner
1687: Notes:
1688: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1690: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1691: @*/
1692: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1693: {
1694: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1695: DMSwarmDataField gfield;
1697: PetscFunctionBegin;
1699: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1700: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1701: if (data) *data = NULL;
1702: PetscFunctionReturn(PETSC_SUCCESS);
1703: }
1705: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
1706: {
1707: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1708: DMSwarmDataField gfield;
1710: PetscFunctionBegin;
1712: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1713: if (blocksize) *blocksize = gfield->bs;
1714: if (type) *type = gfield->petsc_type;
1715: PetscFunctionReturn(PETSC_SUCCESS);
1716: }
1718: /*@
1719: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
1721: Not Collective
1723: Input Parameter:
1724: . dm - a `DMSWARM`
1726: Level: beginner
1728: Notes:
1729: The new point will have all fields initialized to zero.
1731: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
1732: @*/
1733: PetscErrorCode DMSwarmAddPoint(DM dm)
1734: {
1735: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1737: PetscFunctionBegin;
1738: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1739: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1740: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
1741: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1742: PetscFunctionReturn(PETSC_SUCCESS);
1743: }
1745: /*@
1746: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
1748: Not Collective
1750: Input Parameters:
1751: + dm - a `DMSWARM`
1752: - npoints - the number of new points to add
1754: Level: beginner
1756: Notes:
1757: The new point will have all fields initialized to zero.
1759: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
1760: @*/
1761: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1762: {
1763: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1764: PetscInt nlocal;
1766: PetscFunctionBegin;
1767: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1768: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1769: nlocal = nlocal + npoints;
1770: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
1771: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1772: PetscFunctionReturn(PETSC_SUCCESS);
1773: }
1775: /*@
1776: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
1778: Not Collective
1780: Input Parameter:
1781: . dm - a `DMSWARM`
1783: Level: beginner
1785: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
1786: @*/
1787: PetscErrorCode DMSwarmRemovePoint(DM dm)
1788: {
1789: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1791: PetscFunctionBegin;
1792: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1793: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
1794: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1795: PetscFunctionReturn(PETSC_SUCCESS);
1796: }
1798: /*@
1799: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
1801: Not Collective
1803: Input Parameters:
1804: + dm - a `DMSWARM`
1805: - idx - index of point to remove
1807: Level: beginner
1809: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1810: @*/
1811: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
1812: {
1813: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1815: PetscFunctionBegin;
1816: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1817: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
1818: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1819: PetscFunctionReturn(PETSC_SUCCESS);
1820: }
1822: /*@
1823: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
1825: Not Collective
1827: Input Parameters:
1828: + dm - a `DMSWARM`
1829: . pi - the index of the point to copy
1830: - pj - the point index where the copy should be located
1832: Level: beginner
1834: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1835: @*/
1836: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
1837: {
1838: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1840: PetscFunctionBegin;
1841: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1842: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
1843: PetscFunctionReturn(PETSC_SUCCESS);
1844: }
1846: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
1847: {
1848: PetscFunctionBegin;
1849: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
1850: PetscFunctionReturn(PETSC_SUCCESS);
1851: }
1853: /*@
1854: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
1856: Collective
1858: Input Parameters:
1859: + dm - the `DMSWARM`
1860: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1862: Level: advanced
1864: Notes:
1865: The `DM` will be modified to accommodate received points.
1866: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
1867: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
1869: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
1870: @*/
1871: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
1872: {
1873: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1875: PetscFunctionBegin;
1876: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
1877: switch (swarm->migrate_type) {
1878: case DMSWARM_MIGRATE_BASIC:
1879: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
1880: break;
1881: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1882: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
1883: break;
1884: case DMSWARM_MIGRATE_DMCELLEXACT:
1885: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1886: case DMSWARM_MIGRATE_USER:
1887: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
1888: default:
1889: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
1890: }
1891: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
1892: PetscCall(DMClearGlobalVectors(dm));
1893: PetscFunctionReturn(PETSC_SUCCESS);
1894: }
1896: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
1898: /*
1899: DMSwarmCollectViewCreate
1901: * Applies a collection method and gathers point neighbour points into dm
1903: Notes:
1904: Users should call DMSwarmCollectViewDestroy() after
1905: they have finished computations associated with the collected points
1906: */
1908: /*@
1909: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1910: in neighbour ranks into the `DMSWARM`
1912: Collective
1914: Input Parameter:
1915: . dm - the `DMSWARM`
1917: Level: advanced
1919: Notes:
1920: Users should call `DMSwarmCollectViewDestroy()` after
1921: they have finished computations associated with the collected points
1923: Different collect methods are supported. See `DMSwarmSetCollectType()`.
1925: Developer Note:
1926: Create and Destroy routines create new objects that can get destroyed, they do not change the state
1927: of the current object.
1929: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
1930: @*/
1931: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1932: {
1933: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1934: PetscInt ng;
1936: PetscFunctionBegin;
1937: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
1938: PetscCall(DMSwarmGetLocalSize(dm, &ng));
1939: switch (swarm->collect_type) {
1940: case DMSWARM_COLLECT_BASIC:
1941: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
1942: break;
1943: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1944: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1945: case DMSWARM_COLLECT_GENERAL:
1946: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
1947: default:
1948: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
1949: }
1950: swarm->collect_view_active = PETSC_TRUE;
1951: swarm->collect_view_reset_nlocal = ng;
1952: PetscFunctionReturn(PETSC_SUCCESS);
1953: }
1955: /*@
1956: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
1958: Collective
1960: Input Parameters:
1961: . dm - the `DMSWARM`
1963: Notes:
1964: Users should call `DMSwarmCollectViewCreate()` before this function is called.
1966: Level: advanced
1968: Developer Note:
1969: Create and Destroy routines create new objects that can get destroyed, they do not change the state
1970: of the current object.
1972: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
1973: @*/
1974: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1975: {
1976: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1978: PetscFunctionBegin;
1979: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
1980: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
1981: swarm->collect_view_active = PETSC_FALSE;
1982: PetscFunctionReturn(PETSC_SUCCESS);
1983: }
1985: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
1986: {
1987: PetscInt dim;
1989: PetscFunctionBegin;
1990: PetscCall(DMSwarmSetNumSpecies(dm, 1));
1991: PetscCall(DMGetDimension(dm, &dim));
1992: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1993: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1994: PetscFunctionReturn(PETSC_SUCCESS);
1995: }
1997: /*@
1998: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
2000: Collective
2002: Input Parameters:
2003: + dm - the `DMSWARM`
2004: - Npc - The number of particles per cell in the cell `DM`
2006: Level: intermediate
2008: Notes:
2009: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
2010: one particle is in each cell, it is placed at the centroid.
2012: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
2013: @*/
2014: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
2015: {
2016: DM cdm;
2017: DMSwarmCellDM celldm;
2018: PetscRandom rnd;
2019: DMPolytopeType ct;
2020: PetscBool simplex;
2021: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
2022: PetscInt dim, d, cStart, cEnd, c, p, Nfc;
2023: const char **coordFields;
2025: PetscFunctionBeginUser;
2026: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
2027: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
2028: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
2030: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2031: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2032: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2033: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2034: PetscCall(DMGetDimension(cdm, &dim));
2035: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
2036: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
2037: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
2039: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
2040: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
2041: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2042: for (c = cStart; c < cEnd; ++c) {
2043: if (Npc == 1) {
2044: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
2045: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
2046: } else {
2047: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
2048: for (p = 0; p < Npc; ++p) {
2049: const PetscInt n = c * Npc + p;
2050: PetscReal sum = 0.0, refcoords[3];
2052: for (d = 0; d < dim; ++d) {
2053: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
2054: sum += refcoords[d];
2055: }
2056: if (simplex && sum > 0.0)
2057: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
2058: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
2059: }
2060: }
2061: }
2062: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2063: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
2064: PetscCall(PetscRandomDestroy(&rnd));
2065: PetscFunctionReturn(PETSC_SUCCESS);
2066: }
2068: /*@
2069: DMSwarmGetType - Get particular flavor of `DMSWARM`
2071: Collective
2073: Input Parameter:
2074: . sw - the `DMSWARM`
2076: Output Parameter:
2077: . stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2079: Level: advanced
2081: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2082: @*/
2083: PetscErrorCode DMSwarmGetType(DM sw, DMSwarmType *stype)
2084: {
2085: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2087: PetscFunctionBegin;
2089: PetscAssertPointer(stype, 2);
2090: *stype = swarm->swarm_type;
2091: PetscFunctionReturn(PETSC_SUCCESS);
2092: }
2094: /*@
2095: DMSwarmSetType - Set particular flavor of `DMSWARM`
2097: Collective
2099: Input Parameters:
2100: + sw - the `DMSWARM`
2101: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2103: Level: advanced
2105: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2106: @*/
2107: PetscErrorCode DMSwarmSetType(DM sw, DMSwarmType stype)
2108: {
2109: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2111: PetscFunctionBegin;
2113: swarm->swarm_type = stype;
2114: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(sw));
2115: PetscFunctionReturn(PETSC_SUCCESS);
2116: }
2118: static PetscErrorCode DMSwarmCreateRemapDM_Private(DM sw, DM *rdm)
2119: {
2120: PetscFE fe;
2121: DMPolytopeType ct;
2122: PetscInt dim, cStart;
2123: const char *prefix = "remap_";
2125: PetscFunctionBegin;
2126: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), rdm));
2127: PetscCall(DMSetType(*rdm, DMPLEX));
2128: PetscCall(DMPlexSetOptionsPrefix(*rdm, prefix));
2129: PetscCall(DMSetFromOptions(*rdm));
2130: PetscCall(PetscObjectSetName((PetscObject)*rdm, "remap"));
2131: PetscCall(DMViewFromOptions(*rdm, NULL, "-dm_view"));
2133: PetscCall(DMGetDimension(*rdm, &dim));
2134: PetscCall(DMPlexGetHeightStratum(*rdm, 0, &cStart, NULL));
2135: PetscCall(DMPlexGetCellType(*rdm, cStart, &ct));
2136: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
2137: PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
2138: PetscCall(DMSetField(*rdm, 0, NULL, (PetscObject)fe));
2139: PetscCall(DMCreateDS(*rdm));
2140: PetscCall(PetscFEDestroy(&fe));
2141: PetscFunctionReturn(PETSC_SUCCESS);
2142: }
2144: static PetscErrorCode DMSetup_Swarm(DM sw)
2145: {
2146: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2148: PetscFunctionBegin;
2149: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
2150: swarm->issetup = PETSC_TRUE;
2152: if (swarm->remap_type != DMSWARM_REMAP_NONE) {
2153: DMSwarmCellDM celldm;
2154: DM rdm;
2155: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
2156: const char *vfieldnames[1] = {"w_q"};
2158: PetscCall(DMSwarmCreateRemapDM_Private(sw, &rdm));
2159: PetscCall(DMSwarmCellDMCreate(rdm, 1, vfieldnames, 2, fieldnames, &celldm));
2160: PetscCall(DMSwarmAddCellDM(sw, celldm));
2161: PetscCall(DMSwarmCellDMDestroy(&celldm));
2162: PetscCall(DMDestroy(&rdm));
2163: }
2165: if (swarm->swarm_type == DMSWARM_PIC) {
2166: DMSwarmCellDM celldm;
2168: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2169: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "No active cell DM. DMSWARM_PIC requires you call DMSwarmSetCellDM() or DMSwarmAddCellDM()");
2170: if (celldm->dm->ops->locatepointssubdomain) {
2171: /* check methods exists for exact ownership identificiation */
2172: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
2173: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
2174: } else {
2175: /* check methods exist for point location AND rank neighbor identification */
2176: if (celldm->dm->ops->locatepoints) {
2177: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
2178: } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
2180: if (celldm->dm->ops->getneighbors) {
2181: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
2182: } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
2184: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
2185: }
2186: }
2188: PetscCall(DMSwarmFinalizeFieldRegister(sw));
2190: /* check some fields were registered */
2191: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
2192: PetscFunctionReturn(PETSC_SUCCESS);
2193: }
2195: static PetscErrorCode DMDestroy_Swarm(DM dm)
2196: {
2197: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2199: PetscFunctionBegin;
2200: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2201: PetscCall(PetscObjectListDestroy(&swarm->cellDMs));
2202: PetscCall(PetscFree(swarm->activeCellDM));
2203: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
2204: PetscCall(PetscFree(swarm));
2205: PetscFunctionReturn(PETSC_SUCCESS);
2206: }
2208: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
2209: {
2210: DM cdm;
2211: DMSwarmCellDM celldm;
2212: PetscDraw draw;
2213: PetscReal *coords, oldPause, radius = 0.01;
2214: PetscInt Np, p, bs, Nfc;
2215: const char **coordFields;
2217: PetscFunctionBegin;
2218: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
2219: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2220: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2221: PetscCall(PetscDrawGetPause(draw, &oldPause));
2222: PetscCall(PetscDrawSetPause(draw, 0.0));
2223: PetscCall(DMView(cdm, viewer));
2224: PetscCall(PetscDrawSetPause(draw, oldPause));
2226: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2227: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2228: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2229: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2230: PetscCall(DMSwarmGetField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2231: for (p = 0; p < Np; ++p) {
2232: const PetscInt i = p * bs;
2234: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
2235: }
2236: PetscCall(DMSwarmRestoreField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2237: PetscCall(PetscDrawFlush(draw));
2238: PetscCall(PetscDrawPause(draw));
2239: PetscCall(PetscDrawSave(draw));
2240: PetscFunctionReturn(PETSC_SUCCESS);
2241: }
2243: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
2244: {
2245: PetscViewerFormat format;
2246: PetscInt *sizes;
2247: PetscInt dim, Np, maxSize = 17;
2248: MPI_Comm comm;
2249: PetscMPIInt rank, size, p;
2250: const char *name, *cellid;
2252: PetscFunctionBegin;
2253: PetscCall(PetscViewerGetFormat(viewer, &format));
2254: PetscCall(DMGetDimension(dm, &dim));
2255: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2256: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2257: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2258: PetscCallMPI(MPI_Comm_size(comm, &size));
2259: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
2260: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
2261: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
2262: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
2263: else PetscCall(PetscCalloc1(3, &sizes));
2264: if (size < maxSize) {
2265: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
2266: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
2267: for (p = 0; p < size; ++p) {
2268: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
2269: }
2270: } else {
2271: PetscInt locMinMax[2] = {Np, Np};
2273: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
2274: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
2275: }
2276: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2277: PetscCall(PetscFree(sizes));
2278: if (format == PETSC_VIEWER_ASCII_INFO) {
2279: DMSwarmCellDM celldm;
2280: PetscInt *cell;
2282: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
2283: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2284: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2285: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
2286: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&cell));
2287: for (p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%d: %" PetscInt_FMT "\n", p, cell[p]));
2288: PetscCall(PetscViewerFlush(viewer));
2289: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2290: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&cell));
2291: }
2292: PetscFunctionReturn(PETSC_SUCCESS);
2293: }
2295: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
2296: {
2297: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2298: PetscBool iascii, ibinary, isvtk, isdraw, ispython;
2299: #if defined(PETSC_HAVE_HDF5)
2300: PetscBool ishdf5;
2301: #endif
2303: PetscFunctionBegin;
2306: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2307: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
2308: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
2309: #if defined(PETSC_HAVE_HDF5)
2310: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2311: #endif
2312: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2313: PetscCall(PetscObjectHasFunction((PetscObject)viewer, "PetscViewerPythonViewObject_C", &ispython));
2314: if (iascii) {
2315: PetscViewerFormat format;
2317: PetscCall(PetscViewerGetFormat(viewer, &format));
2318: switch (format) {
2319: case PETSC_VIEWER_ASCII_INFO_DETAIL:
2320: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
2321: break;
2322: default:
2323: PetscCall(DMView_Swarm_Ascii(dm, viewer));
2324: }
2325: } else {
2326: #if defined(PETSC_HAVE_HDF5)
2327: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
2328: #endif
2329: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
2330: if (ispython) PetscCall(PetscViewerPythonViewObject(viewer, (PetscObject)dm));
2331: }
2332: PetscFunctionReturn(PETSC_SUCCESS);
2333: }
2335: /*@
2336: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
2337: The cell `DM` is filtered for fields of that cell, and the filtered `DM` is used as the cell `DM` of the new swarm object.
2339: Noncollective
2341: Input Parameters:
2342: + sw - the `DMSWARM`
2343: . cellID - the integer id of the cell to be extracted and filtered
2344: - cellswarm - The `DMSWARM` to receive the cell
2346: Level: beginner
2348: Notes:
2349: This presently only supports `DMSWARM_PIC` type
2351: Should be restored with `DMSwarmRestoreCellSwarm()`
2353: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
2355: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
2356: @*/
2357: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2358: {
2359: DM_Swarm *original = (DM_Swarm *)sw->data;
2360: DMLabel label;
2361: DM dmc, subdmc;
2362: PetscInt *pids, particles, dim;
2363: const char *name;
2365: PetscFunctionBegin;
2366: /* Configure new swarm */
2367: PetscCall(DMSetType(cellswarm, DMSWARM));
2368: PetscCall(DMGetDimension(sw, &dim));
2369: PetscCall(DMSetDimension(cellswarm, dim));
2370: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
2371: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
2372: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
2373: PetscCall(DMSwarmSortGetAccess(sw));
2374: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
2375: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2376: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
2377: PetscCall(DMSwarmSortRestoreAccess(sw));
2378: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2379: PetscCall(DMSwarmGetCellDM(sw, &dmc));
2380: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
2381: PetscCall(DMAddLabel(dmc, label));
2382: PetscCall(DMLabelSetValue(label, cellID, 1));
2383: PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdmc));
2384: PetscCall(PetscObjectGetName((PetscObject)dmc, &name));
2385: PetscCall(PetscObjectSetName((PetscObject)subdmc, name));
2386: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
2387: PetscCall(DMLabelDestroy(&label));
2388: PetscFunctionReturn(PETSC_SUCCESS);
2389: }
2391: /*@
2392: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
2394: Noncollective
2396: Input Parameters:
2397: + sw - the parent `DMSWARM`
2398: . cellID - the integer id of the cell to be copied back into the parent swarm
2399: - cellswarm - the cell swarm object
2401: Level: beginner
2403: Note:
2404: This only supports `DMSWARM_PIC` types of `DMSWARM`s
2406: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
2407: @*/
2408: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2409: {
2410: DM dmc;
2411: PetscInt *pids, particles, p;
2413: PetscFunctionBegin;
2414: PetscCall(DMSwarmSortGetAccess(sw));
2415: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2416: PetscCall(DMSwarmSortRestoreAccess(sw));
2417: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
2418: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
2419: /* Free memory, destroy cell dm */
2420: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
2421: PetscCall(DMDestroy(&dmc));
2422: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2423: PetscFunctionReturn(PETSC_SUCCESS);
2424: }
2426: /*@
2427: DMSwarmComputeMoments - Compute the first three particle moments for a given field
2429: Noncollective
2431: Input Parameters:
2432: + sw - the `DMSWARM`
2433: . coordinate - the coordinate field name
2434: - weight - the weight field name
2436: Output Parameter:
2437: . moments - the field moments
2439: Level: intermediate
2441: Notes:
2442: The `moments` array should be of length bs + 2, where bs is the block size of the coordinate field.
2444: The weight field must be a scalar, having blocksize 1.
2446: .seealso: `DM`, `DMSWARM`, `DMPlexComputeMoments()`
2447: @*/
2448: PetscErrorCode DMSwarmComputeMoments(DM sw, const char coordinate[], const char weight[], PetscReal moments[])
2449: {
2450: const PetscReal *coords;
2451: const PetscReal *w;
2452: PetscReal *mom;
2453: PetscDataType dtc, dtw;
2454: PetscInt bsc, bsw, Np;
2455: MPI_Comm comm;
2457: PetscFunctionBegin;
2459: PetscAssertPointer(coordinate, 2);
2460: PetscAssertPointer(weight, 3);
2461: PetscAssertPointer(moments, 4);
2462: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
2463: PetscCall(DMSwarmGetField(sw, coordinate, &bsc, &dtc, (void **)&coords));
2464: PetscCall(DMSwarmGetField(sw, weight, &bsw, &dtw, (void **)&w));
2465: PetscCheck(dtc == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s must be real, not %s", coordinate, PetscDataTypes[dtc]);
2466: PetscCheck(dtw == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be real, not %s", weight, PetscDataTypes[dtw]);
2467: PetscCheck(bsw == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be a scalar, not blocksize %" PetscInt_FMT, weight, bsw);
2468: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2469: PetscCall(DMGetWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2470: PetscCall(PetscArrayzero(mom, bsc + 2));
2471: for (PetscInt p = 0; p < Np; ++p) {
2472: const PetscReal *c = &coords[p * bsc];
2473: const PetscReal wp = w[p];
2475: mom[0] += wp;
2476: for (PetscInt d = 0; d < bsc; ++d) {
2477: mom[d + 1] += wp * c[d];
2478: mom[d + bsc + 1] += wp * PetscSqr(c[d]);
2479: }
2480: }
2481: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&coords));
2482: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
2483: PetscCallMPI(MPIU_Allreduce(mom, moments, bsc + 2, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
2484: PetscCall(DMRestoreWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2485: PetscFunctionReturn(PETSC_SUCCESS);
2486: }
2488: static PetscErrorCode DMSetFromOptions_Swarm(DM dm, PetscOptionItems *PetscOptionsObject)
2489: {
2490: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2492: PetscFunctionBegin;
2493: PetscOptionsHeadBegin(PetscOptionsObject, "DMSwarm Options");
2494: PetscCall(PetscOptionsEnum("-dm_swarm_remap_type", "Remap algorithm", "DMSwarmSetRemapType", DMSwarmRemapTypeNames, (PetscEnum)swarm->remap_type, (PetscEnum *)&swarm->remap_type, NULL));
2495: PetscOptionsHeadEnd();
2496: PetscFunctionReturn(PETSC_SUCCESS);
2497: }
2499: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
2501: static PetscErrorCode DMInitialize_Swarm(DM sw)
2502: {
2503: PetscFunctionBegin;
2504: sw->ops->view = DMView_Swarm;
2505: sw->ops->load = NULL;
2506: sw->ops->setfromoptions = DMSetFromOptions_Swarm;
2507: sw->ops->clone = DMClone_Swarm;
2508: sw->ops->setup = DMSetup_Swarm;
2509: sw->ops->createlocalsection = NULL;
2510: sw->ops->createsectionpermutation = NULL;
2511: sw->ops->createdefaultconstraints = NULL;
2512: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
2513: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
2514: sw->ops->getlocaltoglobalmapping = NULL;
2515: sw->ops->createfieldis = NULL;
2516: sw->ops->createcoordinatedm = NULL;
2517: sw->ops->getcoloring = NULL;
2518: sw->ops->creatematrix = DMCreateMatrix_Swarm;
2519: sw->ops->createinterpolation = NULL;
2520: sw->ops->createinjection = NULL;
2521: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
2522: sw->ops->refine = NULL;
2523: sw->ops->coarsen = NULL;
2524: sw->ops->refinehierarchy = NULL;
2525: sw->ops->coarsenhierarchy = NULL;
2526: sw->ops->globaltolocalbegin = NULL;
2527: sw->ops->globaltolocalend = NULL;
2528: sw->ops->localtoglobalbegin = NULL;
2529: sw->ops->localtoglobalend = NULL;
2530: sw->ops->destroy = DMDestroy_Swarm;
2531: sw->ops->createsubdm = NULL;
2532: sw->ops->getdimpoints = NULL;
2533: sw->ops->locatepoints = NULL;
2534: PetscFunctionReturn(PETSC_SUCCESS);
2535: }
2537: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
2538: {
2539: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2541: PetscFunctionBegin;
2542: swarm->refct++;
2543: (*newdm)->data = swarm;
2544: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
2545: PetscCall(DMInitialize_Swarm(*newdm));
2546: (*newdm)->dim = dm->dim;
2547: PetscFunctionReturn(PETSC_SUCCESS);
2548: }
2550: /*MC
2551: DMSWARM = "swarm" - A `DM` object used to represent arrays of data (fields) of arbitrary data type.
2552: This implementation was designed for particle methods in which the underlying
2553: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
2555: Level: intermediate
2557: Notes:
2558: User data can be represented by `DMSWARM` through a registering "fields".
2559: To register a field, the user must provide;
2560: (a) a unique name;
2561: (b) the data type (or size in bytes);
2562: (c) the block size of the data.
2564: For example, suppose the application requires a unique id, energy, momentum and density to be stored
2565: on a set of particles. Then the following code could be used
2566: .vb
2567: DMSwarmInitializeFieldRegister(dm)
2568: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
2569: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
2570: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
2571: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
2572: DMSwarmFinalizeFieldRegister(dm)
2573: .ve
2575: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
2576: The only restriction imposed by `DMSWARM` is that all fields contain the same number of points.
2578: To support particle methods, "migration" techniques are provided. These methods migrate data
2579: between ranks.
2581: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
2582: As a `DMSWARM` may internally define and store values of different data types,
2583: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
2584: fields should be used to define a `Vec` object via
2585: `DMSwarmVectorDefineField()`
2586: The specified field can be changed at any time - thereby permitting vectors
2587: compatible with different fields to be created.
2589: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via
2590: `DMSwarmCreateGlobalVectorFromField()`
2591: Here the data defining the field in the `DMSWARM` is shared with a Vec.
2592: This is inherently unsafe if you alter the size of the field at any time between
2593: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
2594: If the local size of the `DMSWARM` does not match the local size of the global vector
2595: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
2597: Additional high-level support is provided for Particle-In-Cell methods.
2598: Please refer to `DMSwarmSetType()`.
2600: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`
2601: M*/
2603: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
2604: {
2605: DM_Swarm *swarm;
2607: PetscFunctionBegin;
2609: PetscCall(PetscNew(&swarm));
2610: dm->data = swarm;
2611: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
2612: PetscCall(DMSwarmInitializeFieldRegister(dm));
2613: dm->dim = 0;
2614: swarm->refct = 1;
2615: swarm->issetup = PETSC_FALSE;
2616: swarm->swarm_type = DMSWARM_BASIC;
2617: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
2618: swarm->collect_type = DMSWARM_COLLECT_BASIC;
2619: swarm->migrate_error_on_missing_point = PETSC_FALSE;
2620: swarm->collect_view_active = PETSC_FALSE;
2621: swarm->collect_view_reset_nlocal = -1;
2622: PetscCall(DMInitialize_Swarm(dm));
2623: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
2624: PetscFunctionReturn(PETSC_SUCCESS);
2625: }
2627: /* Replace dm with the contents of ndm, and then destroy ndm
2628: - Share the DM_Swarm structure
2629: */
2630: PetscErrorCode DMSwarmReplace(DM dm, DM *ndm)
2631: {
2632: DM dmNew = *ndm;
2633: const PetscReal *maxCell, *Lstart, *L;
2634: PetscInt dim;
2636: PetscFunctionBegin;
2637: if (dm == dmNew) {
2638: PetscCall(DMDestroy(ndm));
2639: PetscFunctionReturn(PETSC_SUCCESS);
2640: }
2641: dm->setupcalled = dmNew->setupcalled;
2642: if (!dm->hdr.name) {
2643: const char *name;
2645: PetscCall(PetscObjectGetName((PetscObject)*ndm, &name));
2646: PetscCall(PetscObjectSetName((PetscObject)dm, name));
2647: }
2648: PetscCall(DMGetDimension(dmNew, &dim));
2649: PetscCall(DMSetDimension(dm, dim));
2650: PetscCall(DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L));
2651: PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
2652: PetscCall(DMDestroy_Swarm(dm));
2653: PetscCall(DMInitialize_Swarm(dm));
2654: dm->data = dmNew->data;
2655: ((DM_Swarm *)dmNew->data)->refct++;
2656: PetscCall(DMDestroy(ndm));
2657: PetscFunctionReturn(PETSC_SUCCESS);
2658: }
2660: /*@
2661: DMSwarmDuplicate - Creates a new `DMSWARM` with the same fields and cell `DM`s but no particles
2663: Collective
2665: Input Parameter:
2666: . sw - the `DMSWARM`
2668: Output Parameter:
2669: . nsw - the new `DMSWARM`
2671: Level: beginner
2673: .seealso: `DM`, `DMSWARM`, `DMSwarmCreate()`, `DMClone()`
2674: @*/
2675: PetscErrorCode DMSwarmDuplicate(DM sw, DM *nsw)
2676: {
2677: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2678: DMSwarmDataField *fields;
2679: DMSwarmCellDM celldm, ncelldm;
2680: DMSwarmType stype;
2681: const char *name, **celldmnames;
2682: void *ctx;
2683: PetscInt dim, Nf, Ndm;
2684: PetscBool flg;
2686: PetscFunctionBegin;
2687: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), nsw));
2688: PetscCall(DMSetType(*nsw, DMSWARM));
2689: PetscCall(PetscObjectGetName((PetscObject)sw, &name));
2690: PetscCall(PetscObjectSetName((PetscObject)*nsw, name));
2691: PetscCall(DMGetDimension(sw, &dim));
2692: PetscCall(DMSetDimension(*nsw, dim));
2693: PetscCall(DMSwarmGetType(sw, &stype));
2694: PetscCall(DMSwarmSetType(*nsw, stype));
2695: PetscCall(DMGetApplicationContext(sw, &ctx));
2696: PetscCall(DMSetApplicationContext(*nsw, ctx));
2698: PetscCall(DMSwarmDataBucketGetDMSwarmDataFields(swarm->db, &Nf, &fields));
2699: for (PetscInt f = 0; f < Nf; ++f) {
2700: PetscCall(DMSwarmDataFieldStringInList(fields[f]->name, ((DM_Swarm *)(*nsw)->data)->db->nfields, (const DMSwarmDataField *)((DM_Swarm *)(*nsw)->data)->db->field, &flg));
2701: if (!flg) PetscCall(DMSwarmRegisterPetscDatatypeField(*nsw, fields[f]->name, fields[f]->bs, fields[f]->petsc_type));
2702: }
2704: PetscCall(DMSwarmGetCellDMNames(sw, &Ndm, &celldmnames));
2705: for (PetscInt c = 0; c < Ndm; ++c) {
2706: DM dm;
2707: PetscInt Ncf;
2708: const char **coordfields, **fields;
2710: PetscCall(DMSwarmGetCellDMByName(sw, celldmnames[c], &celldm));
2711: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
2712: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncf, &coordfields));
2713: PetscCall(DMSwarmCellDMGetFields(celldm, &Nf, &fields));
2714: PetscCall(DMSwarmCellDMCreate(dm, Nf, fields, Ncf, coordfields, &ncelldm));
2715: PetscCall(DMSwarmAddCellDM(*nsw, ncelldm));
2716: PetscCall(DMSwarmCellDMDestroy(&ncelldm));
2717: }
2718: PetscCall(PetscFree(celldmnames));
2720: PetscCall(DMSetFromOptions(*nsw));
2721: PetscCall(DMSetUp(*nsw));
2722: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2723: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
2724: PetscCall(DMSwarmSetCellDMActive(*nsw, name));
2725: PetscFunctionReturn(PETSC_SUCCESS);
2726: }