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;
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(PetscViewerHDF5SetTimestep(viewer, seqnum));
46: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
47: PetscCall(VecViewNative(v, viewer));
48: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
49: PetscCall(PetscViewerHDF5PopGroup(viewer));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
54: {
55: Vec coordinates;
56: PetscInt Np;
57: PetscBool isseq;
59: PetscFunctionBegin;
60: PetscCall(DMSwarmGetSize(dm, &Np));
61: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates));
62: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
63: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
64: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
65: PetscCall(VecViewNative(coordinates, viewer));
66: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
67: PetscCall(PetscViewerHDF5PopGroup(viewer));
68: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates));
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
71: #endif
73: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
74: {
75: DM dm;
76: #if defined(PETSC_HAVE_HDF5)
77: PetscBool ishdf5;
78: #endif
80: PetscFunctionBegin;
81: PetscCall(VecGetDM(v, &dm));
82: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
83: #if defined(PETSC_HAVE_HDF5)
84: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
85: if (ishdf5) {
86: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
89: #endif
90: PetscCall(VecViewNative(v, viewer));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /*@C
95: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
96: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
98: Collective
100: Input Parameters:
101: + dm - a `DMSWARM`
102: - fieldname - the textual name given to a registered field
104: Level: beginner
106: Notes:
107: The field with name `fieldname` must be defined as having a data type of `PetscScalar`.
109: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
110: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
112: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
113: @*/
114: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
115: {
116: DM_Swarm *swarm = (DM_Swarm *)dm->data;
117: PetscInt bs, n;
118: PetscScalar *array;
119: PetscDataType type;
121: PetscFunctionBegin;
122: if (!swarm->issetup) PetscCall(DMSetUp(dm));
123: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
124: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
126: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
127: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
128: PetscCall(PetscSNPrintf(swarm->vec_field_name, PETSC_MAX_PATH_LEN - 1, "%s", fieldname));
129: swarm->vec_field_set = PETSC_TRUE;
130: swarm->vec_field_bs = bs;
131: swarm->vec_field_nlocal = n;
132: PetscCall(DMSwarmRestoreField(dm, fieldname, &bs, &type, (void **)&array));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: /* requires DMSwarmDefineFieldVector has been called */
137: static PetscErrorCode DMCreateGlobalVector_Swarm(DM dm, Vec *vec)
138: {
139: DM_Swarm *swarm = (DM_Swarm *)dm->data;
140: Vec x;
141: char name[PETSC_MAX_PATH_LEN];
143: PetscFunctionBegin;
144: if (!swarm->issetup) PetscCall(DMSetUp(dm));
145: PetscCheck(swarm->vec_field_set, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmVectorDefineField first");
146: 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 */
148: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmField_%s", swarm->vec_field_name));
149: PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &x));
150: PetscCall(PetscObjectSetName((PetscObject)x, name));
151: PetscCall(VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE));
152: PetscCall(VecSetBlockSize(x, swarm->vec_field_bs));
153: PetscCall(VecSetDM(x, dm));
154: PetscCall(VecSetFromOptions(x));
155: PetscCall(VecSetDM(x, dm));
156: PetscCall(VecSetOperation(x, VECOP_VIEW, (void (*)(void))VecView_Swarm));
157: *vec = x;
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: /* requires DMSwarmDefineFieldVector has been called */
162: static PetscErrorCode DMCreateLocalVector_Swarm(DM dm, Vec *vec)
163: {
164: DM_Swarm *swarm = (DM_Swarm *)dm->data;
165: Vec x;
166: char name[PETSC_MAX_PATH_LEN];
168: PetscFunctionBegin;
169: if (!swarm->issetup) PetscCall(DMSetUp(dm));
170: PetscCheck(swarm->vec_field_set, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmVectorDefineField first");
171: PetscCheck(swarm->vec_field_nlocal == swarm->db->L, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSwarm sizes have changed since last call to VectorDefineField first");
173: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmField_%s", swarm->vec_field_name));
174: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
175: PetscCall(PetscObjectSetName((PetscObject)x, name));
176: PetscCall(VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE));
177: PetscCall(VecSetBlockSize(x, swarm->vec_field_bs));
178: PetscCall(VecSetDM(x, dm));
179: PetscCall(VecSetFromOptions(x));
180: *vec = x;
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
185: {
186: DM_Swarm *swarm = (DM_Swarm *)dm->data;
187: DMSwarmDataField gfield;
188: PetscInt bs, nlocal, fid = -1, cfid = -2;
189: PetscBool flg;
191: PetscFunctionBegin;
192: /* check vector is an inplace array */
193: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
194: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
195: 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);
196: PetscCall(VecGetLocalSize(*vec, &nlocal));
197: PetscCall(VecGetBlockSize(*vec, &bs));
198: 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");
199: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
200: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
201: PetscCall(VecResetArray(*vec));
202: PetscCall(VecDestroy(vec));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
207: {
208: DM_Swarm *swarm = (DM_Swarm *)dm->data;
209: PetscDataType type;
210: PetscScalar *array;
211: PetscInt bs, n, fid;
212: char name[PETSC_MAX_PATH_LEN];
213: PetscMPIInt size;
214: PetscBool iscuda, iskokkos;
216: PetscFunctionBegin;
217: if (!swarm->issetup) PetscCall(DMSetUp(dm));
218: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
219: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
220: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
222: PetscCallMPI(MPI_Comm_size(comm, &size));
223: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
224: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
225: PetscCall(VecCreate(comm, vec));
226: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
227: PetscCall(VecSetBlockSize(*vec, bs));
228: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
229: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
230: else PetscCall(VecSetType(*vec, VECSTANDARD));
231: PetscCall(VecPlaceArray(*vec, array));
233: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
234: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
236: /* Set guard */
237: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
238: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
240: PetscCall(VecSetDM(*vec, dm));
241: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
247: \hat f = \sum_i f_i \phi_i
249: and a particle function is given by
251: f = \sum_i w_i \delta(x - x_i)
253: then we want to require that
255: M \hat f = M_p f
257: where the particle mass matrix is given by
259: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
261: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
262: his integral. We allow this with the boolean flag.
263: */
264: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
265: {
266: const char *name = "Mass Matrix";
267: MPI_Comm comm;
268: PetscDS prob;
269: PetscSection fsection, globalFSection;
270: PetscHSetIJ ht;
271: PetscLayout rLayout, colLayout;
272: PetscInt *dnz, *onz;
273: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
274: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
275: PetscScalar *elemMat;
276: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
278: PetscFunctionBegin;
279: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
280: PetscCall(DMGetCoordinateDim(dmf, &dim));
281: PetscCall(DMGetDS(dmf, &prob));
282: PetscCall(PetscDSGetNumFields(prob, &Nf));
283: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
284: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
285: PetscCall(DMGetLocalSection(dmf, &fsection));
286: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
287: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
288: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
290: PetscCall(PetscLayoutCreate(comm, &colLayout));
291: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
292: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
293: PetscCall(PetscLayoutSetUp(colLayout));
294: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
295: PetscCall(PetscLayoutDestroy(&colLayout));
297: PetscCall(PetscLayoutCreate(comm, &rLayout));
298: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
299: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
300: PetscCall(PetscLayoutSetUp(rLayout));
301: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
302: PetscCall(PetscLayoutDestroy(&rLayout));
304: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
305: PetscCall(PetscHSetIJCreate(&ht));
307: PetscCall(PetscSynchronizedFlush(comm, NULL));
308: /* count non-zeros */
309: PetscCall(DMSwarmSortGetAccess(dmc));
310: for (field = 0; field < Nf; ++field) {
311: for (cell = cStart; cell < cEnd; ++cell) {
312: PetscInt c, i;
313: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
314: PetscInt numFIndices, numCIndices;
316: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
317: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
318: maxC = PetscMax(maxC, numCIndices);
319: {
320: PetscHashIJKey key;
321: PetscBool missing;
322: for (i = 0; i < numFIndices; ++i) {
323: key.j = findices[i]; /* global column (from Plex) */
324: if (key.j >= 0) {
325: /* Get indices for coarse elements */
326: for (c = 0; c < numCIndices; ++c) {
327: key.i = cindices[c] + rStart; /* global cols (from Swarm) */
328: if (key.i < 0) continue;
329: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
330: if (missing) {
331: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
332: else ++onz[key.i - rStart];
333: } else SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
334: }
335: }
336: }
337: PetscCall(PetscFree(cindices));
338: }
339: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
340: }
341: }
342: PetscCall(PetscHSetIJDestroy(&ht));
343: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
344: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
345: PetscCall(PetscFree2(dnz, onz));
346: PetscCall(PetscMalloc3(maxC * totDim, &elemMat, maxC, &rowIDXs, maxC * dim, &xi));
347: for (field = 0; field < Nf; ++field) {
348: PetscTabulation Tcoarse;
349: PetscObject obj;
350: PetscReal *fieldVals;
351: PetscInt Nc, i;
353: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
354: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
355: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
357: PetscCall(DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&fieldVals));
358: for (cell = cStart; cell < cEnd; ++cell) {
359: PetscInt *findices, *cindices;
360: PetscInt numFIndices, numCIndices;
361: PetscInt p, c;
363: /* TODO: Use DMField instead of assuming affine */
364: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
365: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
366: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
367: for (p = 0; p < numCIndices; ++p) CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &fieldVals[cindices[p] * dim], &xi[p * dim]);
368: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
369: /* Get elemMat entries by multiplying by weight */
370: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
371: for (i = 0; i < numFIndices; ++i) {
372: for (p = 0; p < numCIndices; ++p) {
373: for (c = 0; c < Nc; ++c) {
374: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
375: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
376: }
377: }
378: }
379: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
380: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
381: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
382: PetscCall(PetscFree(cindices));
383: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
384: PetscCall(PetscTabulationDestroy(&Tcoarse));
385: }
386: PetscCall(DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&fieldVals));
387: }
388: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
389: PetscCall(DMSwarmSortRestoreAccess(dmc));
390: PetscCall(PetscFree3(v0, J, invJ));
391: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
392: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /* Returns empty matrix for use with SNES FD */
397: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
398: {
399: Vec field;
400: PetscInt size;
402: PetscFunctionBegin;
403: PetscCall(DMGetGlobalVector(sw, &field));
404: PetscCall(VecGetLocalSize(field, &size));
405: PetscCall(DMRestoreGlobalVector(sw, &field));
406: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
407: PetscCall(MatSetFromOptions(*m));
408: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
409: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
410: PetscCall(MatZeroEntries(*m));
411: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
412: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
413: PetscCall(MatShift(*m, 1.0));
414: PetscCall(MatSetDM(*m, sw));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: /* FEM cols, Particle rows */
419: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
420: {
421: PetscSection gsf;
422: PetscInt m, n;
423: void *ctx;
425: PetscFunctionBegin;
426: PetscCall(DMGetGlobalSection(dmFine, &gsf));
427: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
428: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
429: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
430: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
431: PetscCall(MatSetType(*mass, dmCoarse->mattype));
432: PetscCall(DMGetApplicationContext(dmFine, &ctx));
434: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
435: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
440: {
441: const char *name = "Mass Matrix Square";
442: MPI_Comm comm;
443: PetscDS prob;
444: PetscSection fsection, globalFSection;
445: PetscHSetIJ ht;
446: PetscLayout rLayout, colLayout;
447: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
448: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
449: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
450: PetscScalar *elemMat, *elemMatSq;
451: PetscInt cdim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
453: PetscFunctionBegin;
454: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
455: PetscCall(DMGetCoordinateDim(dmf, &cdim));
456: PetscCall(DMGetDS(dmf, &prob));
457: PetscCall(PetscDSGetNumFields(prob, &Nf));
458: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
459: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
460: PetscCall(DMGetLocalSection(dmf, &fsection));
461: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
462: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
463: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
465: PetscCall(PetscLayoutCreate(comm, &colLayout));
466: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
467: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
468: PetscCall(PetscLayoutSetUp(colLayout));
469: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
470: PetscCall(PetscLayoutDestroy(&colLayout));
472: PetscCall(PetscLayoutCreate(comm, &rLayout));
473: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
474: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
475: PetscCall(PetscLayoutSetUp(rLayout));
476: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
477: PetscCall(PetscLayoutDestroy(&rLayout));
479: PetscCall(DMPlexGetDepth(dmf, &depth));
480: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
481: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
482: PetscCall(PetscMalloc1(maxAdjSize, &adj));
484: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
485: PetscCall(PetscHSetIJCreate(&ht));
486: /* Count nonzeros
487: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
488: */
489: PetscCall(DMSwarmSortGetAccess(dmc));
490: for (cell = cStart; cell < cEnd; ++cell) {
491: PetscInt i;
492: PetscInt *cindices;
493: PetscInt numCIndices;
494: #if 0
495: PetscInt adjSize = maxAdjSize, a, j;
496: #endif
498: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
499: maxC = PetscMax(maxC, numCIndices);
500: /* Diagonal block */
501: for (i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
502: #if 0
503: /* Off-diagonal blocks */
504: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
505: for (a = 0; a < adjSize; ++a) {
506: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
507: const PetscInt ncell = adj[a];
508: PetscInt *ncindices;
509: PetscInt numNCIndices;
511: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
512: {
513: PetscHashIJKey key;
514: PetscBool missing;
516: for (i = 0; i < numCIndices; ++i) {
517: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
518: if (key.i < 0) continue;
519: for (j = 0; j < numNCIndices; ++j) {
520: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
521: if (key.j < 0) continue;
522: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
523: if (missing) {
524: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
525: else ++onz[key.i - rStart];
526: }
527: }
528: }
529: }
530: PetscCall(PetscFree(ncindices));
531: }
532: }
533: #endif
534: PetscCall(PetscFree(cindices));
535: }
536: PetscCall(PetscHSetIJDestroy(&ht));
537: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
538: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
539: PetscCall(PetscFree2(dnz, onz));
540: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
541: /* Fill in values
542: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
543: Start just by producing block diagonal
544: Could loop over adjacent cells
545: Produce neighboring element matrix
546: TODO Determine which columns and rows correspond to shared dual vector
547: Do MatMatMult with rectangular matrices
548: Insert block
549: */
550: for (field = 0; field < Nf; ++field) {
551: PetscTabulation Tcoarse;
552: PetscObject obj;
553: PetscReal *coords;
554: PetscInt Nc, i;
556: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
557: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
558: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
559: PetscCall(DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
560: for (cell = cStart; cell < cEnd; ++cell) {
561: PetscInt *findices, *cindices;
562: PetscInt numFIndices, numCIndices;
563: PetscInt p, c;
565: /* TODO: Use DMField instead of assuming affine */
566: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
567: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
568: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
569: for (p = 0; p < numCIndices; ++p) CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p] * cdim], &xi[p * cdim]);
570: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
571: /* Get elemMat entries by multiplying by weight */
572: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
573: for (i = 0; i < numFIndices; ++i) {
574: for (p = 0; p < numCIndices; ++p) {
575: for (c = 0; c < Nc; ++c) {
576: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
577: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
578: }
579: }
580: }
581: PetscCall(PetscTabulationDestroy(&Tcoarse));
582: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
583: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
584: /* Block diagonal */
585: if (numCIndices) {
586: PetscBLASInt blasn, blask;
587: PetscScalar one = 1.0, zero = 0.0;
589: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
590: PetscCall(PetscBLASIntCast(numFIndices, &blask));
591: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
592: }
593: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
594: /* TODO off-diagonal */
595: PetscCall(PetscFree(cindices));
596: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
597: }
598: PetscCall(DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
599: }
600: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
601: PetscCall(PetscFree(adj));
602: PetscCall(DMSwarmSortRestoreAccess(dmc));
603: PetscCall(PetscFree3(v0, J, invJ));
604: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
605: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /*@C
610: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
612: Collective
614: Input Parameters:
615: + dmCoarse - a `DMSWARM`
616: - dmFine - a `DMPLEX`
618: Output Parameter:
619: . mass - the square of the particle mass matrix
621: Level: advanced
623: Note:
624: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
625: future to compute the full normal equations.
627: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
628: @*/
629: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
630: {
631: PetscInt n;
632: void *ctx;
634: PetscFunctionBegin;
635: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
636: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
637: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
638: PetscCall(MatSetType(*mass, dmCoarse->mattype));
639: PetscCall(DMGetApplicationContext(dmFine, &ctx));
641: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
642: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: /*@C
647: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
649: Collective
651: Input Parameters:
652: + dm - a `DMSWARM`
653: - fieldname - the textual name given to a registered field
655: Output Parameter:
656: . vec - the vector
658: Level: beginner
660: Notes:
661: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
663: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
664: @*/
665: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
666: {
667: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
669: PetscFunctionBegin;
671: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: /*@C
676: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
678: Collective
680: Input Parameters:
681: + dm - a `DMSWARM`
682: - fieldname - the textual name given to a registered field
684: Output Parameter:
685: . vec - the vector
687: Level: beginner
689: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
690: @*/
691: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
692: {
693: PetscFunctionBegin;
695: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: /*@C
700: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
702: Collective
704: Input Parameters:
705: + dm - a `DMSWARM`
706: - fieldname - the textual name given to a registered field
708: Output Parameter:
709: . vec - the vector
711: Level: beginner
713: Note:
714: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
716: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
717: @*/
718: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
719: {
720: MPI_Comm comm = PETSC_COMM_SELF;
722: PetscFunctionBegin;
723: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: /*@C
728: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
730: Collective
732: Input Parameters:
733: + dm - a `DMSWARM`
734: - fieldname - the textual name given to a registered field
736: Output Parameter:
737: . vec - the vector
739: Level: beginner
741: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
742: @*/
743: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
744: {
745: PetscFunctionBegin;
747: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: /*@C
752: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
754: Collective
756: Input Parameter:
757: . dm - a `DMSWARM`
759: Level: beginner
761: Note:
762: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
764: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
765: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
766: @*/
767: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
768: {
769: DM_Swarm *swarm = (DM_Swarm *)dm->data;
771: PetscFunctionBegin;
772: if (!swarm->field_registration_initialized) {
773: swarm->field_registration_initialized = PETSC_TRUE;
774: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
775: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
776: }
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: /*@
781: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
783: Collective
785: Input Parameter:
786: . dm - a `DMSWARM`
788: Level: beginner
790: Note:
791: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
793: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
794: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
795: @*/
796: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
797: {
798: DM_Swarm *swarm = (DM_Swarm *)dm->data;
800: PetscFunctionBegin;
801: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
802: swarm->field_registration_finalized = PETSC_TRUE;
803: PetscFunctionReturn(PETSC_SUCCESS);
804: }
806: /*@
807: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
809: Not Collective
811: Input Parameters:
812: + dm - a `DMSWARM`
813: . nlocal - the length of each registered field
814: - buffer - the length of the buffer used to efficient dynamic re-sizing
816: Level: beginner
818: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
819: @*/
820: PetscErrorCode DMSwarmSetLocalSizes(DM dm, PetscInt nlocal, PetscInt buffer)
821: {
822: DM_Swarm *swarm = (DM_Swarm *)dm->data;
824: PetscFunctionBegin;
825: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
826: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
827: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: /*@
832: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
834: Collective
836: Input Parameters:
837: + dm - a `DMSWARM`
838: - dmcell - the `DM` to attach to the `DMSWARM`
840: Level: beginner
842: Note:
843: The attached `DM` (dmcell) will be queried for point location and
844: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
846: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
847: @*/
848: PetscErrorCode DMSwarmSetCellDM(DM dm, DM dmcell)
849: {
850: DM_Swarm *swarm = (DM_Swarm *)dm->data;
852: PetscFunctionBegin;
853: swarm->dmcell = dmcell;
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: /*@
858: DMSwarmGetCellDM - Fetches the attached cell `DM`
860: Collective
862: Input Parameter:
863: . dm - a `DMSWARM`
865: Output Parameter:
866: . dmcell - the `DM` which was attached to the `DMSWARM`
868: Level: beginner
870: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
871: @*/
872: PetscErrorCode DMSwarmGetCellDM(DM dm, DM *dmcell)
873: {
874: DM_Swarm *swarm = (DM_Swarm *)dm->data;
876: PetscFunctionBegin;
877: *dmcell = swarm->dmcell;
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: /*@
882: DMSwarmGetLocalSize - Retrieves the local length of fields registered
884: Not Collective
886: Input Parameter:
887: . dm - a `DMSWARM`
889: Output Parameter:
890: . nlocal - the length of each registered field
892: Level: beginner
894: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
895: @*/
896: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
897: {
898: DM_Swarm *swarm = (DM_Swarm *)dm->data;
900: PetscFunctionBegin;
901: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
902: PetscFunctionReturn(PETSC_SUCCESS);
903: }
905: /*@
906: DMSwarmGetSize - Retrieves the total length of fields registered
908: Collective
910: Input Parameter:
911: . dm - a `DMSWARM`
913: Output Parameter:
914: . n - the total length of each registered field
916: Level: beginner
918: Note:
919: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
921: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
922: @*/
923: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
924: {
925: DM_Swarm *swarm = (DM_Swarm *)dm->data;
926: PetscInt nlocal;
928: PetscFunctionBegin;
929: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
930: PetscCall(MPIU_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: /*@C
935: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
937: Collective
939: Input Parameters:
940: + dm - a `DMSWARM`
941: . fieldname - the textual name to identify this field
942: . blocksize - the number of each data type
943: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
945: Level: beginner
947: Notes:
948: The textual name for each registered field must be unique.
950: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
951: @*/
952: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
953: {
954: DM_Swarm *swarm = (DM_Swarm *)dm->data;
955: size_t size;
957: PetscFunctionBegin;
958: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
959: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
961: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
962: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
963: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
964: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
965: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
967: PetscCall(PetscDataTypeGetSize(type, &size));
968: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
969: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
970: {
971: DMSwarmDataField gfield;
973: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
974: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
975: }
976: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /*@C
981: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
983: Collective
985: Input Parameters:
986: + dm - a `DMSWARM`
987: . fieldname - the textual name to identify this field
988: - size - the size in bytes of the user struct of each data type
990: Level: beginner
992: Note:
993: The textual name for each registered field must be unique.
995: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
996: @*/
997: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
998: {
999: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1001: PetscFunctionBegin;
1002: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1003: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: /*@C
1008: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1010: Collective
1012: Input Parameters:
1013: + dm - a `DMSWARM`
1014: . fieldname - the textual name to identify this field
1015: . size - the size in bytes of the user data type
1016: - blocksize - the number of each data type
1018: Level: beginner
1020: Note:
1021: The textual name for each registered field must be unique.
1023: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1024: @*/
1025: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1026: {
1027: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1029: PetscFunctionBegin;
1030: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", 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 = PETSC_DATATYPE_UNKNOWN;
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: /*@C
1042: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1044: Not Collective
1046: Input Parameters:
1047: + dm - a `DMSWARM`
1048: - fieldname - the textual name to identify this field
1050: Output Parameters:
1051: + blocksize - the number of each data type
1052: . type - the data type
1053: - data - pointer to raw array
1055: Level: beginner
1057: Notes:
1058: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1060: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1061: @*/
1062: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1063: {
1064: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1065: DMSwarmDataField gfield;
1067: PetscFunctionBegin;
1069: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1070: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1071: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1072: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1073: if (blocksize) *blocksize = gfield->bs;
1074: if (type) *type = gfield->petsc_type;
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: /*@C
1079: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1081: Not Collective
1083: Input Parameters:
1084: + dm - a `DMSWARM`
1085: - fieldname - the textual name to identify this field
1087: Output Parameters:
1088: + blocksize - the number of each data type
1089: . type - the data type
1090: - data - pointer to raw array
1092: Level: beginner
1094: Notes:
1095: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1097: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1098: @*/
1099: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1100: {
1101: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1102: DMSwarmDataField gfield;
1104: PetscFunctionBegin;
1106: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1107: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1108: if (data) *data = NULL;
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: /*@
1113: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
1115: Not Collective
1117: Input Parameter:
1118: . dm - a `DMSWARM`
1120: Level: beginner
1122: Notes:
1123: The new point will have all fields initialized to zero.
1125: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
1126: @*/
1127: PetscErrorCode DMSwarmAddPoint(DM dm)
1128: {
1129: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1131: PetscFunctionBegin;
1132: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1133: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1134: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
1135: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1136: PetscFunctionReturn(PETSC_SUCCESS);
1137: }
1139: /*@
1140: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
1142: Not Collective
1144: Input Parameters:
1145: + dm - a `DMSWARM`
1146: - npoints - the number of new points to add
1148: Level: beginner
1150: Notes:
1151: The new point will have all fields initialized to zero.
1153: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
1154: @*/
1155: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1156: {
1157: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1158: PetscInt nlocal;
1160: PetscFunctionBegin;
1161: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1162: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1163: nlocal = nlocal + npoints;
1164: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
1165: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1166: PetscFunctionReturn(PETSC_SUCCESS);
1167: }
1169: /*@
1170: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
1172: Not Collective
1174: Input Parameter:
1175: . dm - a `DMSWARM`
1177: Level: beginner
1179: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
1180: @*/
1181: PetscErrorCode DMSwarmRemovePoint(DM dm)
1182: {
1183: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1185: PetscFunctionBegin;
1186: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1187: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
1188: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1189: PetscFunctionReturn(PETSC_SUCCESS);
1190: }
1192: /*@
1193: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
1195: Not Collective
1197: Input Parameters:
1198: + dm - a `DMSWARM`
1199: - idx - index of point to remove
1201: Level: beginner
1203: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1204: @*/
1205: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
1206: {
1207: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1209: PetscFunctionBegin;
1210: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1211: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
1212: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1213: PetscFunctionReturn(PETSC_SUCCESS);
1214: }
1216: /*@
1217: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
1219: Not Collective
1221: Input Parameters:
1222: + dm - a `DMSWARM`
1223: . pi - the index of the point to copy
1224: - pj - the point index where the copy should be located
1226: Level: beginner
1228: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1229: @*/
1230: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
1231: {
1232: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1234: PetscFunctionBegin;
1235: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1236: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
1237: PetscFunctionReturn(PETSC_SUCCESS);
1238: }
1240: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
1241: {
1242: PetscFunctionBegin;
1243: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: /*@
1248: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
1250: Collective
1252: Input Parameters:
1253: + dm - the `DMSWARM`
1254: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1256: Level: advanced
1258: Notes:
1259: The `DM` will be modified to accommodate received points.
1260: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
1261: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
1263: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
1264: @*/
1265: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
1266: {
1267: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1269: PetscFunctionBegin;
1270: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
1271: switch (swarm->migrate_type) {
1272: case DMSWARM_MIGRATE_BASIC:
1273: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
1274: break;
1275: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1276: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
1277: break;
1278: case DMSWARM_MIGRATE_DMCELLEXACT:
1279: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1280: case DMSWARM_MIGRATE_USER:
1281: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
1282: default:
1283: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
1284: }
1285: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
1286: PetscCall(DMClearGlobalVectors(dm));
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
1292: /*
1293: DMSwarmCollectViewCreate
1295: * Applies a collection method and gathers point neighbour points into dm
1297: Notes:
1298: Users should call DMSwarmCollectViewDestroy() after
1299: they have finished computations associated with the collected points
1300: */
1302: /*@
1303: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1304: in neighbour ranks into the `DMSWARM`
1306: Collective
1308: Input Parameter:
1309: . dm - the `DMSWARM`
1311: Level: advanced
1313: Notes:
1314: Users should call `DMSwarmCollectViewDestroy()` after
1315: they have finished computations associated with the collected points
1316: Different collect methods are supported. See `DMSwarmSetCollectType()`.
1318: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
1319: @*/
1320: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1321: {
1322: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1323: PetscInt ng;
1325: PetscFunctionBegin;
1326: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
1327: PetscCall(DMSwarmGetLocalSize(dm, &ng));
1328: switch (swarm->collect_type) {
1329: case DMSWARM_COLLECT_BASIC:
1330: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
1331: break;
1332: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1333: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1334: case DMSWARM_COLLECT_GENERAL:
1335: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
1336: default:
1337: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
1338: }
1339: swarm->collect_view_active = PETSC_TRUE;
1340: swarm->collect_view_reset_nlocal = ng;
1341: PetscFunctionReturn(PETSC_SUCCESS);
1342: }
1344: /*@
1345: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
1347: Collective
1349: Input Parameters:
1350: . dm - the `DMSWARM`
1352: Notes:
1353: Users should call `DMSwarmCollectViewCreate()` before this function is called.
1355: Level: advanced
1357: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
1358: @*/
1359: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1360: {
1361: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1363: PetscFunctionBegin;
1364: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
1365: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
1366: swarm->collect_view_active = PETSC_FALSE;
1367: PetscFunctionReturn(PETSC_SUCCESS);
1368: }
1370: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
1371: {
1372: PetscInt dim;
1374: PetscFunctionBegin;
1375: PetscCall(DMSwarmSetNumSpecies(dm, 1));
1376: PetscCall(DMGetDimension(dm, &dim));
1377: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1378: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1379: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_coor, dim, PETSC_DOUBLE));
1380: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_cellid, 1, PETSC_INT));
1381: PetscFunctionReturn(PETSC_SUCCESS);
1382: }
1384: /*@
1385: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
1387: Collective
1389: Input Parameters:
1390: + dm - the `DMSWARM`
1391: - Npc - The number of particles per cell in the cell `DM`
1393: Level: intermediate
1395: Notes:
1396: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
1397: one particle is in each cell, it is placed at the centroid.
1399: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1400: @*/
1401: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1402: {
1403: DM cdm;
1404: PetscRandom rnd;
1405: DMPolytopeType ct;
1406: PetscBool simplex;
1407: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1408: PetscInt dim, d, cStart, cEnd, c, p;
1410: PetscFunctionBeginUser;
1411: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
1412: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
1413: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
1415: PetscCall(DMSwarmGetCellDM(dm, &cdm));
1416: PetscCall(DMGetDimension(cdm, &dim));
1417: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
1418: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
1419: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1421: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
1422: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1423: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1424: for (c = cStart; c < cEnd; ++c) {
1425: if (Npc == 1) {
1426: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
1427: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
1428: } else {
1429: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
1430: for (p = 0; p < Npc; ++p) {
1431: const PetscInt n = c * Npc + p;
1432: PetscReal sum = 0.0, refcoords[3];
1434: for (d = 0; d < dim; ++d) {
1435: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
1436: sum += refcoords[d];
1437: }
1438: if (simplex && sum > 0.0)
1439: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
1440: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
1441: }
1442: }
1443: }
1444: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
1445: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
1446: PetscCall(PetscRandomDestroy(&rnd));
1447: PetscFunctionReturn(PETSC_SUCCESS);
1448: }
1450: /*@
1451: DMSwarmSetType - Set particular flavor of `DMSWARM`
1453: Collective
1455: Input Parameters:
1456: + dm - the `DMSWARM`
1457: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
1459: Level: advanced
1461: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
1462: @*/
1463: PetscErrorCode DMSwarmSetType(DM dm, DMSwarmType stype)
1464: {
1465: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1467: PetscFunctionBegin;
1468: swarm->swarm_type = stype;
1469: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(dm));
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: static PetscErrorCode DMSetup_Swarm(DM dm)
1474: {
1475: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1476: PetscMPIInt rank;
1477: PetscInt p, npoints, *rankval;
1479: PetscFunctionBegin;
1480: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
1481: swarm->issetup = PETSC_TRUE;
1483: if (swarm->swarm_type == DMSWARM_PIC) {
1484: /* check dmcell exists */
1485: PetscCheck(swarm->dmcell, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires you call DMSwarmSetCellDM");
1487: if (swarm->dmcell->ops->locatepointssubdomain) {
1488: /* check methods exists for exact ownership identificiation */
1489: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
1490: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1491: } else {
1492: /* check methods exist for point location AND rank neighbor identification */
1493: if (swarm->dmcell->ops->locatepoints) {
1494: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
1495: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1497: if (swarm->dmcell->ops->getneighbors) {
1498: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
1499: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1501: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1502: }
1503: }
1505: PetscCall(DMSwarmFinalizeFieldRegister(dm));
1507: /* check some fields were registered */
1508: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
1510: /* check local sizes were set */
1511: PetscCheck(swarm->db->L != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Local sizes must be set via DMSwarmSetLocalSizes()");
1513: /* initialize values in pid and rank placeholders */
1514: /* TODO: [pid - use MPI_Scan] */
1515: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1516: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
1517: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1518: for (p = 0; p < npoints; p++) rankval[p] = (PetscInt)rank;
1519: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1520: PetscFunctionReturn(PETSC_SUCCESS);
1521: }
1523: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1525: static PetscErrorCode DMDestroy_Swarm(DM dm)
1526: {
1527: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1529: PetscFunctionBegin;
1530: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
1531: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
1532: if (swarm->sort_context) PetscCall(DMSwarmSortDestroy(&swarm->sort_context));
1533: PetscCall(PetscFree(swarm));
1534: PetscFunctionReturn(PETSC_SUCCESS);
1535: }
1537: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1538: {
1539: DM cdm;
1540: PetscDraw draw;
1541: PetscReal *coords, oldPause, radius = 0.01;
1542: PetscInt Np, p, bs;
1544: PetscFunctionBegin;
1545: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
1546: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1547: PetscCall(DMSwarmGetCellDM(dm, &cdm));
1548: PetscCall(PetscDrawGetPause(draw, &oldPause));
1549: PetscCall(PetscDrawSetPause(draw, 0.0));
1550: PetscCall(DMView(cdm, viewer));
1551: PetscCall(PetscDrawSetPause(draw, oldPause));
1553: PetscCall(DMSwarmGetLocalSize(dm, &Np));
1554: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords));
1555: for (p = 0; p < Np; ++p) {
1556: const PetscInt i = p * bs;
1558: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
1559: }
1560: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords));
1561: PetscCall(PetscDrawFlush(draw));
1562: PetscCall(PetscDrawPause(draw));
1563: PetscCall(PetscDrawSave(draw));
1564: PetscFunctionReturn(PETSC_SUCCESS);
1565: }
1567: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
1568: {
1569: PetscViewerFormat format;
1570: PetscInt *sizes;
1571: PetscInt dim, Np, maxSize = 17;
1572: MPI_Comm comm;
1573: PetscMPIInt rank, size, p;
1574: const char *name;
1576: PetscFunctionBegin;
1577: PetscCall(PetscViewerGetFormat(viewer, &format));
1578: PetscCall(DMGetDimension(dm, &dim));
1579: PetscCall(DMSwarmGetLocalSize(dm, &Np));
1580: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1581: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1582: PetscCallMPI(MPI_Comm_size(comm, &size));
1583: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1584: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
1585: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
1586: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
1587: else PetscCall(PetscCalloc1(3, &sizes));
1588: if (size < maxSize) {
1589: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
1590: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
1591: for (p = 0; p < size; ++p) {
1592: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
1593: }
1594: } else {
1595: PetscInt locMinMax[2] = {Np, Np};
1597: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
1598: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
1599: }
1600: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1601: PetscCall(PetscFree(sizes));
1602: if (format == PETSC_VIEWER_ASCII_INFO) {
1603: PetscInt *cell;
1605: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
1606: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1607: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell));
1608: for (p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%d: %" PetscInt_FMT "\n", p, cell[p]));
1609: PetscCall(PetscViewerFlush(viewer));
1610: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1611: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell));
1612: }
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1617: {
1618: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1619: PetscBool iascii, ibinary, isvtk, isdraw;
1620: #if defined(PETSC_HAVE_HDF5)
1621: PetscBool ishdf5;
1622: #endif
1624: PetscFunctionBegin;
1627: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1628: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
1629: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1630: #if defined(PETSC_HAVE_HDF5)
1631: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1632: #endif
1633: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1634: if (iascii) {
1635: PetscViewerFormat format;
1637: PetscCall(PetscViewerGetFormat(viewer, &format));
1638: switch (format) {
1639: case PETSC_VIEWER_ASCII_INFO_DETAIL:
1640: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
1641: break;
1642: default:
1643: PetscCall(DMView_Swarm_Ascii(dm, viewer));
1644: }
1645: } else {
1646: #if defined(PETSC_HAVE_HDF5)
1647: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
1648: #endif
1649: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
1650: }
1651: PetscFunctionReturn(PETSC_SUCCESS);
1652: }
1654: /*@C
1655: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
1656: 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.
1658: Noncollective
1660: Input Parameters:
1661: + sw - the `DMSWARM`
1662: . cellID - the integer id of the cell to be extracted and filtered
1663: - cellswarm - The `DMSWARM` to receive the cell
1665: Level: beginner
1667: Notes:
1668: This presently only supports `DMSWARM_PIC` type
1670: Should be restored with `DMSwarmRestoreCellSwarm()`
1672: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
1674: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
1675: @*/
1676: PETSC_EXTERN PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1677: {
1678: DM_Swarm *original = (DM_Swarm *)sw->data;
1679: DMLabel label;
1680: DM dmc, subdmc;
1681: PetscInt *pids, particles, dim;
1683: PetscFunctionBegin;
1684: /* Configure new swarm */
1685: PetscCall(DMSetType(cellswarm, DMSWARM));
1686: PetscCall(DMGetDimension(sw, &dim));
1687: PetscCall(DMSetDimension(cellswarm, dim));
1688: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
1689: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
1690: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
1691: PetscCall(DMSwarmSortGetAccess(sw));
1692: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
1693: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
1694: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
1695: PetscCall(DMSwarmSortRestoreAccess(sw));
1696: PetscCall(PetscFree(pids));
1697: PetscCall(DMSwarmGetCellDM(sw, &dmc));
1698: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
1699: PetscCall(DMAddLabel(dmc, label));
1700: PetscCall(DMLabelSetValue(label, cellID, 1));
1701: PetscCall(DMPlexFilter(dmc, label, 1, &subdmc));
1702: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
1703: PetscCall(DMLabelDestroy(&label));
1704: PetscFunctionReturn(PETSC_SUCCESS);
1705: }
1707: /*@C
1708: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
1710: Noncollective
1712: Input Parameters:
1713: + sw - the parent `DMSWARM`
1714: . cellID - the integer id of the cell to be copied back into the parent swarm
1715: - cellswarm - the cell swarm object
1717: Level: beginner
1719: Note:
1720: This only supports `DMSWARM_PIC` types of `DMSWARM`s
1722: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
1723: @*/
1724: PETSC_EXTERN PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1725: {
1726: DM dmc;
1727: PetscInt *pids, particles, p;
1729: PetscFunctionBegin;
1730: PetscCall(DMSwarmSortGetAccess(sw));
1731: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
1732: PetscCall(DMSwarmSortRestoreAccess(sw));
1733: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
1734: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
1735: /* Free memory, destroy cell dm */
1736: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
1737: PetscCall(DMDestroy(&dmc));
1738: PetscCall(PetscFree(pids));
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: }
1742: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
1744: static PetscErrorCode DMInitialize_Swarm(DM sw)
1745: {
1746: PetscFunctionBegin;
1747: sw->dim = 0;
1748: sw->ops->view = DMView_Swarm;
1749: sw->ops->load = NULL;
1750: sw->ops->setfromoptions = NULL;
1751: sw->ops->clone = DMClone_Swarm;
1752: sw->ops->setup = DMSetup_Swarm;
1753: sw->ops->createlocalsection = NULL;
1754: sw->ops->createdefaultconstraints = NULL;
1755: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1756: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
1757: sw->ops->getlocaltoglobalmapping = NULL;
1758: sw->ops->createfieldis = NULL;
1759: sw->ops->createcoordinatedm = NULL;
1760: sw->ops->getcoloring = NULL;
1761: sw->ops->creatematrix = DMCreateMatrix_Swarm;
1762: sw->ops->createinterpolation = NULL;
1763: sw->ops->createinjection = NULL;
1764: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1765: sw->ops->refine = NULL;
1766: sw->ops->coarsen = NULL;
1767: sw->ops->refinehierarchy = NULL;
1768: sw->ops->coarsenhierarchy = NULL;
1769: sw->ops->globaltolocalbegin = NULL;
1770: sw->ops->globaltolocalend = NULL;
1771: sw->ops->localtoglobalbegin = NULL;
1772: sw->ops->localtoglobalend = NULL;
1773: sw->ops->destroy = DMDestroy_Swarm;
1774: sw->ops->createsubdm = NULL;
1775: sw->ops->getdimpoints = NULL;
1776: sw->ops->locatepoints = NULL;
1777: PetscFunctionReturn(PETSC_SUCCESS);
1778: }
1780: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
1781: {
1782: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1784: PetscFunctionBegin;
1785: swarm->refct++;
1786: (*newdm)->data = swarm;
1787: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
1788: PetscCall(DMInitialize_Swarm(*newdm));
1789: (*newdm)->dim = dm->dim;
1790: PetscFunctionReturn(PETSC_SUCCESS);
1791: }
1793: /*MC
1795: DMSWARM = "swarm" - A `DM` object used to represent arrays of data (fields) of arbitrary data type.
1796: This implementation was designed for particle methods in which the underlying
1797: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1799: Level: intermediate
1801: Notes:
1802: User data can be represented by `DMSWARM` through a registering "fields".
1803: To register a field, the user must provide;
1804: (a) a unique name;
1805: (b) the data type (or size in bytes);
1806: (c) the block size of the data.
1808: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1809: on a set of particles. Then the following code could be used
1810: .vb
1811: DMSwarmInitializeFieldRegister(dm)
1812: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1813: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1814: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1815: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1816: DMSwarmFinalizeFieldRegister(dm)
1817: .ve
1819: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
1820: The only restriction imposed by `DMSWARM` is that all fields contain the same number of points.
1822: To support particle methods, "migration" techniques are provided. These methods migrate data
1823: between ranks.
1825: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
1826: As a `DMSWARM` may internally define and store values of different data types,
1827: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
1828: fields should be used to define a `Vec` object via
1829: `DMSwarmVectorDefineField()`
1830: The specified field can be changed at any time - thereby permitting vectors
1831: compatible with different fields to be created.
1833: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via
1834: `DMSwarmCreateGlobalVectorFromField()`
1835: Here the data defining the field in the `DMSWARM` is shared with a Vec.
1836: This is inherently unsafe if you alter the size of the field at any time between
1837: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
1838: If the local size of the `DMSWARM` does not match the local size of the global vector
1839: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
1841: Additional high-level support is provided for Particle-In-Cell methods.
1842: Please refer to `DMSwarmSetType()`.
1844: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`
1845: M*/
1846: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1847: {
1848: DM_Swarm *swarm;
1850: PetscFunctionBegin;
1852: PetscCall(PetscNew(&swarm));
1853: dm->data = swarm;
1854: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
1855: PetscCall(DMSwarmInitializeFieldRegister(dm));
1856: swarm->refct = 1;
1857: swarm->vec_field_set = PETSC_FALSE;
1858: swarm->issetup = PETSC_FALSE;
1859: swarm->swarm_type = DMSWARM_BASIC;
1860: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1861: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1862: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1863: swarm->dmcell = NULL;
1864: swarm->collect_view_active = PETSC_FALSE;
1865: swarm->collect_view_reset_nlocal = -1;
1866: PetscCall(DMInitialize_Swarm(dm));
1867: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
1868: PetscFunctionReturn(PETSC_SUCCESS);
1869: }