Actual source code: swarm.c
1: #include "petscsys.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 *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
23: const char DMSwarmField_pid[] = "DMSwarm_pid";
24: const char DMSwarmField_rank[] = "DMSwarm_rank";
25: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
26: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
28: PetscInt SwarmDataFieldId = -1;
30: #if defined(PETSC_HAVE_HDF5)
31: #include <petscviewerhdf5.h>
33: static PetscErrorCode DMInitialize_Swarm(DM);
34: static PetscErrorCode DMDestroy_Swarm(DM);
36: /* Replace dm with the contents of ndm, and then destroy ndm
37: - Share the DM_Swarm structure
38: */
39: PetscErrorCode DMSwarmReplace_Internal(DM dm, DM *ndm)
40: {
41: DM dmNew = *ndm;
42: const PetscReal *maxCell, *Lstart, *L;
43: PetscInt dim;
45: PetscFunctionBegin;
46: if (dm == dmNew) {
47: PetscCall(DMDestroy(ndm));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
50: dm->setupcalled = dmNew->setupcalled;
51: if (!dm->hdr.name) {
52: const char *name;
54: PetscCall(PetscObjectGetName((PetscObject)*ndm, &name));
55: PetscCall(PetscObjectSetName((PetscObject)dm, name));
56: }
57: PetscCall(DMGetDimension(dmNew, &dim));
58: PetscCall(DMSetDimension(dm, dim));
59: PetscCall(DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L));
60: PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
61: PetscCall(DMDestroy_Swarm(dm));
62: PetscCall(DMInitialize_Swarm(dm));
63: dm->data = dmNew->data;
64: ((DM_Swarm *)dmNew->data)->refct++;
65: PetscCall(DMDestroy(ndm));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
70: {
71: DM dm;
72: PetscReal seqval;
73: PetscInt seqnum, bs;
74: PetscBool isseq, ists;
76: PetscFunctionBegin;
77: PetscCall(VecGetDM(v, &dm));
78: PetscCall(VecGetBlockSize(v, &bs));
79: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particle_fields"));
80: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
81: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
82: PetscCall(PetscViewerHDF5IsTimestepping(viewer, &ists));
83: if (ists) PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
84: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
85: PetscCall(VecViewNative(v, viewer));
86: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
87: PetscCall(PetscViewerHDF5PopGroup(viewer));
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
92: {
93: Vec coordinates;
94: PetscInt Np;
95: PetscBool isseq;
96: const char *coordname;
98: PetscFunctionBegin;
99: PetscCall(DMSwarmGetCoordinateField(dm, &coordname));
100: PetscCall(DMSwarmGetSize(dm, &Np));
101: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, coordname, &coordinates));
102: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
103: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
104: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
105: PetscCall(VecViewNative(coordinates, viewer));
106: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
107: PetscCall(PetscViewerHDF5PopGroup(viewer));
108: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, coordname, &coordinates));
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
111: #endif
113: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
114: {
115: DM dm;
116: #if defined(PETSC_HAVE_HDF5)
117: PetscBool ishdf5;
118: #endif
120: PetscFunctionBegin;
121: PetscCall(VecGetDM(v, &dm));
122: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
123: #if defined(PETSC_HAVE_HDF5)
124: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
125: if (ishdf5) {
126: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
129: #endif
130: PetscCall(VecViewNative(v, viewer));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /*@C
135: DMSwarmVectorGetField - Gets the fields from which to define a `Vec` object
136: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
138: Not collective
140: Input Parameter:
141: . dm - a `DMSWARM`
143: Output Parameters:
144: + Nf - the number of fields
145: - fieldnames - the textual name given to each registered field, or NULL if it has not been set
147: Level: beginner
149: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
150: @*/
151: PetscErrorCode DMSwarmVectorGetField(DM dm, PetscInt *Nf, const char **fieldnames[])
152: {
153: PetscFunctionBegin;
155: if (Nf) {
156: PetscAssertPointer(Nf, 2);
157: *Nf = ((DM_Swarm *)dm->data)->vec_field_num;
158: }
159: if (fieldnames) {
160: PetscAssertPointer(fieldnames, 3);
161: *fieldnames = ((DM_Swarm *)dm->data)->vec_field_names;
162: }
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: /*@
167: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
168: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
170: Collective
172: Input Parameters:
173: + dm - a `DMSWARM`
174: - fieldname - the textual name given to each registered field
176: Level: beginner
178: Notes:
179: The field with name `fieldname` must be defined as having a data type of `PetscScalar`.
181: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
182: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
184: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineFields()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
185: @*/
186: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
187: {
188: PetscFunctionBegin;
189: PetscCall(DMSwarmVectorDefineFields(dm, 1, &fieldname));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: /*@C
194: DMSwarmVectorDefineFields - Sets the fields from which to define a `Vec` object
195: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
197: Collective, No Fortran support
199: Input Parameters:
200: + dm - a `DMSWARM`
201: . Nf - the number of fields
202: - fieldnames - the textual name given to each registered field
204: Level: beginner
206: Notes:
207: Each field with name in `fieldnames` must be defined as having a data type of `PetscScalar`.
209: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
210: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
212: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
213: @*/
214: PetscErrorCode DMSwarmVectorDefineFields(DM dm, PetscInt Nf, const char *fieldnames[])
215: {
216: DM_Swarm *sw = (DM_Swarm *)dm->data;
218: PetscFunctionBegin;
220: if (fieldnames) PetscAssertPointer(fieldnames, 3);
221: if (!sw->issetup) PetscCall(DMSetUp(dm));
222: PetscCheck(Nf >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields must be non-negative, not %" PetscInt_FMT, Nf);
223: for (PetscInt f = 0; f < sw->vec_field_num; ++f) PetscCall(PetscFree(sw->vec_field_names[f]));
224: PetscCall(PetscFree(sw->vec_field_names));
226: sw->vec_field_num = Nf;
227: sw->vec_field_bs = 0;
228: PetscCall(DMSwarmDataBucketGetSizes(sw->db, &sw->vec_field_nlocal, NULL, NULL));
229: PetscCall(PetscMalloc1(Nf, &sw->vec_field_names));
230: for (PetscInt f = 0; f < Nf; ++f) {
231: PetscInt bs;
232: PetscScalar *array;
233: PetscDataType type;
235: PetscCall(DMSwarmGetField(dm, fieldnames[f], &bs, &type, (void **)&array));
236: // Check all fields are of type PETSC_REAL or PETSC_SCALAR
237: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
238: sw->vec_field_bs += bs;
239: PetscCall(DMSwarmRestoreField(dm, fieldnames[f], &bs, &type, (void **)&array));
240: PetscCall(PetscStrallocpy(fieldnames[f], (char **)&sw->vec_field_names[f]));
241: }
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /* requires DMSwarmDefineFieldVector has been called */
246: static PetscErrorCode DMCreateGlobalVector_Swarm(DM dm, Vec *vec)
247: {
248: DM_Swarm *swarm = (DM_Swarm *)dm->data;
249: Vec x;
250: char name[PETSC_MAX_PATH_LEN];
252: PetscFunctionBegin;
253: if (!swarm->issetup) PetscCall(DMSetUp(dm));
254: PetscCheck(swarm->vec_field_num, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmVectorDefineField(s) first");
255: PetscCheck(swarm->vec_field_nlocal == swarm->db->L, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSwarm sizes have changed since last call to VectorDefineField(s) first"); /* Stale data */
257: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
258: for (PetscInt f = 0; f < swarm->vec_field_num; ++f) {
259: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
260: PetscCall(PetscStrlcat(name, swarm->vec_field_names[f], PETSC_MAX_PATH_LEN));
261: }
262: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &x));
263: PetscCall(PetscObjectSetName((PetscObject)x, name));
264: PetscCall(VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE));
265: PetscCall(VecSetBlockSize(x, swarm->vec_field_bs));
266: PetscCall(VecSetDM(x, dm));
267: PetscCall(VecSetFromOptions(x));
268: PetscCall(VecSetDM(x, dm));
269: PetscCall(VecSetOperation(x, VECOP_VIEW, (void (*)(void))VecView_Swarm));
270: *vec = x;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /* requires DMSwarmDefineFieldVector has been called */
275: static PetscErrorCode DMCreateLocalVector_Swarm(DM dm, Vec *vec)
276: {
277: DM_Swarm *swarm = (DM_Swarm *)dm->data;
278: Vec x;
279: char name[PETSC_MAX_PATH_LEN];
281: PetscFunctionBegin;
282: if (!swarm->issetup) PetscCall(DMSetUp(dm));
283: PetscCheck(swarm->vec_field_num, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmVectorDefineField(s) first");
284: PetscCheck(swarm->vec_field_nlocal == swarm->db->L, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSwarm sizes have changed since last call to VectorDefineField(s) first");
286: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
287: for (PetscInt f = 0; f < swarm->vec_field_num; ++f) {
288: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
289: PetscCall(PetscStrlcat(name, swarm->vec_field_names[f], PETSC_MAX_PATH_LEN));
290: }
291: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
292: PetscCall(PetscObjectSetName((PetscObject)x, name));
293: PetscCall(VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE));
294: PetscCall(VecSetBlockSize(x, swarm->vec_field_bs));
295: PetscCall(VecSetDM(x, dm));
296: PetscCall(VecSetFromOptions(x));
297: *vec = x;
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
302: {
303: DM_Swarm *swarm = (DM_Swarm *)dm->data;
304: DMSwarmDataField gfield;
305: PetscInt bs, nlocal, fid = -1, cfid = -2;
306: PetscBool flg;
308: PetscFunctionBegin;
309: /* check vector is an inplace array */
310: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
311: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
312: (void)flg; /* avoid compiler warning */
313: 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);
314: PetscCall(VecGetLocalSize(*vec, &nlocal));
315: PetscCall(VecGetBlockSize(*vec, &bs));
316: 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");
317: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
318: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
319: PetscCall(VecResetArray(*vec));
320: PetscCall(VecDestroy(vec));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
325: {
326: DM_Swarm *swarm = (DM_Swarm *)dm->data;
327: PetscDataType type;
328: PetscScalar *array;
329: PetscInt bs, n, fid;
330: char name[PETSC_MAX_PATH_LEN];
331: PetscMPIInt size;
332: PetscBool iscuda, iskokkos;
334: PetscFunctionBegin;
335: if (!swarm->issetup) PetscCall(DMSetUp(dm));
336: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
337: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
338: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
340: PetscCallMPI(MPI_Comm_size(comm, &size));
341: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
342: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
343: PetscCall(VecCreate(comm, vec));
344: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
345: PetscCall(VecSetBlockSize(*vec, bs));
346: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
347: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
348: else PetscCall(VecSetType(*vec, VECSTANDARD));
349: PetscCall(VecPlaceArray(*vec, array));
351: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
352: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
354: /* Set guard */
355: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
356: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
358: PetscCall(VecSetDM(*vec, dm));
359: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
365: \hat f = \sum_i f_i \phi_i
367: and a particle function is given by
369: f = \sum_i w_i \delta(x - x_i)
371: then we want to require that
373: M \hat f = M_p f
375: where the particle mass matrix is given by
377: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
379: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
380: his integral. We allow this with the boolean flag.
381: */
382: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
383: {
384: const char *name = "Mass Matrix";
385: MPI_Comm comm;
386: PetscDS prob;
387: PetscSection fsection, globalFSection;
388: PetscHSetIJ ht;
389: PetscLayout rLayout, colLayout;
390: PetscInt *dnz, *onz;
391: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
392: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
393: PetscScalar *elemMat;
394: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0, totNc = 0;
395: const char *coordname;
397: PetscFunctionBegin;
398: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
399: PetscCall(DMGetCoordinateDim(dmf, &dim));
400: PetscCall(DMGetDS(dmf, &prob));
401: PetscCall(PetscDSGetNumFields(prob, &Nf));
402: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
403: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
404: PetscCall(DMGetLocalSection(dmf, &fsection));
405: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
406: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
407: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
409: PetscCall(DMSwarmGetCoordinateField(dmc, &coordname));
411: PetscCall(PetscLayoutCreate(comm, &colLayout));
412: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
413: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
414: PetscCall(PetscLayoutSetUp(colLayout));
415: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
416: PetscCall(PetscLayoutDestroy(&colLayout));
418: PetscCall(PetscLayoutCreate(comm, &rLayout));
419: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
420: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
421: PetscCall(PetscLayoutSetUp(rLayout));
422: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
423: PetscCall(PetscLayoutDestroy(&rLayout));
425: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
426: PetscCall(PetscHSetIJCreate(&ht));
428: PetscCall(PetscSynchronizedFlush(comm, NULL));
429: for (field = 0; field < Nf; ++field) {
430: PetscObject obj;
431: PetscClassId id;
432: PetscInt Nc;
434: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
435: PetscCall(PetscObjectGetClassId(obj, &id));
436: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
437: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
438: totNc += Nc;
439: }
440: /* count non-zeros */
441: PetscCall(DMSwarmSortGetAccess(dmc));
442: for (field = 0; field < Nf; ++field) {
443: PetscObject obj;
444: PetscClassId id;
445: PetscInt Nc;
447: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
448: PetscCall(PetscObjectGetClassId(obj, &id));
449: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
450: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
452: for (cell = cStart; cell < cEnd; ++cell) {
453: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
454: PetscInt numFIndices, numCIndices;
456: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
457: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
458: maxC = PetscMax(maxC, numCIndices);
459: {
460: PetscHashIJKey key;
461: PetscBool missing;
462: for (PetscInt i = 0; i < numFIndices; ++i) {
463: key.j = findices[i]; /* global column (from Plex) */
464: if (key.j >= 0) {
465: /* Get indices for coarse elements */
466: for (PetscInt j = 0; j < numCIndices; ++j) {
467: for (PetscInt c = 0; c < Nc; ++c) {
468: // TODO Need field offset on particle here
469: key.i = cindices[j] * totNc + c + rStart; /* global cols (from Swarm) */
470: if (key.i < 0) continue;
471: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
472: if (missing) {
473: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
474: else ++onz[key.i - rStart];
475: } else SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
476: }
477: }
478: }
479: }
480: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
481: }
482: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
483: }
484: }
485: PetscCall(PetscHSetIJDestroy(&ht));
486: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
487: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
488: PetscCall(PetscFree2(dnz, onz));
489: PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
490: for (field = 0; field < Nf; ++field) {
491: PetscTabulation Tcoarse;
492: PetscObject obj;
493: PetscClassId id;
494: PetscReal *fieldVals;
495: PetscInt Nc, bs;
497: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
498: PetscCall(PetscObjectGetClassId(obj, &id));
499: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
500: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
502: PetscCall(DMSwarmGetField(dmc, coordname, &bs, NULL, (void **)&fieldVals));
503: for (cell = cStart; cell < cEnd; ++cell) {
504: PetscInt *findices, *cindices;
505: PetscInt numFIndices, numCIndices;
507: /* TODO: Use DMField instead of assuming affine */
508: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
509: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
510: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
511: for (PetscInt j = 0; j < numCIndices; ++j) CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &fieldVals[cindices[j] * bs], &xi[j * dim]);
512: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
513: else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
514: /* Get elemMat entries by multiplying by weight */
515: PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
516: for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
517: for (PetscInt j = 0; j < numCIndices; ++j) {
518: for (PetscInt c = 0; c < Nc; ++c) {
519: // TODO Need field offset on particle and field here
520: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
521: elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
522: }
523: }
524: }
525: for (PetscInt j = 0; j < numCIndices; ++j)
526: // TODO Need field offset on particle here
527: for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
528: if (0) PetscCall(DMPrintCellMatrix(cell, name, numCIndices * Nc, numFIndices, elemMat));
529: PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
530: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
531: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
532: PetscCall(PetscTabulationDestroy(&Tcoarse));
533: }
534: PetscCall(DMSwarmRestoreField(dmc, coordname, &bs, NULL, (void **)&fieldVals));
535: }
536: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
537: PetscCall(DMSwarmSortRestoreAccess(dmc));
538: PetscCall(PetscFree3(v0, J, invJ));
539: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
540: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /* Returns empty matrix for use with SNES FD */
545: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
546: {
547: Vec field;
548: PetscInt size;
550: PetscFunctionBegin;
551: PetscCall(DMGetGlobalVector(sw, &field));
552: PetscCall(VecGetLocalSize(field, &size));
553: PetscCall(DMRestoreGlobalVector(sw, &field));
554: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
555: PetscCall(MatSetFromOptions(*m));
556: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
557: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
558: PetscCall(MatZeroEntries(*m));
559: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
560: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
561: PetscCall(MatShift(*m, 1.0));
562: PetscCall(MatSetDM(*m, sw));
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: /* FEM cols, Particle rows */
567: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
568: {
569: DM_Swarm *swarm = (DM_Swarm *)dmCoarse->data;
570: PetscSection gsf;
571: PetscInt m, n, Np;
572: void *ctx;
574: PetscFunctionBegin;
575: PetscCheck(swarm->vec_field_num, PetscObjectComm((PetscObject)dmCoarse), PETSC_ERR_USER, "Must call DMSwarmVectorDefineField(s) first");
576: PetscCall(DMGetGlobalSection(dmFine, &gsf));
577: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
578: PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
579: n = Np * swarm->vec_field_bs;
580: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
581: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
582: PetscCall(MatSetType(*mass, dmCoarse->mattype));
583: PetscCall(DMGetApplicationContext(dmFine, &ctx));
585: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
586: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
591: {
592: const char *name = "Mass Matrix Square";
593: MPI_Comm comm;
594: PetscDS prob;
595: PetscSection fsection, globalFSection;
596: PetscHSetIJ ht;
597: PetscLayout rLayout, colLayout;
598: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
599: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
600: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
601: PetscScalar *elemMat, *elemMatSq;
602: PetscInt cdim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
603: const char *coordname;
605: PetscFunctionBegin;
606: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
607: PetscCall(DMGetCoordinateDim(dmf, &cdim));
608: PetscCall(DMGetDS(dmf, &prob));
609: PetscCall(PetscDSGetNumFields(prob, &Nf));
610: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
611: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
612: PetscCall(DMGetLocalSection(dmf, &fsection));
613: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
614: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
615: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
617: PetscCall(DMSwarmGetCoordinateField(dmc, &coordname));
619: PetscCall(PetscLayoutCreate(comm, &colLayout));
620: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
621: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
622: PetscCall(PetscLayoutSetUp(colLayout));
623: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
624: PetscCall(PetscLayoutDestroy(&colLayout));
626: PetscCall(PetscLayoutCreate(comm, &rLayout));
627: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
628: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
629: PetscCall(PetscLayoutSetUp(rLayout));
630: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
631: PetscCall(PetscLayoutDestroy(&rLayout));
633: PetscCall(DMPlexGetDepth(dmf, &depth));
634: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
635: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
636: PetscCall(PetscMalloc1(maxAdjSize, &adj));
638: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
639: PetscCall(PetscHSetIJCreate(&ht));
640: /* Count nonzeros
641: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
642: */
643: PetscCall(DMSwarmSortGetAccess(dmc));
644: for (cell = cStart; cell < cEnd; ++cell) {
645: PetscInt i;
646: PetscInt *cindices;
647: PetscInt numCIndices;
648: #if 0
649: PetscInt adjSize = maxAdjSize, a, j;
650: #endif
652: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
653: maxC = PetscMax(maxC, numCIndices);
654: /* Diagonal block */
655: for (i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
656: #if 0
657: /* Off-diagonal blocks */
658: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
659: for (a = 0; a < adjSize; ++a) {
660: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
661: const PetscInt ncell = adj[a];
662: PetscInt *ncindices;
663: PetscInt numNCIndices;
665: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
666: {
667: PetscHashIJKey key;
668: PetscBool missing;
670: for (i = 0; i < numCIndices; ++i) {
671: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
672: if (key.i < 0) continue;
673: for (j = 0; j < numNCIndices; ++j) {
674: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
675: if (key.j < 0) continue;
676: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
677: if (missing) {
678: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
679: else ++onz[key.i - rStart];
680: }
681: }
682: }
683: }
684: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
685: }
686: }
687: #endif
688: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
689: }
690: PetscCall(PetscHSetIJDestroy(&ht));
691: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
692: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
693: PetscCall(PetscFree2(dnz, onz));
694: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
695: /* Fill in values
696: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
697: Start just by producing block diagonal
698: Could loop over adjacent cells
699: Produce neighboring element matrix
700: TODO Determine which columns and rows correspond to shared dual vector
701: Do MatMatMult with rectangular matrices
702: Insert block
703: */
704: for (field = 0; field < Nf; ++field) {
705: PetscTabulation Tcoarse;
706: PetscObject obj;
707: PetscReal *coords;
708: PetscInt Nc, i;
710: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
711: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
712: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
713: PetscCall(DMSwarmGetField(dmc, coordname, NULL, NULL, (void **)&coords));
714: for (cell = cStart; cell < cEnd; ++cell) {
715: PetscInt *findices, *cindices;
716: PetscInt numFIndices, numCIndices;
717: PetscInt p, c;
719: /* TODO: Use DMField instead of assuming affine */
720: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
721: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
722: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
723: for (p = 0; p < numCIndices; ++p) CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p] * cdim], &xi[p * cdim]);
724: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
725: /* Get elemMat entries by multiplying by weight */
726: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
727: for (i = 0; i < numFIndices; ++i) {
728: for (p = 0; p < numCIndices; ++p) {
729: for (c = 0; c < Nc; ++c) {
730: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
731: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
732: }
733: }
734: }
735: PetscCall(PetscTabulationDestroy(&Tcoarse));
736: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
737: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
738: /* Block diagonal */
739: if (numCIndices) {
740: PetscBLASInt blasn, blask;
741: PetscScalar one = 1.0, zero = 0.0;
743: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
744: PetscCall(PetscBLASIntCast(numFIndices, &blask));
745: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
746: }
747: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
748: /* TODO off-diagonal */
749: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
750: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
751: }
752: PetscCall(DMSwarmRestoreField(dmc, coordname, NULL, NULL, (void **)&coords));
753: }
754: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
755: PetscCall(PetscFree(adj));
756: PetscCall(DMSwarmSortRestoreAccess(dmc));
757: PetscCall(PetscFree3(v0, J, invJ));
758: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
759: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
760: PetscFunctionReturn(PETSC_SUCCESS);
761: }
763: /*@
764: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
766: Collective
768: Input Parameters:
769: + dmCoarse - a `DMSWARM`
770: - dmFine - a `DMPLEX`
772: Output Parameter:
773: . mass - the square of the particle mass matrix
775: Level: advanced
777: Note:
778: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
779: future to compute the full normal equations.
781: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
782: @*/
783: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
784: {
785: PetscInt n;
786: void *ctx;
788: PetscFunctionBegin;
789: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
790: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
791: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
792: PetscCall(MatSetType(*mass, dmCoarse->mattype));
793: PetscCall(DMGetApplicationContext(dmFine, &ctx));
795: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
796: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }
800: /*@
801: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
803: Collective
805: Input Parameters:
806: + dm - a `DMSWARM`
807: - fieldname - the textual name given to a registered field
809: Output Parameter:
810: . vec - the vector
812: Level: beginner
814: Note:
815: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
817: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
818: @*/
819: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
820: {
821: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
823: PetscFunctionBegin;
825: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }
829: /*@
830: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
832: Collective
834: Input Parameters:
835: + dm - a `DMSWARM`
836: - fieldname - the textual name given to a registered field
838: Output Parameter:
839: . vec - the vector
841: Level: beginner
843: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
844: @*/
845: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
846: {
847: PetscFunctionBegin;
849: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }
853: /*@
854: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
856: Collective
858: Input Parameters:
859: + dm - a `DMSWARM`
860: - fieldname - the textual name given to a registered field
862: Output Parameter:
863: . vec - the vector
865: Level: beginner
867: Note:
868: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
870: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
871: @*/
872: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
873: {
874: MPI_Comm comm = PETSC_COMM_SELF;
876: PetscFunctionBegin;
877: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: /*@
882: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
884: Collective
886: Input Parameters:
887: + dm - a `DMSWARM`
888: - fieldname - the textual name given to a registered field
890: Output Parameter:
891: . vec - the vector
893: Level: beginner
895: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
896: @*/
897: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
898: {
899: PetscFunctionBegin;
901: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
902: PetscFunctionReturn(PETSC_SUCCESS);
903: }
905: /*@
906: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
908: Collective
910: Input Parameter:
911: . dm - a `DMSWARM`
913: Level: beginner
915: Note:
916: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
918: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
919: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
920: @*/
921: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
922: {
923: DM_Swarm *swarm = (DM_Swarm *)dm->data;
925: PetscFunctionBegin;
926: if (!swarm->field_registration_initialized) {
927: swarm->field_registration_initialized = PETSC_TRUE;
928: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
929: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
930: }
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: /*@
935: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
937: Collective
939: Input Parameter:
940: . dm - a `DMSWARM`
942: Level: beginner
944: Note:
945: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
947: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
948: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
949: @*/
950: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
951: {
952: DM_Swarm *swarm = (DM_Swarm *)dm->data;
954: PetscFunctionBegin;
955: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
956: swarm->field_registration_finalized = PETSC_TRUE;
957: PetscFunctionReturn(PETSC_SUCCESS);
958: }
960: /*@
961: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
963: Not Collective
965: Input Parameters:
966: + dm - a `DMSWARM`
967: . nlocal - the length of each registered field
968: - buffer - the length of the buffer used to efficient dynamic re-sizing
970: Level: beginner
972: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
973: @*/
974: PetscErrorCode DMSwarmSetLocalSizes(DM dm, PetscInt nlocal, PetscInt buffer)
975: {
976: DM_Swarm *swarm = (DM_Swarm *)dm->data;
978: PetscFunctionBegin;
979: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
980: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
981: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
982: PetscFunctionReturn(PETSC_SUCCESS);
983: }
985: /*@
986: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
988: Collective
990: Input Parameters:
991: + dm - a `DMSWARM`
992: - dmcell - the `DM` to attach to the `DMSWARM`
994: Level: beginner
996: Note:
997: The attached `DM` (dmcell) will be queried for point location and
998: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1000: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
1001: @*/
1002: PetscErrorCode DMSwarmSetCellDM(DM dm, DM dmcell)
1003: {
1004: PetscFunctionBegin;
1007: PetscCall(DMSwarmPushCellDM(dm, dmcell, 0, NULL, DMSwarmPICField_coor));
1008: PetscFunctionReturn(PETSC_SUCCESS);
1009: }
1011: /*@
1012: DMSwarmGetCellDM - Fetches the attached cell `DM`
1014: Collective
1016: Input Parameter:
1017: . dm - a `DMSWARM`
1019: Output Parameter:
1020: . dmcell - the `DM` which was attached to the `DMSWARM`
1022: Level: beginner
1024: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1025: @*/
1026: PetscErrorCode DMSwarmGetCellDM(DM dm, DM *dmcell)
1027: {
1028: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1030: PetscFunctionBegin;
1031: *dmcell = swarm->cellinfo ? swarm->cellinfo[0].dm : NULL;
1032: PetscFunctionReturn(PETSC_SUCCESS);
1033: }
1035: static PetscErrorCode CellDMInfoDestroy(CellDMInfo *info)
1036: {
1037: PetscFunctionBegin;
1038: for (PetscInt f = 0; f < (*info)->Nf; ++f) PetscCall(PetscFree((*info)->dmFields[f]));
1039: PetscCall(PetscFree((*info)->dmFields));
1040: PetscCall(PetscFree((*info)->coordField));
1041: PetscCall(DMDestroy(&(*info)->dm));
1042: PetscCall(PetscFree(*info));
1043: PetscFunctionReturn(PETSC_SUCCESS);
1044: }
1046: /*@
1047: DMSwarmPushCellDM - Attaches a `DM` to a `DMSWARM`
1049: Collective
1051: Input Parameters:
1052: + sw - a `DMSWARM`
1053: . dmcell - the `DM` to attach to the `DMSWARM`
1054: . Nf - the number of swarm fields defining the `DM`
1055: . dmFields - an array of field names for the fields defining the `DM`
1056: - coordField - the name for the swarm field to use for particle coordinates on the cell `DM`
1058: Level: beginner
1060: Note:
1061: The attached `DM` (dmcell) will be queried for point location and
1062: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1064: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPopCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1065: @*/
1066: PetscErrorCode DMSwarmPushCellDM(DM sw, DM dmcell, PetscInt Nf, const char *dmFields[], const char coordField[])
1067: {
1068: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1069: CellDMInfo info;
1070: PetscBool rebin = swarm->cellinfo ? PETSC_TRUE : PETSC_FALSE;
1072: PetscFunctionBegin;
1075: if (Nf) PetscAssertPointer(dmFields, 4);
1076: PetscAssertPointer(coordField, 5);
1077: PetscCall(PetscNew(&info));
1078: PetscCall(PetscObjectReference((PetscObject)dmcell));
1079: info->dm = dmcell;
1080: info->Nf = Nf;
1081: info->next = swarm->cellinfo;
1082: swarm->cellinfo = info;
1083: // Define the DM fields
1084: PetscCall(PetscMalloc1(info->Nf, &info->dmFields));
1085: for (PetscInt f = 0; f < info->Nf; ++f) PetscCall(PetscStrallocpy(dmFields[f], &info->dmFields[f]));
1086: if (info->Nf) PetscCall(DMSwarmVectorDefineFields(sw, info->Nf, (const char **)info->dmFields));
1087: // Set the coordinate field
1088: PetscCall(PetscStrallocpy(coordField, &info->coordField));
1089: if (info->coordField) PetscCall(DMSwarmSetCoordinateField(sw, info->coordField));
1090: // Rebin the cells and set cell_id field
1091: if (rebin) PetscCall(DMSwarmMigrate(sw, PETSC_FALSE));
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: /*@
1096: DMSwarmPopCellDM - Discard the current cell `DM` and restore the previous one, if it exists
1098: Collective
1100: Input Parameter:
1101: . sw - a `DMSWARM`
1103: Level: beginner
1105: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPushCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1106: @*/
1107: PetscErrorCode DMSwarmPopCellDM(DM sw)
1108: {
1109: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1110: CellDMInfo info = swarm->cellinfo;
1112: PetscFunctionBegin;
1114: if (!swarm->cellinfo) PetscFunctionReturn(PETSC_SUCCESS);
1115: if (info->next) {
1116: CellDMInfo newinfo = info->next;
1118: swarm->cellinfo = info->next;
1119: // Define the DM fields
1120: PetscCall(DMSwarmVectorDefineFields(sw, newinfo->Nf, (const char **)newinfo->dmFields));
1121: // Set the coordinate field
1122: PetscCall(DMSwarmSetCoordinateField(sw, newinfo->coordField));
1123: // Rebin the cells and set cell_id field
1124: PetscCall(DMSwarmMigrate(sw, PETSC_FALSE));
1125: }
1126: PetscCall(CellDMInfoDestroy(&info));
1127: PetscFunctionReturn(PETSC_SUCCESS);
1128: }
1130: /*@
1131: DMSwarmGetLocalSize - Retrieves the local length of fields registered
1133: Not Collective
1135: Input Parameter:
1136: . dm - a `DMSWARM`
1138: Output Parameter:
1139: . nlocal - the length of each registered field
1141: Level: beginner
1143: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
1144: @*/
1145: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
1146: {
1147: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1149: PetscFunctionBegin;
1150: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
1151: PetscFunctionReturn(PETSC_SUCCESS);
1152: }
1154: /*@
1155: DMSwarmGetSize - Retrieves the total length of fields registered
1157: Collective
1159: Input Parameter:
1160: . dm - a `DMSWARM`
1162: Output Parameter:
1163: . n - the total length of each registered field
1165: Level: beginner
1167: Note:
1168: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
1170: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
1171: @*/
1172: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
1173: {
1174: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1175: PetscInt nlocal;
1177: PetscFunctionBegin;
1178: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1179: PetscCallMPI(MPIU_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1180: PetscFunctionReturn(PETSC_SUCCESS);
1181: }
1183: /*@
1184: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
1186: Collective
1188: Input Parameters:
1189: + dm - a `DMSWARM`
1190: . fieldname - the textual name to identify this field
1191: . blocksize - the number of each data type
1192: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
1194: Level: beginner
1196: Notes:
1197: The textual name for each registered field must be unique.
1199: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1200: @*/
1201: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1202: {
1203: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1204: size_t size;
1206: PetscFunctionBegin;
1207: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1208: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1210: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1211: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1212: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1213: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1214: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1216: PetscCall(PetscDataTypeGetSize(type, &size));
1217: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1218: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1219: {
1220: DMSwarmDataField gfield;
1222: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1223: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1224: }
1225: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1226: PetscFunctionReturn(PETSC_SUCCESS);
1227: }
1229: /*@C
1230: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
1232: Collective
1234: Input Parameters:
1235: + dm - a `DMSWARM`
1236: . fieldname - the textual name to identify this field
1237: - size - the size in bytes of the user struct of each data type
1239: Level: beginner
1241: Note:
1242: The textual name for each registered field must be unique.
1244: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1245: @*/
1246: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1247: {
1248: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1250: PetscFunctionBegin;
1251: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1252: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: /*@C
1257: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1259: Collective
1261: Input Parameters:
1262: + dm - a `DMSWARM`
1263: . fieldname - the textual name to identify this field
1264: . size - the size in bytes of the user data type
1265: - blocksize - the number of each data type
1267: Level: beginner
1269: Note:
1270: The textual name for each registered field must be unique.
1272: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1273: @*/
1274: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1275: {
1276: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1278: PetscFunctionBegin;
1279: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1280: {
1281: DMSwarmDataField gfield;
1283: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1284: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1285: }
1286: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: /*@C
1291: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1293: Not Collective, No Fortran Support
1295: Input Parameters:
1296: + dm - a `DMSWARM`
1297: - fieldname - the textual name to identify this field
1299: Output Parameters:
1300: + blocksize - the number of each data type
1301: . type - the data type
1302: - data - pointer to raw array
1304: Level: beginner
1306: Notes:
1307: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1309: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1310: @*/
1311: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1312: {
1313: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1314: DMSwarmDataField gfield;
1316: PetscFunctionBegin;
1318: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1319: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1320: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1321: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1322: if (blocksize) *blocksize = gfield->bs;
1323: if (type) *type = gfield->petsc_type;
1324: PetscFunctionReturn(PETSC_SUCCESS);
1325: }
1327: /*@C
1328: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1330: Not Collective, No Fortran Support
1332: Input Parameters:
1333: + dm - a `DMSWARM`
1334: - fieldname - the textual name to identify this field
1336: Output Parameters:
1337: + blocksize - the number of each data type
1338: . type - the data type
1339: - data - pointer to raw array
1341: Level: beginner
1343: Notes:
1344: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1346: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1347: @*/
1348: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1349: {
1350: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1351: DMSwarmDataField gfield;
1353: PetscFunctionBegin;
1355: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1356: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1357: if (data) *data = NULL;
1358: PetscFunctionReturn(PETSC_SUCCESS);
1359: }
1361: /*@
1362: DMSwarmGetCoordinateField - Get the name of the field holding particle coordinates
1364: Not Collective
1366: Input Parameter:
1367: . sw - a `DMSWARM`
1369: Output Parameter:
1370: . fieldname - the name of the coordinate field
1372: Level: intermediate
1374: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCoordinateField()`, `DMSwarmGetField()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1375: @*/
1376: PetscErrorCode DMSwarmGetCoordinateField(DM sw, const char *fieldname[])
1377: {
1378: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1380: PetscFunctionBegin;
1382: PetscAssertPointer(fieldname, 2);
1383: *fieldname = swarm->coord_name;
1384: PetscFunctionReturn(PETSC_SUCCESS);
1385: }
1387: /*@
1388: DMSwarmSetCoordinateField - Set the name of the field holding particle coordinates
1390: Not Collective
1392: Input Parameters:
1393: + sw - a `DMSWARM`
1394: - fieldname - the name of the coordinate field
1396: Level: intermediate
1398: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCoordinateField()`, `DMSwarmGetField()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1399: @*/
1400: PetscErrorCode DMSwarmSetCoordinateField(DM sw, const char fieldname[])
1401: {
1402: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1404: PetscFunctionBegin;
1406: PetscAssertPointer(fieldname, 2);
1407: PetscCall(PetscFree(swarm->coord_name));
1408: PetscCall(PetscStrallocpy(fieldname, (char **)&swarm->coord_name));
1409: PetscFunctionReturn(PETSC_SUCCESS);
1410: }
1412: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
1413: {
1414: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1415: DMSwarmDataField gfield;
1417: PetscFunctionBegin;
1419: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1420: if (blocksize) *blocksize = gfield->bs;
1421: if (type) *type = gfield->petsc_type;
1422: PetscFunctionReturn(PETSC_SUCCESS);
1423: }
1425: /*@
1426: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
1428: Not Collective
1430: Input Parameter:
1431: . dm - a `DMSWARM`
1433: Level: beginner
1435: Notes:
1436: The new point will have all fields initialized to zero.
1438: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
1439: @*/
1440: PetscErrorCode DMSwarmAddPoint(DM dm)
1441: {
1442: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1444: PetscFunctionBegin;
1445: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1446: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1447: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
1448: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1449: PetscFunctionReturn(PETSC_SUCCESS);
1450: }
1452: /*@
1453: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
1455: Not Collective
1457: Input Parameters:
1458: + dm - a `DMSWARM`
1459: - npoints - the number of new points to add
1461: Level: beginner
1463: Notes:
1464: The new point will have all fields initialized to zero.
1466: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
1467: @*/
1468: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1469: {
1470: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1471: PetscInt nlocal;
1473: PetscFunctionBegin;
1474: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1475: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1476: nlocal = nlocal + npoints;
1477: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
1478: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1479: PetscFunctionReturn(PETSC_SUCCESS);
1480: }
1482: /*@
1483: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
1485: Not Collective
1487: Input Parameter:
1488: . dm - a `DMSWARM`
1490: Level: beginner
1492: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
1493: @*/
1494: PetscErrorCode DMSwarmRemovePoint(DM dm)
1495: {
1496: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1498: PetscFunctionBegin;
1499: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1500: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
1501: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1502: PetscFunctionReturn(PETSC_SUCCESS);
1503: }
1505: /*@
1506: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
1508: Not Collective
1510: Input Parameters:
1511: + dm - a `DMSWARM`
1512: - idx - index of point to remove
1514: Level: beginner
1516: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1517: @*/
1518: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
1519: {
1520: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1522: PetscFunctionBegin;
1523: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1524: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
1525: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: /*@
1530: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
1532: Not Collective
1534: Input Parameters:
1535: + dm - a `DMSWARM`
1536: . pi - the index of the point to copy
1537: - pj - the point index where the copy should be located
1539: Level: beginner
1541: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1542: @*/
1543: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
1544: {
1545: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1547: PetscFunctionBegin;
1548: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1549: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
1550: PetscFunctionReturn(PETSC_SUCCESS);
1551: }
1553: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
1554: {
1555: PetscFunctionBegin;
1556: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: /*@
1561: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
1563: Collective
1565: Input Parameters:
1566: + dm - the `DMSWARM`
1567: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1569: Level: advanced
1571: Notes:
1572: The `DM` will be modified to accommodate received points.
1573: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
1574: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
1576: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
1577: @*/
1578: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
1579: {
1580: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1582: PetscFunctionBegin;
1583: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
1584: switch (swarm->migrate_type) {
1585: case DMSWARM_MIGRATE_BASIC:
1586: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
1587: break;
1588: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1589: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
1590: break;
1591: case DMSWARM_MIGRATE_DMCELLEXACT:
1592: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1593: case DMSWARM_MIGRATE_USER:
1594: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
1595: default:
1596: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
1597: }
1598: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
1599: PetscCall(DMClearGlobalVectors(dm));
1600: PetscFunctionReturn(PETSC_SUCCESS);
1601: }
1603: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
1605: /*
1606: DMSwarmCollectViewCreate
1608: * Applies a collection method and gathers point neighbour points into dm
1610: Notes:
1611: Users should call DMSwarmCollectViewDestroy() after
1612: they have finished computations associated with the collected points
1613: */
1615: /*@
1616: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1617: in neighbour ranks into the `DMSWARM`
1619: Collective
1621: Input Parameter:
1622: . dm - the `DMSWARM`
1624: Level: advanced
1626: Notes:
1627: Users should call `DMSwarmCollectViewDestroy()` after
1628: they have finished computations associated with the collected points
1630: Different collect methods are supported. See `DMSwarmSetCollectType()`.
1632: Developer Note:
1633: Create and Destroy routines create new objects that can get destroyed, they do not change the state
1634: of the current object.
1636: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
1637: @*/
1638: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1639: {
1640: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1641: PetscInt ng;
1643: PetscFunctionBegin;
1644: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
1645: PetscCall(DMSwarmGetLocalSize(dm, &ng));
1646: switch (swarm->collect_type) {
1647: case DMSWARM_COLLECT_BASIC:
1648: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
1649: break;
1650: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1651: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1652: case DMSWARM_COLLECT_GENERAL:
1653: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
1654: default:
1655: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
1656: }
1657: swarm->collect_view_active = PETSC_TRUE;
1658: swarm->collect_view_reset_nlocal = ng;
1659: PetscFunctionReturn(PETSC_SUCCESS);
1660: }
1662: /*@
1663: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
1665: Collective
1667: Input Parameters:
1668: . dm - the `DMSWARM`
1670: Notes:
1671: Users should call `DMSwarmCollectViewCreate()` before this function is called.
1673: Level: advanced
1675: Developer Note:
1676: Create and Destroy routines create new objects that can get destroyed, they do not change the state
1677: of the current object.
1679: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
1680: @*/
1681: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1682: {
1683: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1685: PetscFunctionBegin;
1686: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
1687: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
1688: swarm->collect_view_active = PETSC_FALSE;
1689: PetscFunctionReturn(PETSC_SUCCESS);
1690: }
1692: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
1693: {
1694: PetscInt dim;
1696: PetscFunctionBegin;
1697: PetscCall(DMSwarmSetNumSpecies(dm, 1));
1698: PetscCall(DMGetDimension(dm, &dim));
1699: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1700: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1701: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_coor, dim, PETSC_DOUBLE));
1702: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_cellid, 1, PETSC_INT));
1703: PetscCall(DMSwarmSetCoordinateField(dm, DMSwarmPICField_coor));
1704: PetscFunctionReturn(PETSC_SUCCESS);
1705: }
1707: /*@
1708: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
1710: Collective
1712: Input Parameters:
1713: + dm - the `DMSWARM`
1714: - Npc - The number of particles per cell in the cell `DM`
1716: Level: intermediate
1718: Notes:
1719: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
1720: one particle is in each cell, it is placed at the centroid.
1722: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1723: @*/
1724: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1725: {
1726: DM cdm;
1727: PetscRandom rnd;
1728: DMPolytopeType ct;
1729: PetscBool simplex;
1730: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1731: PetscInt dim, d, cStart, cEnd, c, p;
1732: const char *coordname;
1734: PetscFunctionBeginUser;
1735: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
1736: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
1737: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
1739: PetscCall(DMSwarmGetCoordinateField(dm, &coordname));
1740: PetscCall(DMSwarmGetCellDM(dm, &cdm));
1741: PetscCall(DMGetDimension(cdm, &dim));
1742: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
1743: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
1744: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1746: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
1747: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1748: PetscCall(DMSwarmGetField(dm, coordname, NULL, NULL, (void **)&coords));
1749: for (c = cStart; c < cEnd; ++c) {
1750: if (Npc == 1) {
1751: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
1752: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
1753: } else {
1754: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
1755: for (p = 0; p < Npc; ++p) {
1756: const PetscInt n = c * Npc + p;
1757: PetscReal sum = 0.0, refcoords[3];
1759: for (d = 0; d < dim; ++d) {
1760: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
1761: sum += refcoords[d];
1762: }
1763: if (simplex && sum > 0.0)
1764: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
1765: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
1766: }
1767: }
1768: }
1769: PetscCall(DMSwarmRestoreField(dm, coordname, NULL, NULL, (void **)&coords));
1770: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
1771: PetscCall(PetscRandomDestroy(&rnd));
1772: PetscFunctionReturn(PETSC_SUCCESS);
1773: }
1775: /*@
1776: DMSwarmSetType - Set particular flavor of `DMSWARM`
1778: Collective
1780: Input Parameters:
1781: + dm - the `DMSWARM`
1782: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
1784: Level: advanced
1786: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
1787: @*/
1788: PetscErrorCode DMSwarmSetType(DM dm, DMSwarmType stype)
1789: {
1790: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1792: PetscFunctionBegin;
1793: swarm->swarm_type = stype;
1794: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(dm));
1795: PetscFunctionReturn(PETSC_SUCCESS);
1796: }
1798: static PetscErrorCode DMSetup_Swarm(DM dm)
1799: {
1800: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1801: PetscMPIInt rank;
1802: PetscInt p, npoints, *rankval;
1804: PetscFunctionBegin;
1805: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
1806: swarm->issetup = PETSC_TRUE;
1808: if (swarm->swarm_type == DMSWARM_PIC) {
1809: /* check dmcell exists */
1810: PetscCheck(swarm->cellinfo && swarm->cellinfo[0].dm, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires you call DMSwarmSetCellDM() or DMSwarmPushCellDM()");
1812: if (swarm->cellinfo[0].dm->ops->locatepointssubdomain) {
1813: /* check methods exists for exact ownership identificiation */
1814: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
1815: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1816: } else {
1817: /* check methods exist for point location AND rank neighbor identification */
1818: if (swarm->cellinfo[0].dm->ops->locatepoints) {
1819: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
1820: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1822: if (swarm->cellinfo[0].dm->ops->getneighbors) {
1823: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
1824: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1826: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1827: }
1828: }
1830: PetscCall(DMSwarmFinalizeFieldRegister(dm));
1832: /* check some fields were registered */
1833: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
1835: /* check local sizes were set */
1836: PetscCheck(swarm->db->L != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Local sizes must be set via DMSwarmSetLocalSizes()");
1838: /* initialize values in pid and rank placeholders */
1839: /* TODO: [pid - use MPI_Scan] */
1840: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1841: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
1842: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1843: for (p = 0; p < npoints; p++) rankval[p] = rank;
1844: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1845: PetscFunctionReturn(PETSC_SUCCESS);
1846: }
1848: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1850: static PetscErrorCode DMDestroy_Swarm(DM dm)
1851: {
1852: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1853: CellDMInfo info = swarm->cellinfo;
1855: PetscFunctionBegin;
1856: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
1857: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
1858: for (PetscInt f = 0; f < swarm->vec_field_num; ++f) PetscCall(PetscFree(swarm->vec_field_names[f]));
1859: PetscCall(PetscFree(swarm->vec_field_names));
1860: PetscCall(PetscFree(swarm->coord_name));
1861: if (swarm->sort_context) PetscCall(DMSwarmSortDestroy(&swarm->sort_context));
1862: while (info) {
1863: CellDMInfo tmp = info;
1865: info = info->next;
1866: PetscCall(CellDMInfoDestroy(&tmp));
1867: }
1868: PetscCall(PetscFree(swarm));
1869: PetscFunctionReturn(PETSC_SUCCESS);
1870: }
1872: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1873: {
1874: DM cdm;
1875: PetscDraw draw;
1876: PetscReal *coords, oldPause, radius = 0.01;
1877: PetscInt Np, p, bs;
1878: const char *coordname;
1880: PetscFunctionBegin;
1881: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
1882: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1883: PetscCall(DMSwarmGetCellDM(dm, &cdm));
1884: PetscCall(PetscDrawGetPause(draw, &oldPause));
1885: PetscCall(PetscDrawSetPause(draw, 0.0));
1886: PetscCall(DMView(cdm, viewer));
1887: PetscCall(PetscDrawSetPause(draw, oldPause));
1889: PetscCall(DMSwarmGetCoordinateField(dm, &coordname));
1890: PetscCall(DMSwarmGetLocalSize(dm, &Np));
1891: PetscCall(DMSwarmGetField(dm, coordname, &bs, NULL, (void **)&coords));
1892: for (p = 0; p < Np; ++p) {
1893: const PetscInt i = p * bs;
1895: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
1896: }
1897: PetscCall(DMSwarmRestoreField(dm, coordname, &bs, NULL, (void **)&coords));
1898: PetscCall(PetscDrawFlush(draw));
1899: PetscCall(PetscDrawPause(draw));
1900: PetscCall(PetscDrawSave(draw));
1901: PetscFunctionReturn(PETSC_SUCCESS);
1902: }
1904: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
1905: {
1906: PetscViewerFormat format;
1907: PetscInt *sizes;
1908: PetscInt dim, Np, maxSize = 17;
1909: MPI_Comm comm;
1910: PetscMPIInt rank, size, p;
1911: const char *name;
1913: PetscFunctionBegin;
1914: PetscCall(PetscViewerGetFormat(viewer, &format));
1915: PetscCall(DMGetDimension(dm, &dim));
1916: PetscCall(DMSwarmGetLocalSize(dm, &Np));
1917: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1918: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1919: PetscCallMPI(MPI_Comm_size(comm, &size));
1920: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1921: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
1922: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
1923: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
1924: else PetscCall(PetscCalloc1(3, &sizes));
1925: if (size < maxSize) {
1926: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
1927: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
1928: for (p = 0; p < size; ++p) {
1929: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
1930: }
1931: } else {
1932: PetscInt locMinMax[2] = {Np, Np};
1934: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
1935: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
1936: }
1937: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1938: PetscCall(PetscFree(sizes));
1939: if (format == PETSC_VIEWER_ASCII_INFO) {
1940: PetscInt *cell;
1942: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
1943: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1944: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell));
1945: for (p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%d: %" PetscInt_FMT "\n", p, cell[p]));
1946: PetscCall(PetscViewerFlush(viewer));
1947: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1948: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell));
1949: }
1950: PetscFunctionReturn(PETSC_SUCCESS);
1951: }
1953: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1954: {
1955: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1956: PetscBool iascii, ibinary, isvtk, isdraw, ispython;
1957: #if defined(PETSC_HAVE_HDF5)
1958: PetscBool ishdf5;
1959: #endif
1961: PetscFunctionBegin;
1964: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1965: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
1966: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1967: #if defined(PETSC_HAVE_HDF5)
1968: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1969: #endif
1970: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1971: PetscCall(PetscObjectHasFunction((PetscObject)viewer, "PetscViewerPythonViewObject_C", &ispython));
1972: if (iascii) {
1973: PetscViewerFormat format;
1975: PetscCall(PetscViewerGetFormat(viewer, &format));
1976: switch (format) {
1977: case PETSC_VIEWER_ASCII_INFO_DETAIL:
1978: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
1979: break;
1980: default:
1981: PetscCall(DMView_Swarm_Ascii(dm, viewer));
1982: }
1983: } else {
1984: #if defined(PETSC_HAVE_HDF5)
1985: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
1986: #endif
1987: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
1988: if (ispython) PetscCall(PetscViewerPythonViewObject(viewer, (PetscObject)dm));
1989: }
1990: PetscFunctionReturn(PETSC_SUCCESS);
1991: }
1993: /*@
1994: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
1995: 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.
1997: Noncollective
1999: Input Parameters:
2000: + sw - the `DMSWARM`
2001: . cellID - the integer id of the cell to be extracted and filtered
2002: - cellswarm - The `DMSWARM` to receive the cell
2004: Level: beginner
2006: Notes:
2007: This presently only supports `DMSWARM_PIC` type
2009: Should be restored with `DMSwarmRestoreCellSwarm()`
2011: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
2013: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
2014: @*/
2015: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2016: {
2017: DM_Swarm *original = (DM_Swarm *)sw->data;
2018: DMLabel label;
2019: DM dmc, subdmc;
2020: PetscInt *pids, particles, dim;
2022: PetscFunctionBegin;
2023: /* Configure new swarm */
2024: PetscCall(DMSetType(cellswarm, DMSWARM));
2025: PetscCall(DMGetDimension(sw, &dim));
2026: PetscCall(DMSetDimension(cellswarm, dim));
2027: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
2028: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
2029: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
2030: PetscCall(DMSwarmSortGetAccess(sw));
2031: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
2032: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2033: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
2034: PetscCall(DMSwarmSortRestoreAccess(sw));
2035: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2036: PetscCall(DMSwarmGetCellDM(sw, &dmc));
2037: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
2038: PetscCall(DMAddLabel(dmc, label));
2039: PetscCall(DMLabelSetValue(label, cellID, 1));
2040: PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdmc));
2041: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
2042: PetscCall(DMLabelDestroy(&label));
2043: PetscFunctionReturn(PETSC_SUCCESS);
2044: }
2046: /*@
2047: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
2049: Noncollective
2051: Input Parameters:
2052: + sw - the parent `DMSWARM`
2053: . cellID - the integer id of the cell to be copied back into the parent swarm
2054: - cellswarm - the cell swarm object
2056: Level: beginner
2058: Note:
2059: This only supports `DMSWARM_PIC` types of `DMSWARM`s
2061: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
2062: @*/
2063: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2064: {
2065: DM dmc;
2066: PetscInt *pids, particles, p;
2068: PetscFunctionBegin;
2069: PetscCall(DMSwarmSortGetAccess(sw));
2070: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2071: PetscCall(DMSwarmSortRestoreAccess(sw));
2072: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
2073: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
2074: /* Free memory, destroy cell dm */
2075: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
2076: PetscCall(DMDestroy(&dmc));
2077: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2078: PetscFunctionReturn(PETSC_SUCCESS);
2079: }
2081: /*@
2082: DMSwarmComputeMoments - Compute the first three particle moments for a given field
2084: Noncollective
2086: Input Parameters:
2087: + sw - the `DMSWARM`
2088: . coordinate - the coordinate field name
2089: - weight - the weight field name
2091: Output Parameter:
2092: . moments - the field moments
2094: Level: intermediate
2096: Notes:
2097: The `moments` array should be of length bs + 2, where bs is the block size of the coordinate field.
2099: The weight field must be a scalar, having blocksize 1.
2101: .seealso: `DM`, `DMSWARM`, `DMPlexComputeMoments()`
2102: @*/
2103: PetscErrorCode DMSwarmComputeMoments(DM sw, const char coordinate[], const char weight[], PetscReal moments[])
2104: {
2105: const PetscReal *coords;
2106: const PetscReal *w;
2107: PetscReal *mom;
2108: PetscDataType dtc, dtw;
2109: PetscInt bsc, bsw, Np;
2110: MPI_Comm comm;
2112: PetscFunctionBegin;
2114: PetscAssertPointer(coordinate, 2);
2115: PetscAssertPointer(weight, 3);
2116: PetscAssertPointer(moments, 4);
2117: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
2118: PetscCall(DMSwarmGetField(sw, coordinate, &bsc, &dtc, (void **)&coords));
2119: PetscCall(DMSwarmGetField(sw, weight, &bsw, &dtw, (void **)&w));
2120: PetscCheck(dtc == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s must be real, not %s", coordinate, PetscDataTypes[dtc]);
2121: PetscCheck(dtw == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be real, not %s", weight, PetscDataTypes[dtw]);
2122: PetscCheck(bsw == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be a scalar, not blocksize %" PetscInt_FMT, weight, bsw);
2123: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2124: PetscCall(DMGetWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2125: PetscCall(PetscArrayzero(mom, bsc + 2));
2126: for (PetscInt p = 0; p < Np; ++p) {
2127: const PetscReal *c = &coords[p * bsc];
2128: const PetscReal wp = w[p];
2130: mom[0] += wp;
2131: for (PetscInt d = 0; d < bsc; ++d) {
2132: mom[d + 1] += wp * c[d];
2133: mom[d + bsc + 1] += wp * PetscSqr(c[d]);
2134: }
2135: }
2136: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&coords));
2137: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
2138: PetscCallMPI(MPIU_Allreduce(mom, moments, bsc + 2, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
2139: PetscCall(DMRestoreWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2140: PetscFunctionReturn(PETSC_SUCCESS);
2141: }
2143: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
2145: static PetscErrorCode DMInitialize_Swarm(DM sw)
2146: {
2147: PetscFunctionBegin;
2148: sw->ops->view = DMView_Swarm;
2149: sw->ops->load = NULL;
2150: sw->ops->setfromoptions = NULL;
2151: sw->ops->clone = DMClone_Swarm;
2152: sw->ops->setup = DMSetup_Swarm;
2153: sw->ops->createlocalsection = NULL;
2154: sw->ops->createsectionpermutation = NULL;
2155: sw->ops->createdefaultconstraints = NULL;
2156: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
2157: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
2158: sw->ops->getlocaltoglobalmapping = NULL;
2159: sw->ops->createfieldis = NULL;
2160: sw->ops->createcoordinatedm = NULL;
2161: sw->ops->getcoloring = NULL;
2162: sw->ops->creatematrix = DMCreateMatrix_Swarm;
2163: sw->ops->createinterpolation = NULL;
2164: sw->ops->createinjection = NULL;
2165: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
2166: sw->ops->refine = NULL;
2167: sw->ops->coarsen = NULL;
2168: sw->ops->refinehierarchy = NULL;
2169: sw->ops->coarsenhierarchy = NULL;
2170: sw->ops->globaltolocalbegin = NULL;
2171: sw->ops->globaltolocalend = NULL;
2172: sw->ops->localtoglobalbegin = NULL;
2173: sw->ops->localtoglobalend = NULL;
2174: sw->ops->destroy = DMDestroy_Swarm;
2175: sw->ops->createsubdm = NULL;
2176: sw->ops->getdimpoints = NULL;
2177: sw->ops->locatepoints = NULL;
2178: PetscFunctionReturn(PETSC_SUCCESS);
2179: }
2181: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
2182: {
2183: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2185: PetscFunctionBegin;
2186: swarm->refct++;
2187: (*newdm)->data = swarm;
2188: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
2189: PetscCall(DMInitialize_Swarm(*newdm));
2190: (*newdm)->dim = dm->dim;
2191: PetscFunctionReturn(PETSC_SUCCESS);
2192: }
2194: /*MC
2195: DMSWARM = "swarm" - A `DM` object used to represent arrays of data (fields) of arbitrary data type.
2196: This implementation was designed for particle methods in which the underlying
2197: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
2199: Level: intermediate
2201: Notes:
2202: User data can be represented by `DMSWARM` through a registering "fields".
2203: To register a field, the user must provide;
2204: (a) a unique name;
2205: (b) the data type (or size in bytes);
2206: (c) the block size of the data.
2208: For example, suppose the application requires a unique id, energy, momentum and density to be stored
2209: on a set of particles. Then the following code could be used
2210: .vb
2211: DMSwarmInitializeFieldRegister(dm)
2212: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
2213: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
2214: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
2215: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
2216: DMSwarmFinalizeFieldRegister(dm)
2217: .ve
2219: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
2220: The only restriction imposed by `DMSWARM` is that all fields contain the same number of points.
2222: To support particle methods, "migration" techniques are provided. These methods migrate data
2223: between ranks.
2225: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
2226: As a `DMSWARM` may internally define and store values of different data types,
2227: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
2228: fields should be used to define a `Vec` object via
2229: `DMSwarmVectorDefineField()`
2230: The specified field can be changed at any time - thereby permitting vectors
2231: compatible with different fields to be created.
2233: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via
2234: `DMSwarmCreateGlobalVectorFromField()`
2235: Here the data defining the field in the `DMSWARM` is shared with a Vec.
2236: This is inherently unsafe if you alter the size of the field at any time between
2237: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
2238: If the local size of the `DMSWARM` does not match the local size of the global vector
2239: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
2241: Additional high-level support is provided for Particle-In-Cell methods.
2242: Please refer to `DMSwarmSetType()`.
2244: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`
2245: M*/
2247: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
2248: {
2249: DM_Swarm *swarm;
2251: PetscFunctionBegin;
2253: PetscCall(PetscNew(&swarm));
2254: dm->data = swarm;
2255: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
2256: PetscCall(DMSwarmInitializeFieldRegister(dm));
2257: dm->dim = 0;
2258: swarm->refct = 1;
2259: swarm->issetup = PETSC_FALSE;
2260: swarm->swarm_type = DMSWARM_BASIC;
2261: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
2262: swarm->collect_type = DMSWARM_COLLECT_BASIC;
2263: swarm->migrate_error_on_missing_point = PETSC_FALSE;
2264: swarm->cellinfo = NULL;
2265: swarm->collect_view_active = PETSC_FALSE;
2266: swarm->collect_view_reset_nlocal = -1;
2267: PetscCall(DMInitialize_Swarm(dm));
2268: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
2269: PetscFunctionReturn(PETSC_SUCCESS);
2270: }