Actual source code: swarm.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petsc/private/hashsetij.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petscviewer.h>
6: #include <petscdraw.h>
7: #include <petscdmplex.h>
8: #include <petscblaslapack.h>
9: #include "../src/dm/impls/swarm/data_bucket.h"
10: #include <petscdmlabel.h>
11: #include <petscsection.h>
13: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
14: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
15: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
17: const char *DMSwarmTypeNames[] = {"basic", "pic", NULL};
18: const char *DMSwarmMigrateTypeNames[] = {"basic", "dmcellnscatter", "dmcellexact", "user", NULL};
19: const char *DMSwarmCollectTypeNames[] = {"basic", "boundingbox", "general", "user", NULL};
20: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
22: const char DMSwarmField_pid[] = "DMSwarm_pid";
23: const char DMSwarmField_rank[] = "DMSwarm_rank";
24: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
25: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
27: PetscInt SwarmDataFieldId = -1;
29: #if defined(PETSC_HAVE_HDF5)
30: #include <petscviewerhdf5.h>
32: static PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
33: {
34: DM dm;
35: PetscReal seqval;
36: PetscInt seqnum, bs;
37: PetscBool isseq, ists;
39: PetscFunctionBegin;
40: PetscCall(VecGetDM(v, &dm));
41: PetscCall(VecGetBlockSize(v, &bs));
42: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particle_fields"));
43: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
44: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
45: PetscCall(PetscViewerHDF5IsTimestepping(viewer, &ists));
46: if (ists) PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
47: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
48: PetscCall(VecViewNative(v, viewer));
49: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
50: PetscCall(PetscViewerHDF5PopGroup(viewer));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
55: {
56: Vec coordinates;
57: PetscInt Np;
58: PetscBool isseq;
60: PetscFunctionBegin;
61: PetscCall(DMSwarmGetSize(dm, &Np));
62: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates));
63: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
64: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
65: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
66: PetscCall(VecViewNative(coordinates, viewer));
67: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
68: PetscCall(PetscViewerHDF5PopGroup(viewer));
69: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
72: #endif
74: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
75: {
76: DM dm;
77: #if defined(PETSC_HAVE_HDF5)
78: PetscBool ishdf5;
79: #endif
81: PetscFunctionBegin;
82: PetscCall(VecGetDM(v, &dm));
83: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
84: #if defined(PETSC_HAVE_HDF5)
85: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
86: if (ishdf5) {
87: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
90: #endif
91: PetscCall(VecViewNative(v, viewer));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: /*@
96: DMSwarmVectorGetField - Gets the field from which to define a `Vec` object
97: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
99: Not collective
101: Input Parameter:
102: . dm - a `DMSWARM`
104: Output Parameter:
105: . fieldname - the textual name given to a registered field, or the empty string if it has not been set
107: Level: beginner
109: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()` `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
110: @*/
111: PetscErrorCode DMSwarmVectorGetField(DM dm, const char *fieldname[])
112: {
113: PetscFunctionBegin;
115: PetscAssertPointer(fieldname, 2);
116: *fieldname = ((DM_Swarm *)dm->data)->vec_field_name;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
122: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
124: Collective
126: Input Parameters:
127: + dm - a `DMSWARM`
128: - fieldname - the textual name given to a registered field
130: Level: beginner
132: Notes:
133: The field with name `fieldname` must be defined as having a data type of `PetscScalar`.
135: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
136: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
138: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
139: @*/
140: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
141: {
142: DM_Swarm *swarm = (DM_Swarm *)dm->data;
143: PetscInt bs, n;
144: PetscScalar *array;
145: PetscDataType type;
147: PetscFunctionBegin;
149: if (fieldname) PetscAssertPointer(fieldname, 2);
150: if (!swarm->issetup) PetscCall(DMSetUp(dm));
151: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
152: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
154: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
155: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
156: PetscCall(PetscSNPrintf(swarm->vec_field_name, PETSC_MAX_PATH_LEN - 1, "%s", fieldname));
157: swarm->vec_field_set = PETSC_TRUE;
158: swarm->vec_field_bs = bs;
159: swarm->vec_field_nlocal = n;
160: PetscCall(DMSwarmRestoreField(dm, fieldname, &bs, &type, (void **)&array));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /* requires DMSwarmDefineFieldVector has been called */
165: static PetscErrorCode DMCreateGlobalVector_Swarm(DM dm, Vec *vec)
166: {
167: DM_Swarm *swarm = (DM_Swarm *)dm->data;
168: Vec x;
169: char name[PETSC_MAX_PATH_LEN];
171: PetscFunctionBegin;
172: if (!swarm->issetup) PetscCall(DMSetUp(dm));
173: PetscCheck(swarm->vec_field_set, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmVectorDefineField first");
174: PetscCheck(swarm->vec_field_nlocal == swarm->db->L, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
176: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmField_%s", swarm->vec_field_name));
177: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &x));
178: PetscCall(PetscObjectSetName((PetscObject)x, name));
179: PetscCall(VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE));
180: PetscCall(VecSetBlockSize(x, swarm->vec_field_bs));
181: PetscCall(VecSetDM(x, dm));
182: PetscCall(VecSetFromOptions(x));
183: PetscCall(VecSetDM(x, dm));
184: PetscCall(VecSetOperation(x, VECOP_VIEW, (void (*)(void))VecView_Swarm));
185: *vec = x;
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: /* requires DMSwarmDefineFieldVector has been called */
190: static PetscErrorCode DMCreateLocalVector_Swarm(DM dm, Vec *vec)
191: {
192: DM_Swarm *swarm = (DM_Swarm *)dm->data;
193: Vec x;
194: char name[PETSC_MAX_PATH_LEN];
196: PetscFunctionBegin;
197: if (!swarm->issetup) PetscCall(DMSetUp(dm));
198: PetscCheck(swarm->vec_field_set, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmVectorDefineField first");
199: PetscCheck(swarm->vec_field_nlocal == swarm->db->L, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSwarm sizes have changed since last call to VectorDefineField first");
201: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmField_%s", swarm->vec_field_name));
202: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
203: PetscCall(PetscObjectSetName((PetscObject)x, name));
204: PetscCall(VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE));
205: PetscCall(VecSetBlockSize(x, swarm->vec_field_bs));
206: PetscCall(VecSetDM(x, dm));
207: PetscCall(VecSetFromOptions(x));
208: *vec = x;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
213: {
214: DM_Swarm *swarm = (DM_Swarm *)dm->data;
215: DMSwarmDataField gfield;
216: PetscInt bs, nlocal, fid = -1, cfid = -2;
217: PetscBool flg;
219: PetscFunctionBegin;
220: /* check vector is an inplace array */
221: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
222: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
223: (void)flg; /* avoid compiler warning */
224: 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);
225: PetscCall(VecGetLocalSize(*vec, &nlocal));
226: PetscCall(VecGetBlockSize(*vec, &bs));
227: 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");
228: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
229: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
230: PetscCall(VecResetArray(*vec));
231: PetscCall(VecDestroy(vec));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
236: {
237: DM_Swarm *swarm = (DM_Swarm *)dm->data;
238: PetscDataType type;
239: PetscScalar *array;
240: PetscInt bs, n, fid;
241: char name[PETSC_MAX_PATH_LEN];
242: PetscMPIInt size;
243: PetscBool iscuda, iskokkos;
245: PetscFunctionBegin;
246: if (!swarm->issetup) PetscCall(DMSetUp(dm));
247: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
248: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
249: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
251: PetscCallMPI(MPI_Comm_size(comm, &size));
252: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
253: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
254: PetscCall(VecCreate(comm, vec));
255: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
256: PetscCall(VecSetBlockSize(*vec, bs));
257: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
258: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
259: else PetscCall(VecSetType(*vec, VECSTANDARD));
260: PetscCall(VecPlaceArray(*vec, array));
262: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
263: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
265: /* Set guard */
266: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
267: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
269: PetscCall(VecSetDM(*vec, dm));
270: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
276: \hat f = \sum_i f_i \phi_i
278: and a particle function is given by
280: f = \sum_i w_i \delta(x - x_i)
282: then we want to require that
284: M \hat f = M_p f
286: where the particle mass matrix is given by
288: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
290: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
291: his integral. We allow this with the boolean flag.
292: */
293: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
294: {
295: const char *name = "Mass Matrix";
296: MPI_Comm comm;
297: PetscDS prob;
298: PetscSection fsection, globalFSection;
299: PetscHSetIJ ht;
300: PetscLayout rLayout, colLayout;
301: PetscInt *dnz, *onz;
302: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
303: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
304: PetscScalar *elemMat;
305: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0, totNc = 0;
307: PetscFunctionBegin;
308: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
309: PetscCall(DMGetCoordinateDim(dmf, &dim));
310: PetscCall(DMGetDS(dmf, &prob));
311: PetscCall(PetscDSGetNumFields(prob, &Nf));
312: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
313: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
314: PetscCall(DMGetLocalSection(dmf, &fsection));
315: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
316: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
317: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
319: PetscCall(PetscLayoutCreate(comm, &colLayout));
320: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
321: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
322: PetscCall(PetscLayoutSetUp(colLayout));
323: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
324: PetscCall(PetscLayoutDestroy(&colLayout));
326: PetscCall(PetscLayoutCreate(comm, &rLayout));
327: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
328: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
329: PetscCall(PetscLayoutSetUp(rLayout));
330: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
331: PetscCall(PetscLayoutDestroy(&rLayout));
333: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
334: PetscCall(PetscHSetIJCreate(&ht));
336: PetscCall(PetscSynchronizedFlush(comm, NULL));
337: for (field = 0; field < Nf; ++field) {
338: PetscObject obj;
339: PetscClassId id;
340: PetscInt Nc;
342: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
343: PetscCall(PetscObjectGetClassId(obj, &id));
344: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
345: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
346: totNc += Nc;
347: }
348: /* count non-zeros */
349: PetscCall(DMSwarmSortGetAccess(dmc));
350: for (field = 0; field < Nf; ++field) {
351: PetscObject obj;
352: PetscClassId id;
353: PetscInt Nc;
355: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
356: PetscCall(PetscObjectGetClassId(obj, &id));
357: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
358: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
360: for (cell = cStart; cell < cEnd; ++cell) {
361: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
362: PetscInt numFIndices, numCIndices;
364: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
365: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
366: maxC = PetscMax(maxC, numCIndices);
367: {
368: PetscHashIJKey key;
369: PetscBool missing;
370: for (PetscInt i = 0; i < numFIndices; ++i) {
371: key.j = findices[i]; /* global column (from Plex) */
372: if (key.j >= 0) {
373: /* Get indices for coarse elements */
374: for (PetscInt j = 0; j < numCIndices; ++j) {
375: for (PetscInt c = 0; c < Nc; ++c) {
376: // TODO Need field offset on particle here
377: key.i = cindices[j] * totNc + c + rStart; /* global cols (from Swarm) */
378: if (key.i < 0) continue;
379: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
380: if (missing) {
381: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
382: else ++onz[key.i - rStart];
383: } else SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
384: }
385: }
386: }
387: }
388: PetscCall(PetscFree(cindices));
389: }
390: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
391: }
392: }
393: PetscCall(PetscHSetIJDestroy(&ht));
394: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
395: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
396: PetscCall(PetscFree2(dnz, onz));
397: PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
398: for (field = 0; field < Nf; ++field) {
399: PetscTabulation Tcoarse;
400: PetscObject obj;
401: PetscClassId id;
402: PetscReal *fieldVals;
403: PetscInt Nc;
405: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
406: PetscCall(PetscObjectGetClassId(obj, &id));
407: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
408: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
410: PetscCall(DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&fieldVals));
411: for (cell = cStart; cell < cEnd; ++cell) {
412: PetscInt *findices, *cindices;
413: PetscInt numFIndices, numCIndices;
415: /* TODO: Use DMField instead of assuming affine */
416: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
417: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
418: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
419: for (PetscInt j = 0; j < numCIndices; ++j) CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &fieldVals[cindices[j] * dim], &xi[j * dim]);
420: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
421: else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
422: /* Get elemMat entries by multiplying by weight */
423: PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
424: for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
425: for (PetscInt j = 0; j < numCIndices; ++j) {
426: for (PetscInt c = 0; c < Nc; ++c) {
427: // TODO Need field offset on particle and field here
428: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
429: elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
430: }
431: }
432: }
433: for (PetscInt j = 0; j < numCIndices; ++j)
434: // TODO Need field offset on particle here
435: for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
436: if (0) PetscCall(DMPrintCellMatrix(cell, name, numCIndices * Nc, numFIndices, elemMat));
437: PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
438: PetscCall(PetscFree(cindices));
439: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
440: PetscCall(PetscTabulationDestroy(&Tcoarse));
441: }
442: PetscCall(DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&fieldVals));
443: }
444: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
445: PetscCall(DMSwarmSortRestoreAccess(dmc));
446: PetscCall(PetscFree3(v0, J, invJ));
447: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
448: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /* Returns empty matrix for use with SNES FD */
453: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
454: {
455: Vec field;
456: PetscInt size;
458: PetscFunctionBegin;
459: PetscCall(DMGetGlobalVector(sw, &field));
460: PetscCall(VecGetLocalSize(field, &size));
461: PetscCall(DMRestoreGlobalVector(sw, &field));
462: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
463: PetscCall(MatSetFromOptions(*m));
464: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
465: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
466: PetscCall(MatZeroEntries(*m));
467: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
468: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
469: PetscCall(MatShift(*m, 1.0));
470: PetscCall(MatSetDM(*m, sw));
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: /* FEM cols, Particle rows */
475: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
476: {
477: PetscSection gsf;
478: PetscInt m, n, Np, bs = 1;
479: void *ctx;
480: PetscBool set = ((DM_Swarm *)dmCoarse->data)->vec_field_set;
481: char *name = ((DM_Swarm *)dmCoarse->data)->vec_field_name;
483: PetscFunctionBegin;
484: PetscCall(DMGetGlobalSection(dmFine, &gsf));
485: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
486: PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
487: // TODO Include all fields
488: if (set) PetscCall(DMSwarmGetFieldInfo(dmCoarse, name, &bs, NULL));
489: n = Np * bs;
490: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
491: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
492: PetscCall(MatSetType(*mass, dmCoarse->mattype));
493: PetscCall(DMGetApplicationContext(dmFine, &ctx));
495: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
496: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
501: {
502: const char *name = "Mass Matrix Square";
503: MPI_Comm comm;
504: PetscDS prob;
505: PetscSection fsection, globalFSection;
506: PetscHSetIJ ht;
507: PetscLayout rLayout, colLayout;
508: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
509: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
510: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
511: PetscScalar *elemMat, *elemMatSq;
512: PetscInt cdim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
514: PetscFunctionBegin;
515: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
516: PetscCall(DMGetCoordinateDim(dmf, &cdim));
517: PetscCall(DMGetDS(dmf, &prob));
518: PetscCall(PetscDSGetNumFields(prob, &Nf));
519: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
520: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
521: PetscCall(DMGetLocalSection(dmf, &fsection));
522: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
523: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
524: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
526: PetscCall(PetscLayoutCreate(comm, &colLayout));
527: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
528: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
529: PetscCall(PetscLayoutSetUp(colLayout));
530: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
531: PetscCall(PetscLayoutDestroy(&colLayout));
533: PetscCall(PetscLayoutCreate(comm, &rLayout));
534: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
535: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
536: PetscCall(PetscLayoutSetUp(rLayout));
537: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
538: PetscCall(PetscLayoutDestroy(&rLayout));
540: PetscCall(DMPlexGetDepth(dmf, &depth));
541: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
542: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
543: PetscCall(PetscMalloc1(maxAdjSize, &adj));
545: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
546: PetscCall(PetscHSetIJCreate(&ht));
547: /* Count nonzeros
548: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
549: */
550: PetscCall(DMSwarmSortGetAccess(dmc));
551: for (cell = cStart; cell < cEnd; ++cell) {
552: PetscInt i;
553: PetscInt *cindices;
554: PetscInt numCIndices;
555: #if 0
556: PetscInt adjSize = maxAdjSize, a, j;
557: #endif
559: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
560: maxC = PetscMax(maxC, numCIndices);
561: /* Diagonal block */
562: for (i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
563: #if 0
564: /* Off-diagonal blocks */
565: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
566: for (a = 0; a < adjSize; ++a) {
567: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
568: const PetscInt ncell = adj[a];
569: PetscInt *ncindices;
570: PetscInt numNCIndices;
572: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
573: {
574: PetscHashIJKey key;
575: PetscBool missing;
577: for (i = 0; i < numCIndices; ++i) {
578: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
579: if (key.i < 0) continue;
580: for (j = 0; j < numNCIndices; ++j) {
581: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
582: if (key.j < 0) continue;
583: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
584: if (missing) {
585: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
586: else ++onz[key.i - rStart];
587: }
588: }
589: }
590: }
591: PetscCall(PetscFree(ncindices));
592: }
593: }
594: #endif
595: PetscCall(PetscFree(cindices));
596: }
597: PetscCall(PetscHSetIJDestroy(&ht));
598: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
599: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
600: PetscCall(PetscFree2(dnz, onz));
601: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
602: /* Fill in values
603: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
604: Start just by producing block diagonal
605: Could loop over adjacent cells
606: Produce neighboring element matrix
607: TODO Determine which columns and rows correspond to shared dual vector
608: Do MatMatMult with rectangular matrices
609: Insert block
610: */
611: for (field = 0; field < Nf; ++field) {
612: PetscTabulation Tcoarse;
613: PetscObject obj;
614: PetscReal *coords;
615: PetscInt Nc, i;
617: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
618: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
619: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
620: PetscCall(DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
621: for (cell = cStart; cell < cEnd; ++cell) {
622: PetscInt *findices, *cindices;
623: PetscInt numFIndices, numCIndices;
624: PetscInt p, c;
626: /* TODO: Use DMField instead of assuming affine */
627: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
628: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
629: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
630: for (p = 0; p < numCIndices; ++p) CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p] * cdim], &xi[p * cdim]);
631: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
632: /* Get elemMat entries by multiplying by weight */
633: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
634: for (i = 0; i < numFIndices; ++i) {
635: for (p = 0; p < numCIndices; ++p) {
636: for (c = 0; c < Nc; ++c) {
637: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
638: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
639: }
640: }
641: }
642: PetscCall(PetscTabulationDestroy(&Tcoarse));
643: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
644: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
645: /* Block diagonal */
646: if (numCIndices) {
647: PetscBLASInt blasn, blask;
648: PetscScalar one = 1.0, zero = 0.0;
650: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
651: PetscCall(PetscBLASIntCast(numFIndices, &blask));
652: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
653: }
654: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
655: /* TODO off-diagonal */
656: PetscCall(PetscFree(cindices));
657: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
658: }
659: PetscCall(DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
660: }
661: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
662: PetscCall(PetscFree(adj));
663: PetscCall(DMSwarmSortRestoreAccess(dmc));
664: PetscCall(PetscFree3(v0, J, invJ));
665: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
666: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: /*@
671: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
673: Collective
675: Input Parameters:
676: + dmCoarse - a `DMSWARM`
677: - dmFine - a `DMPLEX`
679: Output Parameter:
680: . mass - the square of the particle mass matrix
682: Level: advanced
684: Note:
685: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
686: future to compute the full normal equations.
688: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
689: @*/
690: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
691: {
692: PetscInt n;
693: void *ctx;
695: PetscFunctionBegin;
696: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
697: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
698: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
699: PetscCall(MatSetType(*mass, dmCoarse->mattype));
700: PetscCall(DMGetApplicationContext(dmFine, &ctx));
702: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
703: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: /*@
708: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
710: Collective
712: Input Parameters:
713: + dm - a `DMSWARM`
714: - fieldname - the textual name given to a registered field
716: Output Parameter:
717: . vec - the vector
719: Level: beginner
721: Note:
722: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
724: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
725: @*/
726: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
727: {
728: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
730: PetscFunctionBegin;
732: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: /*@
737: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
739: Collective
741: Input Parameters:
742: + dm - a `DMSWARM`
743: - fieldname - the textual name given to a registered field
745: Output Parameter:
746: . vec - the vector
748: Level: beginner
750: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
751: @*/
752: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
753: {
754: PetscFunctionBegin;
756: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: /*@
761: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
763: Collective
765: Input Parameters:
766: + dm - a `DMSWARM`
767: - fieldname - the textual name given to a registered field
769: Output Parameter:
770: . vec - the vector
772: Level: beginner
774: Note:
775: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
777: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
778: @*/
779: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
780: {
781: MPI_Comm comm = PETSC_COMM_SELF;
783: PetscFunctionBegin;
784: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
785: PetscFunctionReturn(PETSC_SUCCESS);
786: }
788: /*@
789: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
791: Collective
793: Input Parameters:
794: + dm - a `DMSWARM`
795: - fieldname - the textual name given to a registered field
797: Output Parameter:
798: . vec - the vector
800: Level: beginner
802: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
803: @*/
804: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
805: {
806: PetscFunctionBegin;
808: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
809: PetscFunctionReturn(PETSC_SUCCESS);
810: }
812: /*@
813: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
815: Collective
817: Input Parameter:
818: . dm - a `DMSWARM`
820: Level: beginner
822: Note:
823: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
825: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
826: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
827: @*/
828: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
829: {
830: DM_Swarm *swarm = (DM_Swarm *)dm->data;
832: PetscFunctionBegin;
833: if (!swarm->field_registration_initialized) {
834: swarm->field_registration_initialized = PETSC_TRUE;
835: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
836: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
837: }
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: /*@
842: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
844: Collective
846: Input Parameter:
847: . dm - a `DMSWARM`
849: Level: beginner
851: Note:
852: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
854: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
855: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
856: @*/
857: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
858: {
859: DM_Swarm *swarm = (DM_Swarm *)dm->data;
861: PetscFunctionBegin;
862: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
863: swarm->field_registration_finalized = PETSC_TRUE;
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@
868: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
870: Not Collective
872: Input Parameters:
873: + dm - a `DMSWARM`
874: . nlocal - the length of each registered field
875: - buffer - the length of the buffer used to efficient dynamic re-sizing
877: Level: beginner
879: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
880: @*/
881: PetscErrorCode DMSwarmSetLocalSizes(DM dm, PetscInt nlocal, PetscInt buffer)
882: {
883: DM_Swarm *swarm = (DM_Swarm *)dm->data;
885: PetscFunctionBegin;
886: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
887: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
888: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: /*@
893: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
895: Collective
897: Input Parameters:
898: + dm - a `DMSWARM`
899: - dmcell - the `DM` to attach to the `DMSWARM`
901: Level: beginner
903: Note:
904: The attached `DM` (dmcell) will be queried for point location and
905: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
907: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
908: @*/
909: PetscErrorCode DMSwarmSetCellDM(DM dm, DM dmcell)
910: {
911: DM_Swarm *swarm = (DM_Swarm *)dm->data;
913: PetscFunctionBegin;
914: swarm->dmcell = dmcell;
915: PetscFunctionReturn(PETSC_SUCCESS);
916: }
918: /*@
919: DMSwarmGetCellDM - Fetches the attached cell `DM`
921: Collective
923: Input Parameter:
924: . dm - a `DMSWARM`
926: Output Parameter:
927: . dmcell - the `DM` which was attached to the `DMSWARM`
929: Level: beginner
931: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
932: @*/
933: PetscErrorCode DMSwarmGetCellDM(DM dm, DM *dmcell)
934: {
935: DM_Swarm *swarm = (DM_Swarm *)dm->data;
937: PetscFunctionBegin;
938: *dmcell = swarm->dmcell;
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }
942: /*@
943: DMSwarmGetLocalSize - Retrieves the local length of fields registered
945: Not Collective
947: Input Parameter:
948: . dm - a `DMSWARM`
950: Output Parameter:
951: . nlocal - the length of each registered field
953: Level: beginner
955: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
956: @*/
957: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
958: {
959: DM_Swarm *swarm = (DM_Swarm *)dm->data;
961: PetscFunctionBegin;
962: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: /*@
967: DMSwarmGetSize - Retrieves the total length of fields registered
969: Collective
971: Input Parameter:
972: . dm - a `DMSWARM`
974: Output Parameter:
975: . n - the total length of each registered field
977: Level: beginner
979: Note:
980: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
982: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
983: @*/
984: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
985: {
986: DM_Swarm *swarm = (DM_Swarm *)dm->data;
987: PetscInt nlocal;
989: PetscFunctionBegin;
990: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
991: PetscCallMPI(MPIU_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
992: PetscFunctionReturn(PETSC_SUCCESS);
993: }
995: /*@
996: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
998: Collective
1000: Input Parameters:
1001: + dm - a `DMSWARM`
1002: . fieldname - the textual name to identify this field
1003: . blocksize - the number of each data type
1004: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
1006: Level: beginner
1008: Notes:
1009: The textual name for each registered field must be unique.
1011: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1012: @*/
1013: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1014: {
1015: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1016: size_t size;
1018: PetscFunctionBegin;
1019: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1020: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1022: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1023: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1024: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1025: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1026: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1028: PetscCall(PetscDataTypeGetSize(type, &size));
1029: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1030: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1031: {
1032: DMSwarmDataField gfield;
1034: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1035: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1036: }
1037: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: /*@C
1042: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
1044: Collective
1046: Input Parameters:
1047: + dm - a `DMSWARM`
1048: . fieldname - the textual name to identify this field
1049: - size - the size in bytes of the user struct of each data type
1051: Level: beginner
1053: Note:
1054: The textual name for each registered field must be unique.
1056: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1057: @*/
1058: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1059: {
1060: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1062: PetscFunctionBegin;
1063: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1064: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1065: PetscFunctionReturn(PETSC_SUCCESS);
1066: }
1068: /*@C
1069: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1071: Collective
1073: Input Parameters:
1074: + dm - a `DMSWARM`
1075: . fieldname - the textual name to identify this field
1076: . size - the size in bytes of the user data type
1077: - blocksize - the number of each data type
1079: Level: beginner
1081: Note:
1082: The textual name for each registered field must be unique.
1084: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1085: @*/
1086: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1087: {
1088: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1090: PetscFunctionBegin;
1091: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1092: {
1093: DMSwarmDataField gfield;
1095: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1096: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1097: }
1098: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: /*@C
1103: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1105: Not Collective, No Fortran Support
1107: Input Parameters:
1108: + dm - a `DMSWARM`
1109: - fieldname - the textual name to identify this field
1111: Output Parameters:
1112: + blocksize - the number of each data type
1113: . type - the data type
1114: - data - pointer to raw array
1116: Level: beginner
1118: Notes:
1119: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1121: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1122: @*/
1123: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1124: {
1125: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1126: DMSwarmDataField gfield;
1128: PetscFunctionBegin;
1130: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1131: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1132: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1133: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1134: if (blocksize) *blocksize = gfield->bs;
1135: if (type) *type = gfield->petsc_type;
1136: PetscFunctionReturn(PETSC_SUCCESS);
1137: }
1139: /*@C
1140: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1142: Not Collective, No Fortran Support
1144: Input Parameters:
1145: + dm - a `DMSWARM`
1146: - fieldname - the textual name to identify this field
1148: Output Parameters:
1149: + blocksize - the number of each data type
1150: . type - the data type
1151: - data - pointer to raw array
1153: Level: beginner
1155: Notes:
1156: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1158: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1159: @*/
1160: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1161: {
1162: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1163: DMSwarmDataField gfield;
1165: PetscFunctionBegin;
1167: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1168: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1169: if (data) *data = NULL;
1170: PetscFunctionReturn(PETSC_SUCCESS);
1171: }
1173: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
1174: {
1175: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1176: DMSwarmDataField gfield;
1178: PetscFunctionBegin;
1180: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1181: if (blocksize) *blocksize = gfield->bs;
1182: if (type) *type = gfield->petsc_type;
1183: PetscFunctionReturn(PETSC_SUCCESS);
1184: }
1186: /*@
1187: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
1189: Not Collective
1191: Input Parameter:
1192: . dm - a `DMSWARM`
1194: Level: beginner
1196: Notes:
1197: The new point will have all fields initialized to zero.
1199: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
1200: @*/
1201: PetscErrorCode DMSwarmAddPoint(DM dm)
1202: {
1203: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1205: PetscFunctionBegin;
1206: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1207: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1208: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
1209: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: /*@
1214: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
1216: Not Collective
1218: Input Parameters:
1219: + dm - a `DMSWARM`
1220: - npoints - the number of new points to add
1222: Level: beginner
1224: Notes:
1225: The new point will have all fields initialized to zero.
1227: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
1228: @*/
1229: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1230: {
1231: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1232: PetscInt nlocal;
1234: PetscFunctionBegin;
1235: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1236: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1237: nlocal = nlocal + npoints;
1238: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
1239: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1240: PetscFunctionReturn(PETSC_SUCCESS);
1241: }
1243: /*@
1244: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
1246: Not Collective
1248: Input Parameter:
1249: . dm - a `DMSWARM`
1251: Level: beginner
1253: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
1254: @*/
1255: PetscErrorCode DMSwarmRemovePoint(DM dm)
1256: {
1257: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1259: PetscFunctionBegin;
1260: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1261: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
1262: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1263: PetscFunctionReturn(PETSC_SUCCESS);
1264: }
1266: /*@
1267: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
1269: Not Collective
1271: Input Parameters:
1272: + dm - a `DMSWARM`
1273: - idx - index of point to remove
1275: Level: beginner
1277: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1278: @*/
1279: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
1280: {
1281: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1283: PetscFunctionBegin;
1284: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1285: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
1286: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: /*@
1291: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
1293: Not Collective
1295: Input Parameters:
1296: + dm - a `DMSWARM`
1297: . pi - the index of the point to copy
1298: - pj - the point index where the copy should be located
1300: Level: beginner
1302: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1303: @*/
1304: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
1305: {
1306: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1308: PetscFunctionBegin;
1309: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1310: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
1311: PetscFunctionReturn(PETSC_SUCCESS);
1312: }
1314: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
1315: {
1316: PetscFunctionBegin;
1317: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: /*@
1322: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
1324: Collective
1326: Input Parameters:
1327: + dm - the `DMSWARM`
1328: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1330: Level: advanced
1332: Notes:
1333: The `DM` will be modified to accommodate received points.
1334: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
1335: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
1337: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
1338: @*/
1339: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
1340: {
1341: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1343: PetscFunctionBegin;
1344: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
1345: switch (swarm->migrate_type) {
1346: case DMSWARM_MIGRATE_BASIC:
1347: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
1348: break;
1349: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1350: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
1351: break;
1352: case DMSWARM_MIGRATE_DMCELLEXACT:
1353: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1354: case DMSWARM_MIGRATE_USER:
1355: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
1356: default:
1357: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
1358: }
1359: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
1360: PetscCall(DMClearGlobalVectors(dm));
1361: PetscFunctionReturn(PETSC_SUCCESS);
1362: }
1364: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
1366: /*
1367: DMSwarmCollectViewCreate
1369: * Applies a collection method and gathers point neighbour points into dm
1371: Notes:
1372: Users should call DMSwarmCollectViewDestroy() after
1373: they have finished computations associated with the collected points
1374: */
1376: /*@
1377: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1378: in neighbour ranks into the `DMSWARM`
1380: Collective
1382: Input Parameter:
1383: . dm - the `DMSWARM`
1385: Level: advanced
1387: Notes:
1388: Users should call `DMSwarmCollectViewDestroy()` after
1389: they have finished computations associated with the collected points
1391: Different collect methods are supported. See `DMSwarmSetCollectType()`.
1393: Developer Note:
1394: Create and Destroy routines create new objects that can get destroyed, they do not change the state
1395: of the current object.
1397: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
1398: @*/
1399: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1400: {
1401: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1402: PetscInt ng;
1404: PetscFunctionBegin;
1405: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
1406: PetscCall(DMSwarmGetLocalSize(dm, &ng));
1407: switch (swarm->collect_type) {
1408: case DMSWARM_COLLECT_BASIC:
1409: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
1410: break;
1411: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1412: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1413: case DMSWARM_COLLECT_GENERAL:
1414: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
1415: default:
1416: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
1417: }
1418: swarm->collect_view_active = PETSC_TRUE;
1419: swarm->collect_view_reset_nlocal = ng;
1420: PetscFunctionReturn(PETSC_SUCCESS);
1421: }
1423: /*@
1424: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
1426: Collective
1428: Input Parameters:
1429: . dm - the `DMSWARM`
1431: Notes:
1432: Users should call `DMSwarmCollectViewCreate()` before this function is called.
1434: Level: advanced
1436: Developer Note:
1437: Create and Destroy routines create new objects that can get destroyed, they do not change the state
1438: of the current object.
1440: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
1441: @*/
1442: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1443: {
1444: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1446: PetscFunctionBegin;
1447: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
1448: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
1449: swarm->collect_view_active = PETSC_FALSE;
1450: PetscFunctionReturn(PETSC_SUCCESS);
1451: }
1453: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
1454: {
1455: PetscInt dim;
1457: PetscFunctionBegin;
1458: PetscCall(DMSwarmSetNumSpecies(dm, 1));
1459: PetscCall(DMGetDimension(dm, &dim));
1460: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1461: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1462: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_coor, dim, PETSC_DOUBLE));
1463: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_cellid, 1, PETSC_INT));
1464: PetscFunctionReturn(PETSC_SUCCESS);
1465: }
1467: /*@
1468: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
1470: Collective
1472: Input Parameters:
1473: + dm - the `DMSWARM`
1474: - Npc - The number of particles per cell in the cell `DM`
1476: Level: intermediate
1478: Notes:
1479: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
1480: one particle is in each cell, it is placed at the centroid.
1482: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1483: @*/
1484: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1485: {
1486: DM cdm;
1487: PetscRandom rnd;
1488: DMPolytopeType ct;
1489: PetscBool simplex;
1490: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1491: PetscInt dim, d, cStart, cEnd, c, p;
1493: PetscFunctionBeginUser;
1494: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
1495: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
1496: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
1498: PetscCall(DMSwarmGetCellDM(dm, &cdm));
1499: PetscCall(DMGetDimension(cdm, &dim));
1500: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
1501: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
1502: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1504: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
1505: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1506: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1507: for (c = cStart; c < cEnd; ++c) {
1508: if (Npc == 1) {
1509: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
1510: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
1511: } else {
1512: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
1513: for (p = 0; p < Npc; ++p) {
1514: const PetscInt n = c * Npc + p;
1515: PetscReal sum = 0.0, refcoords[3];
1517: for (d = 0; d < dim; ++d) {
1518: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
1519: sum += refcoords[d];
1520: }
1521: if (simplex && sum > 0.0)
1522: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
1523: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
1524: }
1525: }
1526: }
1527: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1528: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
1529: PetscCall(PetscRandomDestroy(&rnd));
1530: PetscFunctionReturn(PETSC_SUCCESS);
1531: }
1533: /*@
1534: DMSwarmSetType - Set particular flavor of `DMSWARM`
1536: Collective
1538: Input Parameters:
1539: + dm - the `DMSWARM`
1540: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
1542: Level: advanced
1544: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
1545: @*/
1546: PetscErrorCode DMSwarmSetType(DM dm, DMSwarmType stype)
1547: {
1548: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1550: PetscFunctionBegin;
1551: swarm->swarm_type = stype;
1552: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(dm));
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: static PetscErrorCode DMSetup_Swarm(DM dm)
1557: {
1558: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1559: PetscMPIInt rank;
1560: PetscInt p, npoints, *rankval;
1562: PetscFunctionBegin;
1563: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
1564: swarm->issetup = PETSC_TRUE;
1566: if (swarm->swarm_type == DMSWARM_PIC) {
1567: /* check dmcell exists */
1568: PetscCheck(swarm->dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires you call DMSwarmSetCellDM");
1570: if (swarm->dmcell->ops->locatepointssubdomain) {
1571: /* check methods exists for exact ownership identificiation */
1572: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
1573: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1574: } else {
1575: /* check methods exist for point location AND rank neighbor identification */
1576: if (swarm->dmcell->ops->locatepoints) {
1577: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
1578: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1580: if (swarm->dmcell->ops->getneighbors) {
1581: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
1582: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1584: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1585: }
1586: }
1588: PetscCall(DMSwarmFinalizeFieldRegister(dm));
1590: /* check some fields were registered */
1591: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
1593: /* check local sizes were set */
1594: PetscCheck(swarm->db->L != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Local sizes must be set via DMSwarmSetLocalSizes()");
1596: /* initialize values in pid and rank placeholders */
1597: /* TODO: [pid - use MPI_Scan] */
1598: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1599: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
1600: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1601: for (p = 0; p < npoints; p++) rankval[p] = (PetscInt)rank;
1602: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1603: PetscFunctionReturn(PETSC_SUCCESS);
1604: }
1606: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1608: static PetscErrorCode DMDestroy_Swarm(DM dm)
1609: {
1610: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1612: PetscFunctionBegin;
1613: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
1614: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
1615: if (swarm->sort_context) PetscCall(DMSwarmSortDestroy(&swarm->sort_context));
1616: PetscCall(PetscFree(swarm));
1617: PetscFunctionReturn(PETSC_SUCCESS);
1618: }
1620: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1621: {
1622: DM cdm;
1623: PetscDraw draw;
1624: PetscReal *coords, oldPause, radius = 0.01;
1625: PetscInt Np, p, bs;
1627: PetscFunctionBegin;
1628: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
1629: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1630: PetscCall(DMSwarmGetCellDM(dm, &cdm));
1631: PetscCall(PetscDrawGetPause(draw, &oldPause));
1632: PetscCall(PetscDrawSetPause(draw, 0.0));
1633: PetscCall(DMView(cdm, viewer));
1634: PetscCall(PetscDrawSetPause(draw, oldPause));
1636: PetscCall(DMSwarmGetLocalSize(dm, &Np));
1637: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords));
1638: for (p = 0; p < Np; ++p) {
1639: const PetscInt i = p * bs;
1641: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
1642: }
1643: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords));
1644: PetscCall(PetscDrawFlush(draw));
1645: PetscCall(PetscDrawPause(draw));
1646: PetscCall(PetscDrawSave(draw));
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
1651: {
1652: PetscViewerFormat format;
1653: PetscInt *sizes;
1654: PetscInt dim, Np, maxSize = 17;
1655: MPI_Comm comm;
1656: PetscMPIInt rank, size, p;
1657: const char *name;
1659: PetscFunctionBegin;
1660: PetscCall(PetscViewerGetFormat(viewer, &format));
1661: PetscCall(DMGetDimension(dm, &dim));
1662: PetscCall(DMSwarmGetLocalSize(dm, &Np));
1663: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1664: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1665: PetscCallMPI(MPI_Comm_size(comm, &size));
1666: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1667: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
1668: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
1669: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
1670: else PetscCall(PetscCalloc1(3, &sizes));
1671: if (size < maxSize) {
1672: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
1673: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
1674: for (p = 0; p < size; ++p) {
1675: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
1676: }
1677: } else {
1678: PetscInt locMinMax[2] = {Np, Np};
1680: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
1681: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
1682: }
1683: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1684: PetscCall(PetscFree(sizes));
1685: if (format == PETSC_VIEWER_ASCII_INFO) {
1686: PetscInt *cell;
1688: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
1689: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1690: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell));
1691: for (p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%d: %" PetscInt_FMT "\n", p, cell[p]));
1692: PetscCall(PetscViewerFlush(viewer));
1693: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1694: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell));
1695: }
1696: PetscFunctionReturn(PETSC_SUCCESS);
1697: }
1699: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1700: {
1701: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1702: PetscBool iascii, ibinary, isvtk, isdraw;
1703: #if defined(PETSC_HAVE_HDF5)
1704: PetscBool ishdf5;
1705: #endif
1707: PetscFunctionBegin;
1710: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1711: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
1712: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1713: #if defined(PETSC_HAVE_HDF5)
1714: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1715: #endif
1716: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1717: if (iascii) {
1718: PetscViewerFormat format;
1720: PetscCall(PetscViewerGetFormat(viewer, &format));
1721: switch (format) {
1722: case PETSC_VIEWER_ASCII_INFO_DETAIL:
1723: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
1724: break;
1725: default:
1726: PetscCall(DMView_Swarm_Ascii(dm, viewer));
1727: }
1728: } else {
1729: #if defined(PETSC_HAVE_HDF5)
1730: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
1731: #endif
1732: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
1733: }
1734: PetscFunctionReturn(PETSC_SUCCESS);
1735: }
1737: /*@
1738: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
1739: 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.
1741: Noncollective
1743: Input Parameters:
1744: + sw - the `DMSWARM`
1745: . cellID - the integer id of the cell to be extracted and filtered
1746: - cellswarm - The `DMSWARM` to receive the cell
1748: Level: beginner
1750: Notes:
1751: This presently only supports `DMSWARM_PIC` type
1753: Should be restored with `DMSwarmRestoreCellSwarm()`
1755: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
1757: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
1758: @*/
1759: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1760: {
1761: DM_Swarm *original = (DM_Swarm *)sw->data;
1762: DMLabel label;
1763: DM dmc, subdmc;
1764: PetscInt *pids, particles, dim;
1766: PetscFunctionBegin;
1767: /* Configure new swarm */
1768: PetscCall(DMSetType(cellswarm, DMSWARM));
1769: PetscCall(DMGetDimension(sw, &dim));
1770: PetscCall(DMSetDimension(cellswarm, dim));
1771: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
1772: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
1773: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
1774: PetscCall(DMSwarmSortGetAccess(sw));
1775: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
1776: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
1777: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
1778: PetscCall(DMSwarmSortRestoreAccess(sw));
1779: PetscCall(PetscFree(pids));
1780: PetscCall(DMSwarmGetCellDM(sw, &dmc));
1781: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
1782: PetscCall(DMAddLabel(dmc, label));
1783: PetscCall(DMLabelSetValue(label, cellID, 1));
1784: PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdmc));
1785: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
1786: PetscCall(DMLabelDestroy(&label));
1787: PetscFunctionReturn(PETSC_SUCCESS);
1788: }
1790: /*@
1791: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
1793: Noncollective
1795: Input Parameters:
1796: + sw - the parent `DMSWARM`
1797: . cellID - the integer id of the cell to be copied back into the parent swarm
1798: - cellswarm - the cell swarm object
1800: Level: beginner
1802: Note:
1803: This only supports `DMSWARM_PIC` types of `DMSWARM`s
1805: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
1806: @*/
1807: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1808: {
1809: DM dmc;
1810: PetscInt *pids, particles, p;
1812: PetscFunctionBegin;
1813: PetscCall(DMSwarmSortGetAccess(sw));
1814: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
1815: PetscCall(DMSwarmSortRestoreAccess(sw));
1816: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
1817: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
1818: /* Free memory, destroy cell dm */
1819: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
1820: PetscCall(DMDestroy(&dmc));
1821: PetscCall(PetscFree(pids));
1822: PetscFunctionReturn(PETSC_SUCCESS);
1823: }
1825: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
1827: static PetscErrorCode DMInitialize_Swarm(DM sw)
1828: {
1829: PetscFunctionBegin;
1830: sw->dim = 0;
1831: sw->ops->view = DMView_Swarm;
1832: sw->ops->load = NULL;
1833: sw->ops->setfromoptions = NULL;
1834: sw->ops->clone = DMClone_Swarm;
1835: sw->ops->setup = DMSetup_Swarm;
1836: sw->ops->createlocalsection = NULL;
1837: sw->ops->createsectionpermutation = NULL;
1838: sw->ops->createdefaultconstraints = NULL;
1839: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1840: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
1841: sw->ops->getlocaltoglobalmapping = NULL;
1842: sw->ops->createfieldis = NULL;
1843: sw->ops->createcoordinatedm = NULL;
1844: sw->ops->getcoloring = NULL;
1845: sw->ops->creatematrix = DMCreateMatrix_Swarm;
1846: sw->ops->createinterpolation = NULL;
1847: sw->ops->createinjection = NULL;
1848: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1849: sw->ops->refine = NULL;
1850: sw->ops->coarsen = NULL;
1851: sw->ops->refinehierarchy = NULL;
1852: sw->ops->coarsenhierarchy = NULL;
1853: sw->ops->globaltolocalbegin = NULL;
1854: sw->ops->globaltolocalend = NULL;
1855: sw->ops->localtoglobalbegin = NULL;
1856: sw->ops->localtoglobalend = NULL;
1857: sw->ops->destroy = DMDestroy_Swarm;
1858: sw->ops->createsubdm = NULL;
1859: sw->ops->getdimpoints = NULL;
1860: sw->ops->locatepoints = NULL;
1861: PetscFunctionReturn(PETSC_SUCCESS);
1862: }
1864: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
1865: {
1866: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1868: PetscFunctionBegin;
1869: swarm->refct++;
1870: (*newdm)->data = swarm;
1871: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
1872: PetscCall(DMInitialize_Swarm(*newdm));
1873: (*newdm)->dim = dm->dim;
1874: PetscFunctionReturn(PETSC_SUCCESS);
1875: }
1877: /*MC
1878: DMSWARM = "swarm" - A `DM` object used to represent arrays of data (fields) of arbitrary data type.
1879: This implementation was designed for particle methods in which the underlying
1880: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1882: Level: intermediate
1884: Notes:
1885: User data can be represented by `DMSWARM` through a registering "fields".
1886: To register a field, the user must provide;
1887: (a) a unique name;
1888: (b) the data type (or size in bytes);
1889: (c) the block size of the data.
1891: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1892: on a set of particles. Then the following code could be used
1893: .vb
1894: DMSwarmInitializeFieldRegister(dm)
1895: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1896: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1897: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1898: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1899: DMSwarmFinalizeFieldRegister(dm)
1900: .ve
1902: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
1903: The only restriction imposed by `DMSWARM` is that all fields contain the same number of points.
1905: To support particle methods, "migration" techniques are provided. These methods migrate data
1906: between ranks.
1908: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
1909: As a `DMSWARM` may internally define and store values of different data types,
1910: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
1911: fields should be used to define a `Vec` object via
1912: `DMSwarmVectorDefineField()`
1913: The specified field can be changed at any time - thereby permitting vectors
1914: compatible with different fields to be created.
1916: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via
1917: `DMSwarmCreateGlobalVectorFromField()`
1918: Here the data defining the field in the `DMSWARM` is shared with a Vec.
1919: This is inherently unsafe if you alter the size of the field at any time between
1920: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
1921: If the local size of the `DMSWARM` does not match the local size of the global vector
1922: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
1924: Additional high-level support is provided for Particle-In-Cell methods.
1925: Please refer to `DMSwarmSetType()`.
1927: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`
1928: M*/
1930: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1931: {
1932: DM_Swarm *swarm;
1934: PetscFunctionBegin;
1936: PetscCall(PetscNew(&swarm));
1937: dm->data = swarm;
1938: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
1939: PetscCall(DMSwarmInitializeFieldRegister(dm));
1940: swarm->refct = 1;
1941: swarm->vec_field_set = PETSC_FALSE;
1942: swarm->issetup = PETSC_FALSE;
1943: swarm->swarm_type = DMSWARM_BASIC;
1944: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1945: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1946: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1947: swarm->dmcell = NULL;
1948: swarm->collect_view_active = PETSC_FALSE;
1949: swarm->collect_view_reset_nlocal = -1;
1950: PetscCall(DMInitialize_Swarm(dm));
1951: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
1952: PetscFunctionReturn(PETSC_SUCCESS);
1953: }