Actual source code: swarm.c
1: #include "petscdmswarm.h"
2: #define PETSCDM_DLL
3: #include <petsc/private/dmswarmimpl.h>
4: #include <petsc/private/hashsetij.h>
5: #include <petsc/private/petscfeimpl.h>
6: #include <petscviewer.h>
7: #include <petscdraw.h>
8: #include <petscdmplex.h>
9: #include <petscblaslapack.h>
10: #include "../src/dm/impls/swarm/data_bucket.h"
11: #include <petscdmlabel.h>
12: #include <petscsection.h>
14: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
15: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
16: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
18: const char *DMSwarmTypeNames[] = {"basic", "pic", NULL};
19: const char *DMSwarmMigrateTypeNames[] = {"basic", "dmcellnscatter", "dmcellexact", "user", NULL};
20: const char *DMSwarmCollectTypeNames[] = {"basic", "boundingbox", "general", "user", NULL};
21: const char *DMSwarmRemapTypeNames[] = {"none", "pfak", "colella", "DMSwarmRemapType", "DMSWARM_REMAP_", NULL};
22: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
24: const char DMSwarmField_pid[] = "DMSwarm_pid";
25: const char DMSwarmField_rank[] = "DMSwarm_rank";
26: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
28: PetscInt SwarmDataFieldId = -1;
30: #if defined(PETSC_HAVE_HDF5)
31: #include <petscviewerhdf5.h>
33: static PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
34: {
35: DM dm;
36: PetscReal seqval;
37: PetscInt seqnum, bs;
38: PetscBool isseq, ists;
40: PetscFunctionBegin;
41: PetscCall(VecGetDM(v, &dm));
42: PetscCall(VecGetBlockSize(v, &bs));
43: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particle_fields"));
44: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
45: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
46: PetscCall(PetscViewerHDF5IsTimestepping(viewer, &ists));
47: if (ists) PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
48: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
49: PetscCall(VecViewNative(v, viewer));
50: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
51: PetscCall(PetscViewerHDF5PopGroup(viewer));
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
56: {
57: DMSwarmCellDM celldm;
58: Vec coordinates;
59: PetscInt Np, Nfc;
60: PetscBool isseq;
61: const char **coordFields;
63: PetscFunctionBegin;
64: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
65: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
66: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
67: PetscCall(DMSwarmGetSize(dm, &Np));
68: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, coordFields[0], &coordinates));
69: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
70: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
71: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
72: PetscCall(VecViewNative(coordinates, viewer));
73: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
74: PetscCall(PetscViewerHDF5PopGroup(viewer));
75: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, coordFields[0], &coordinates));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
78: #endif
80: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
81: {
82: DM dm;
83: #if defined(PETSC_HAVE_HDF5)
84: PetscBool ishdf5;
85: #endif
87: PetscFunctionBegin;
88: PetscCall(VecGetDM(v, &dm));
89: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
90: #if defined(PETSC_HAVE_HDF5)
91: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
92: if (ishdf5) {
93: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
96: #endif
97: PetscCall(VecViewNative(v, viewer));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /*@C
102: DMSwarmVectorGetField - Gets the fields from which to define a `Vec` object
103: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
105: Not collective
107: Input Parameter:
108: . sw - a `DMSWARM`
110: Output Parameters:
111: + Nf - the number of fields
112: - fieldnames - the textual name given to each registered field, or NULL if it has not been set
114: Level: beginner
116: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
117: @*/
118: PetscErrorCode DMSwarmVectorGetField(DM sw, PetscInt *Nf, const char **fieldnames[])
119: {
120: DMSwarmCellDM celldm;
122: PetscFunctionBegin;
124: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
125: PetscCall(DMSwarmCellDMGetFields(celldm, Nf, fieldnames));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*@
130: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
131: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
133: Collective
135: Input Parameters:
136: + dm - a `DMSWARM`
137: - fieldname - the textual name given to each registered field
139: Level: beginner
141: Notes:
142: The field with name `fieldname` must be defined as having a data type of `PetscScalar`.
144: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
145: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
147: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineFields()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
148: @*/
149: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
150: {
151: PetscFunctionBegin;
152: PetscCall(DMSwarmVectorDefineFields(dm, 1, &fieldname));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: /*@C
157: DMSwarmVectorDefineFields - Sets the fields from which to define a `Vec` object
158: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
160: Collective, No Fortran support
162: Input Parameters:
163: + sw - a `DMSWARM`
164: . Nf - the number of fields
165: - fieldnames - the textual name given to each registered field
167: Level: beginner
169: Notes:
170: Each field with name in `fieldnames` must be defined as having a data type of `PetscScalar`.
172: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
173: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
175: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
176: @*/
177: PetscErrorCode DMSwarmVectorDefineFields(DM sw, PetscInt Nf, const char *fieldnames[])
178: {
179: DM_Swarm *swarm = (DM_Swarm *)sw->data;
180: DMSwarmCellDM celldm;
182: PetscFunctionBegin;
184: if (fieldnames) PetscAssertPointer(fieldnames, 3);
185: if (!swarm->issetup) PetscCall(DMSetUp(sw));
186: PetscCheck(Nf >= 0, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields must be non-negative, not %" PetscInt_FMT, Nf);
187: // Create a dummy cell DM if none has been specified (I think we should not support this mode)
188: if (!swarm->activeCellDM) {
189: DM dm;
190: DMSwarmCellDM celldm;
192: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), &dm));
193: PetscCall(DMSetType(dm, DMSHELL));
194: PetscCall(PetscObjectSetName((PetscObject)dm, "dummy"));
195: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 0, NULL, &celldm));
196: PetscCall(DMDestroy(&dm));
197: PetscCall(DMSwarmAddCellDM(sw, celldm));
198: PetscCall(DMSwarmCellDMDestroy(&celldm));
199: PetscCall(DMSwarmSetCellDMActive(sw, "dummy"));
200: }
201: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
202: for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscFree(celldm->dmFields[f]));
203: PetscCall(PetscFree(celldm->dmFields));
205: celldm->Nf = Nf;
206: PetscCall(PetscMalloc1(Nf, &celldm->dmFields));
207: for (PetscInt f = 0; f < Nf; ++f) {
208: PetscDataType type;
210: // Check all fields are of type PETSC_REAL or PETSC_SCALAR
211: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], NULL, &type));
212: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
213: PetscCall(PetscStrallocpy(fieldnames[f], (char **)&celldm->dmFields[f]));
214: }
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /* requires DMSwarmDefineFieldVector has been called */
219: static PetscErrorCode DMCreateGlobalVector_Swarm(DM sw, Vec *vec)
220: {
221: DM_Swarm *swarm = (DM_Swarm *)sw->data;
222: DMSwarmCellDM celldm;
223: Vec x;
224: char name[PETSC_MAX_PATH_LEN];
225: PetscInt bs = 0, n;
227: PetscFunctionBegin;
228: if (!swarm->issetup) PetscCall(DMSetUp(sw));
229: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
230: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
231: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
233: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
234: for (PetscInt f = 0; f < celldm->Nf; ++f) {
235: PetscInt fbs;
236: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
237: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
238: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
239: bs += fbs;
240: }
241: PetscCall(VecCreate(PetscObjectComm((PetscObject)sw), &x));
242: PetscCall(PetscObjectSetName((PetscObject)x, name));
243: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
244: PetscCall(VecSetBlockSize(x, bs));
245: PetscCall(VecSetDM(x, sw));
246: PetscCall(VecSetFromOptions(x));
247: PetscCall(VecSetOperation(x, VECOP_VIEW, (void (*)(void))VecView_Swarm));
248: *vec = x;
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /* requires DMSwarmDefineFieldVector has been called */
253: static PetscErrorCode DMCreateLocalVector_Swarm(DM sw, Vec *vec)
254: {
255: DM_Swarm *swarm = (DM_Swarm *)sw->data;
256: DMSwarmCellDM celldm;
257: Vec x;
258: char name[PETSC_MAX_PATH_LEN];
259: PetscInt bs = 0, n;
261: PetscFunctionBegin;
262: if (!swarm->issetup) PetscCall(DMSetUp(sw));
263: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
264: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
265: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
267: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
268: for (PetscInt f = 0; f < celldm->Nf; ++f) {
269: PetscInt fbs;
270: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
271: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
272: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
273: bs += fbs;
274: }
275: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
276: PetscCall(PetscObjectSetName((PetscObject)x, name));
277: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
278: PetscCall(VecSetBlockSize(x, bs));
279: PetscCall(VecSetDM(x, sw));
280: PetscCall(VecSetFromOptions(x));
281: *vec = x;
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
286: {
287: DM_Swarm *swarm = (DM_Swarm *)dm->data;
288: DMSwarmDataField gfield;
289: PetscInt bs, nlocal, fid = -1, cfid = -2;
290: PetscBool flg;
292: PetscFunctionBegin;
293: /* check vector is an inplace array */
294: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
295: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
296: (void)flg; /* avoid compiler warning */
297: PetscCheck(cfid == fid, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)! %" PetscInt_FMT " != %" PetscInt_FMT, fieldname, cfid, fid);
298: PetscCall(VecGetLocalSize(*vec, &nlocal));
299: PetscCall(VecGetBlockSize(*vec, &bs));
300: PetscCheck(nlocal / bs == swarm->db->L, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
301: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
302: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
303: PetscCall(VecResetArray(*vec));
304: PetscCall(VecDestroy(vec));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
309: {
310: DM_Swarm *swarm = (DM_Swarm *)dm->data;
311: PetscDataType type;
312: PetscScalar *array;
313: PetscInt bs, n, fid;
314: char name[PETSC_MAX_PATH_LEN];
315: PetscMPIInt size;
316: PetscBool iscuda, iskokkos, iship;
318: PetscFunctionBegin;
319: if (!swarm->issetup) PetscCall(DMSetUp(dm));
320: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
321: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
322: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
324: PetscCallMPI(MPI_Comm_size(comm, &size));
325: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
326: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
327: PetscCall(PetscStrcmp(dm->vectype, VECHIP, &iship));
328: PetscCall(VecCreate(comm, vec));
329: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
330: PetscCall(VecSetBlockSize(*vec, bs));
331: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
332: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
333: else if (iship) PetscCall(VecSetType(*vec, VECHIP));
334: else PetscCall(VecSetType(*vec, VECSTANDARD));
335: PetscCall(VecPlaceArray(*vec, array));
337: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
338: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
340: /* Set guard */
341: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
342: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
344: PetscCall(VecSetDM(*vec, dm));
345: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: static PetscErrorCode DMSwarmDestroyVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], Vec *vec)
350: {
351: DM_Swarm *swarm = (DM_Swarm *)sw->data;
352: const PetscScalar *array;
353: PetscInt bs, n, id = 0, cid = -2;
354: PetscBool flg;
356: PetscFunctionBegin;
357: for (PetscInt f = 0; f < Nf; ++f) {
358: PetscInt fid;
360: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
361: id += fid;
362: }
363: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cid, flg));
364: (void)flg; /* avoid compiler warning */
365: PetscCheck(cid == id, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)! %" PetscInt_FMT " != %" PetscInt_FMT, fieldnames[0], cid, id);
366: PetscCall(VecGetLocalSize(*vec, &n));
367: PetscCall(VecGetBlockSize(*vec, &bs));
368: n /= bs;
369: PetscCheck(n == swarm->db->L, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
370: PetscCall(VecGetArrayRead(*vec, &array));
371: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
372: PetscScalar *farray;
373: PetscDataType ftype;
374: PetscInt fbs;
376: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
377: PetscCheck(off + fbs <= bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid blocksize %" PetscInt_FMT " < %" PetscInt_FMT, bs, off + fbs);
378: for (PetscInt i = 0; i < n; ++i) {
379: for (PetscInt b = 0; b < fbs; ++b) farray[i * fbs + b] = array[i * bs + off + b];
380: }
381: off += fbs;
382: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
383: }
384: PetscCall(VecRestoreArrayRead(*vec, &array));
385: PetscCall(VecDestroy(vec));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: static PetscErrorCode DMSwarmCreateVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], MPI_Comm comm, Vec *vec)
390: {
391: DM_Swarm *swarm = (DM_Swarm *)sw->data;
392: PetscScalar *array;
393: PetscInt n, bs = 0, id = 0;
394: char name[PETSC_MAX_PATH_LEN];
396: PetscFunctionBegin;
397: if (!swarm->issetup) PetscCall(DMSetUp(sw));
398: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
399: for (PetscInt f = 0; f < Nf; ++f) {
400: PetscDataType ftype;
401: PetscInt fbs;
403: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
404: PetscCheck(ftype == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
405: bs += fbs;
406: }
408: PetscCall(VecCreate(comm, vec));
409: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
410: PetscCall(VecSetBlockSize(*vec, bs));
411: PetscCall(VecSetType(*vec, sw->vectype));
413: PetscCall(VecGetArrayWrite(*vec, &array));
414: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
415: PetscScalar *farray;
416: PetscDataType ftype;
417: PetscInt fbs;
419: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
420: for (PetscInt i = 0; i < n; ++i) {
421: for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = farray[i * fbs + b];
422: }
423: off += fbs;
424: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
425: }
426: PetscCall(VecRestoreArrayWrite(*vec, &array));
428: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
429: for (PetscInt f = 0; f < Nf; ++f) {
430: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
431: PetscCall(PetscStrlcat(name, fieldnames[f], PETSC_MAX_PATH_LEN));
432: }
433: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
435: for (PetscInt f = 0; f < Nf; ++f) {
436: PetscInt fid;
438: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
439: id += fid;
440: }
441: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, id));
443: PetscCall(VecSetDM(*vec, sw));
444: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
450: \hat f = \sum_i f_i \phi_i
452: and a particle function is given by
454: f = \sum_i w_i \delta(x - x_i)
456: then we want to require that
458: M \hat f = M_p f
460: where the particle mass matrix is given by
462: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
464: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
465: his integral. We allow this with the boolean flag.
466: */
467: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
468: {
469: const char *name = "Mass Matrix";
470: MPI_Comm comm;
471: DMSwarmCellDM celldm;
472: PetscDS prob;
473: PetscSection fsection, globalFSection;
474: PetscHSetIJ ht;
475: PetscLayout rLayout, colLayout;
476: PetscInt *dnz, *onz;
477: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
478: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
479: PetscScalar *elemMat;
480: PetscInt dim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0, totNc = 0;
481: const char **coordFields;
482: PetscReal **coordVals;
483: PetscInt *bs;
485: PetscFunctionBegin;
486: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
487: PetscCall(DMGetCoordinateDim(dmf, &dim));
488: PetscCall(DMGetDS(dmf, &prob));
489: PetscCall(PetscDSGetNumFields(prob, &Nf));
490: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
491: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
492: PetscCall(DMGetLocalSection(dmf, &fsection));
493: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
494: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
495: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
497: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
498: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
499: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
501: PetscCall(PetscLayoutCreate(comm, &colLayout));
502: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
503: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
504: PetscCall(PetscLayoutSetUp(colLayout));
505: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
506: PetscCall(PetscLayoutDestroy(&colLayout));
508: PetscCall(PetscLayoutCreate(comm, &rLayout));
509: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
510: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
511: PetscCall(PetscLayoutSetUp(rLayout));
512: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
513: PetscCall(PetscLayoutDestroy(&rLayout));
515: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
516: PetscCall(PetscHSetIJCreate(&ht));
518: PetscCall(PetscSynchronizedFlush(comm, NULL));
519: for (PetscInt field = 0; field < Nf; ++field) {
520: PetscObject obj;
521: PetscClassId id;
522: PetscInt Nc;
524: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
525: PetscCall(PetscObjectGetClassId(obj, &id));
526: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
527: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
528: totNc += Nc;
529: }
530: /* count non-zeros */
531: PetscCall(DMSwarmSortGetAccess(dmc));
532: for (PetscInt field = 0; field < Nf; ++field) {
533: PetscObject obj;
534: PetscClassId id;
535: PetscInt Nc;
537: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
538: PetscCall(PetscObjectGetClassId(obj, &id));
539: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
540: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
542: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
543: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
544: PetscInt numFIndices, numCIndices;
546: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
547: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
548: maxC = PetscMax(maxC, numCIndices);
549: {
550: PetscHashIJKey key;
551: PetscBool missing;
552: for (PetscInt i = 0; i < numFIndices; ++i) {
553: key.j = findices[i]; /* global column (from Plex) */
554: if (key.j >= 0) {
555: /* Get indices for coarse elements */
556: for (PetscInt j = 0; j < numCIndices; ++j) {
557: for (PetscInt c = 0; c < Nc; ++c) {
558: // TODO Need field offset on particle here
559: key.i = cindices[j] * totNc + c + rStart; /* global cols (from Swarm) */
560: if (key.i < 0) continue;
561: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
562: PetscCheck(missing, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
563: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
564: else ++onz[key.i - rStart];
565: }
566: }
567: }
568: }
569: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
570: }
571: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
572: }
573: }
574: PetscCall(PetscHSetIJDestroy(&ht));
575: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
576: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
577: PetscCall(PetscFree2(dnz, onz));
578: PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
579: for (PetscInt field = 0; field < Nf; ++field) {
580: PetscTabulation Tcoarse;
581: PetscObject obj;
582: PetscClassId id;
583: PetscInt Nc;
585: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
586: PetscCall(PetscObjectGetClassId(obj, &id));
587: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
588: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
590: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
591: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
592: PetscInt *findices, *cindices;
593: PetscInt numFIndices, numCIndices;
595: /* TODO: Use DMField instead of assuming affine */
596: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
597: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
598: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
599: for (PetscInt j = 0; j < numCIndices; ++j) {
600: PetscReal xr[8];
601: PetscInt off = 0;
603: for (PetscInt i = 0; i < Nfc; ++i) {
604: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[j] * bs[i] + b];
605: }
606: PetscCheck(off == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, dim);
607: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, xr, &xi[j * dim]);
608: }
609: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
610: else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
611: /* Get elemMat entries by multiplying by weight */
612: PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
613: for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
614: for (PetscInt j = 0; j < numCIndices; ++j) {
615: for (PetscInt c = 0; c < Nc; ++c) {
616: // TODO Need field offset on particle and field here
617: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
618: elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
619: }
620: }
621: }
622: for (PetscInt j = 0; j < numCIndices; ++j)
623: // TODO Need field offset on particle here
624: for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
625: if (0) PetscCall(DMPrintCellMatrix(cell, name, numCIndices * Nc, numFIndices, elemMat));
626: PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
627: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
628: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
629: PetscCall(PetscTabulationDestroy(&Tcoarse));
630: }
631: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
632: }
633: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
634: PetscCall(DMSwarmSortRestoreAccess(dmc));
635: PetscCall(PetscFree3(v0, J, invJ));
636: PetscCall(PetscFree2(coordVals, bs));
637: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
638: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: /* Returns empty matrix for use with SNES FD */
643: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
644: {
645: Vec field;
646: PetscInt size;
648: PetscFunctionBegin;
649: PetscCall(DMGetGlobalVector(sw, &field));
650: PetscCall(VecGetLocalSize(field, &size));
651: PetscCall(DMRestoreGlobalVector(sw, &field));
652: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
653: PetscCall(MatSetFromOptions(*m));
654: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
655: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
656: PetscCall(MatZeroEntries(*m));
657: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
658: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
659: PetscCall(MatShift(*m, 1.0));
660: PetscCall(MatSetDM(*m, sw));
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: /* FEM cols, Particle rows */
665: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
666: {
667: DMSwarmCellDM celldm;
668: PetscSection gsf;
669: PetscInt m, n, Np, bs;
670: void *ctx;
672: PetscFunctionBegin;
673: PetscCall(DMSwarmGetCellDMActive(dmCoarse, &celldm));
674: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)dmCoarse), PETSC_ERR_USER, "Active cell DM does not define any fields");
675: PetscCall(DMGetGlobalSection(dmFine, &gsf));
676: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
677: PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
678: PetscCall(DMSwarmCellDMGetBlockSize(celldm, dmCoarse, &bs));
679: n = Np * bs;
680: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
681: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
682: PetscCall(MatSetType(*mass, dmCoarse->mattype));
683: PetscCall(DMGetApplicationContext(dmFine, &ctx));
685: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
686: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
691: {
692: const char *name = "Mass Matrix Square";
693: MPI_Comm comm;
694: DMSwarmCellDM celldm;
695: PetscDS prob;
696: PetscSection fsection, globalFSection;
697: PetscHSetIJ ht;
698: PetscLayout rLayout, colLayout;
699: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
700: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
701: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
702: PetscScalar *elemMat, *elemMatSq;
703: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0;
704: const char **coordFields;
705: PetscReal **coordVals;
706: PetscInt *bs;
708: PetscFunctionBegin;
709: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
710: PetscCall(DMGetCoordinateDim(dmf, &cdim));
711: PetscCall(DMGetDS(dmf, &prob));
712: PetscCall(PetscDSGetNumFields(prob, &Nf));
713: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
714: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
715: PetscCall(DMGetLocalSection(dmf, &fsection));
716: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
717: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
718: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
720: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
721: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
722: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
724: PetscCall(PetscLayoutCreate(comm, &colLayout));
725: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
726: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
727: PetscCall(PetscLayoutSetUp(colLayout));
728: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
729: PetscCall(PetscLayoutDestroy(&colLayout));
731: PetscCall(PetscLayoutCreate(comm, &rLayout));
732: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
733: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
734: PetscCall(PetscLayoutSetUp(rLayout));
735: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
736: PetscCall(PetscLayoutDestroy(&rLayout));
738: PetscCall(DMPlexGetDepth(dmf, &depth));
739: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
740: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
741: PetscCall(PetscMalloc1(maxAdjSize, &adj));
743: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
744: PetscCall(PetscHSetIJCreate(&ht));
745: /* Count nonzeros
746: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
747: */
748: PetscCall(DMSwarmSortGetAccess(dmc));
749: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
750: PetscInt *cindices;
751: PetscInt numCIndices;
752: #if 0
753: PetscInt adjSize = maxAdjSize, a, j;
754: #endif
756: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
757: maxC = PetscMax(maxC, numCIndices);
758: /* Diagonal block */
759: for (PetscInt i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
760: #if 0
761: /* Off-diagonal blocks */
762: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
763: for (a = 0; a < adjSize; ++a) {
764: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
765: const PetscInt ncell = adj[a];
766: PetscInt *ncindices;
767: PetscInt numNCIndices;
769: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
770: {
771: PetscHashIJKey key;
772: PetscBool missing;
774: for (i = 0; i < numCIndices; ++i) {
775: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
776: if (key.i < 0) continue;
777: for (j = 0; j < numNCIndices; ++j) {
778: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
779: if (key.j < 0) continue;
780: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
781: if (missing) {
782: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
783: else ++onz[key.i - rStart];
784: }
785: }
786: }
787: }
788: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
789: }
790: }
791: #endif
792: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
793: }
794: PetscCall(PetscHSetIJDestroy(&ht));
795: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
796: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
797: PetscCall(PetscFree2(dnz, onz));
798: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
799: /* Fill in values
800: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
801: Start just by producing block diagonal
802: Could loop over adjacent cells
803: Produce neighboring element matrix
804: TODO Determine which columns and rows correspond to shared dual vector
805: Do MatMatMult with rectangular matrices
806: Insert block
807: */
808: for (PetscInt field = 0; field < Nf; ++field) {
809: PetscTabulation Tcoarse;
810: PetscObject obj;
811: PetscInt Nc;
813: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
814: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
815: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
816: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
817: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
818: PetscInt *findices, *cindices;
819: PetscInt numFIndices, numCIndices;
821: /* TODO: Use DMField instead of assuming affine */
822: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
823: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
824: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
825: for (PetscInt p = 0; p < numCIndices; ++p) {
826: PetscReal xr[8];
827: PetscInt off = 0;
829: for (PetscInt i = 0; i < Nfc; ++i) {
830: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[p] * bs[i] + b];
831: }
832: PetscCheck(off == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, cdim);
833: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[p * cdim]);
834: }
835: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
836: /* Get elemMat entries by multiplying by weight */
837: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
838: for (PetscInt i = 0; i < numFIndices; ++i) {
839: for (PetscInt p = 0; p < numCIndices; ++p) {
840: for (PetscInt c = 0; c < Nc; ++c) {
841: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
842: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
843: }
844: }
845: }
846: PetscCall(PetscTabulationDestroy(&Tcoarse));
847: for (PetscInt p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
848: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
849: /* Block diagonal */
850: if (numCIndices) {
851: PetscBLASInt blasn, blask;
852: PetscScalar one = 1.0, zero = 0.0;
854: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
855: PetscCall(PetscBLASIntCast(numFIndices, &blask));
856: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
857: }
858: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
859: /* TODO off-diagonal */
860: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
861: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
862: }
863: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
864: }
865: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
866: PetscCall(PetscFree(adj));
867: PetscCall(DMSwarmSortRestoreAccess(dmc));
868: PetscCall(PetscFree3(v0, J, invJ));
869: PetscCall(PetscFree2(coordVals, bs));
870: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
871: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
872: PetscFunctionReturn(PETSC_SUCCESS);
873: }
875: /*@
876: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
878: Collective
880: Input Parameters:
881: + dmCoarse - a `DMSWARM`
882: - dmFine - a `DMPLEX`
884: Output Parameter:
885: . mass - the square of the particle mass matrix
887: Level: advanced
889: Note:
890: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
891: future to compute the full normal equations.
893: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
894: @*/
895: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
896: {
897: PetscInt n;
898: void *ctx;
900: PetscFunctionBegin;
901: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
902: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
903: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
904: PetscCall(MatSetType(*mass, dmCoarse->mattype));
905: PetscCall(DMGetApplicationContext(dmFine, &ctx));
907: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
908: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: /* This creates a "gradient matrix" between a finite element and particle space, which is meant to enforce a weak divergence condition on the particle space. We are looking for a finite element field that has the same divergence as our particle field, so that
914: \int_X \psi_i \nabla \cdot \hat f = \int_X \psi_i \nabla \cdot f
916: and then integrate by parts
918: \int_X \nabla \psi_i \cdot \hat f = \int_X \nabla \psi_i \cdot f
920: where \psi is from a scalar FE space. If a finite element interpolant is given by
922: \hat f^c = \sum_i f_i \phi^c_i
924: and a particle function is given by
926: f^c = \sum_p f^c_p \delta(x - x_p)
928: then we want to require that
930: D_f \hat f = D_p f
932: where the gradient matrices are given by
934: (D_f)_{i(jc)} = \int \partial_c \psi_i \phi_j
935: (D_p)_{i(jc)} = \int \partial_c \psi_i \delta(x - x_j)
937: Thus we need two finite element spaces, a scalar and a vector. The vector space holds the representer for the
938: vector particle field. The scalar space holds the output of D_p or D_f, which is the weak divergence of the field.
940: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
941: his integral. We allow this with the boolean flag.
942: */
943: static PetscErrorCode DMSwarmComputeGradientMatrix_Private(DM sw, DM dm, Mat derv, PetscBool useDeltaFunction, void *ctx)
944: {
945: const char *name = "Derivative Matrix";
946: MPI_Comm comm;
947: DMSwarmCellDM celldm;
948: PetscDS ds;
949: PetscSection fsection, globalFSection;
950: PetscLayout rLayout;
951: PetscInt locRows, rStart, *rowIDXs;
952: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
953: PetscScalar *elemMat;
954: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxNpc = 0, totNc = 0;
955: const char **coordFields;
956: PetscReal **coordVals;
957: PetscInt *bs;
959: PetscFunctionBegin;
960: PetscCall(PetscObjectGetComm((PetscObject)derv, &comm));
961: PetscCall(DMGetCoordinateDim(dm, &cdim));
962: PetscCall(DMGetDS(dm, &ds));
963: PetscCall(PetscDSGetNumFields(ds, &Nf));
964: PetscCheck(Nf == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
965: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
966: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
967: PetscCall(DMGetLocalSection(dm, &fsection));
968: PetscCall(DMGetGlobalSection(dm, &globalFSection));
969: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
970: PetscCall(MatGetLocalSize(derv, &locRows, NULL));
972: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
973: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
974: PetscCheck(Nfc == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
975: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
977: PetscCall(PetscLayoutCreate(comm, &rLayout));
978: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
979: PetscCall(PetscLayoutSetBlockSize(rLayout, cdim));
980: PetscCall(PetscLayoutSetUp(rLayout));
981: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
982: PetscCall(PetscLayoutDestroy(&rLayout));
984: for (PetscInt field = 0; field < Nf; ++field) {
985: PetscObject obj;
986: PetscInt Nc;
988: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
989: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
990: totNc += Nc;
991: }
992: PetscCheck(totNc == 1, comm, PETSC_ERR_ARG_WRONG, "The number of field components %" PetscInt_FMT " != 1", totNc);
993: /* count non-zeros */
994: PetscCall(DMSwarmSortGetAccess(sw));
995: for (PetscInt field = 0; field < Nf; ++field) {
996: PetscObject obj;
997: PetscInt Nc;
999: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
1000: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
1001: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1002: PetscInt *pind;
1003: PetscInt Npc;
1005: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1006: maxNpc = PetscMax(maxNpc, Npc);
1007: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1008: }
1009: }
1010: PetscCall(PetscMalloc3(maxNpc * cdim * totDim, &elemMat, maxNpc * cdim, &rowIDXs, maxNpc * cdim, &xi));
1011: for (PetscInt field = 0; field < Nf; ++field) {
1012: PetscTabulation Tcoarse;
1013: PetscFE fe;
1015: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1016: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1017: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1018: PetscInt *findices, *pind;
1019: PetscInt numFIndices, Npc;
1021: /* TODO: Use DMField instead of assuming affine */
1022: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
1023: PetscCall(DMPlexGetClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1024: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1025: for (PetscInt j = 0; j < Npc; ++j) {
1026: PetscReal xr[8];
1027: PetscInt off = 0;
1029: for (PetscInt i = 0; i < Nfc; ++i) {
1030: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][pind[j] * bs[i] + b];
1031: }
1032: PetscCheck(off == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, cdim);
1033: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[j * cdim]);
1034: }
1035: PetscCall(PetscFECreateTabulation(fe, 1, Npc, xi, 1, &Tcoarse));
1036: /* Get elemMat entries by multiplying by weight */
1037: PetscCall(PetscArrayzero(elemMat, Npc * cdim * totDim));
1038: for (PetscInt i = 0; i < numFIndices; ++i) {
1039: for (PetscInt j = 0; j < Npc; ++j) {
1040: /* D[((p*pdim + i)*Nc + c)*cdim + d] is the value at point p for basis function i, component c, derviative d */
1041: for (PetscInt d = 0; d < cdim; ++d) {
1042: xi[d] = 0.;
1043: for (PetscInt e = 0; e < cdim; ++e) xi[d] += invJ[e * cdim + d] * Tcoarse->T[1][(j * numFIndices + i) * cdim + e];
1044: elemMat[(j * cdim + d) * numFIndices + i] += xi[d] * (useDeltaFunction ? 1.0 : detJ);
1045: }
1046: }
1047: }
1048: for (PetscInt j = 0; j < Npc; ++j)
1049: for (PetscInt d = 0; d < cdim; ++d) rowIDXs[j * cdim + d] = pind[j] * cdim + d + rStart;
1050: if (0) PetscCall(DMPrintCellMatrix(cell, name, Npc * cdim, numFIndices, elemMat));
1051: PetscCall(MatSetValues(derv, Npc * cdim, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
1052: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1053: PetscCall(DMPlexRestoreClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1054: PetscCall(PetscTabulationDestroy(&Tcoarse));
1055: }
1056: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1057: }
1058: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
1059: PetscCall(DMSwarmSortRestoreAccess(sw));
1060: PetscCall(PetscFree3(v0, J, invJ));
1061: PetscCall(PetscFree2(coordVals, bs));
1062: PetscCall(MatAssemblyBegin(derv, MAT_FINAL_ASSEMBLY));
1063: PetscCall(MatAssemblyEnd(derv, MAT_FINAL_ASSEMBLY));
1064: PetscFunctionReturn(PETSC_SUCCESS);
1065: }
1067: /* FEM cols: this is a scalar space
1068: Particle rows: this is a vector space that contracts with the derivative
1069: */
1070: static PetscErrorCode DMCreateGradientMatrix_Swarm(DM sw, DM dm, Mat *derv)
1071: {
1072: DMSwarmCellDM celldm;
1073: PetscSection gs;
1074: PetscInt cdim, m, n, Np, bs;
1075: void *ctx;
1076: MPI_Comm comm;
1078: PetscFunctionBegin;
1079: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1080: PetscCall(DMGetCoordinateDim(dm, &cdim));
1081: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1082: PetscCheck(celldm->Nf, comm, PETSC_ERR_USER, "Active cell DM does not define any fields");
1083: PetscCall(DMGetGlobalSection(dm, &gs));
1084: PetscCall(PetscSectionGetConstrainedStorageSize(gs, &n));
1085: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1086: PetscCall(DMSwarmCellDMGetBlockSize(celldm, sw, &bs));
1087: PetscCheck(cdim == bs, comm, PETSC_ERR_ARG_WRONG, "Coordinate dimension %" PetscInt_FMT " != %" PetscInt_FMT " swarm field block size", cdim, bs);
1088: m = Np * bs;
1089: PetscCall(MatCreate(PetscObjectComm((PetscObject)sw), derv));
1090: PetscCall(PetscObjectSetName((PetscObject)*derv, "Swarm Derivative Matrix"));
1091: PetscCall(MatSetSizes(*derv, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
1092: PetscCall(MatSetType(*derv, sw->mattype));
1093: PetscCall(DMGetApplicationContext(dm, &ctx));
1095: PetscCall(DMSwarmComputeGradientMatrix_Private(sw, dm, *derv, PETSC_TRUE, ctx));
1096: PetscCall(MatViewFromOptions(*derv, NULL, "-gradient_mat_view"));
1097: PetscFunctionReturn(PETSC_SUCCESS);
1098: }
1100: /*@
1101: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
1103: Collective
1105: Input Parameters:
1106: + dm - a `DMSWARM`
1107: - fieldname - the textual name given to a registered field
1109: Output Parameter:
1110: . vec - the vector
1112: Level: beginner
1114: Note:
1115: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
1117: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
1118: @*/
1119: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1120: {
1121: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1123: PetscFunctionBegin;
1125: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
1126: PetscFunctionReturn(PETSC_SUCCESS);
1127: }
1129: /*@
1130: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
1132: Collective
1134: Input Parameters:
1135: + dm - a `DMSWARM`
1136: - fieldname - the textual name given to a registered field
1138: Output Parameter:
1139: . vec - the vector
1141: Level: beginner
1143: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1144: @*/
1145: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1146: {
1147: PetscFunctionBegin;
1149: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1150: PetscFunctionReturn(PETSC_SUCCESS);
1151: }
1153: /*@
1154: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
1156: Collective
1158: Input Parameters:
1159: + dm - a `DMSWARM`
1160: - fieldname - the textual name given to a registered field
1162: Output Parameter:
1163: . vec - the vector
1165: Level: beginner
1167: Note:
1168: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1170: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1171: @*/
1172: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1173: {
1174: MPI_Comm comm = PETSC_COMM_SELF;
1176: PetscFunctionBegin;
1177: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
1178: PetscFunctionReturn(PETSC_SUCCESS);
1179: }
1181: /*@
1182: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
1184: Collective
1186: Input Parameters:
1187: + dm - a `DMSWARM`
1188: - fieldname - the textual name given to a registered field
1190: Output Parameter:
1191: . vec - the vector
1193: Level: beginner
1195: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
1196: @*/
1197: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1198: {
1199: PetscFunctionBegin;
1201: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: /*@
1206: DMSwarmCreateGlobalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1208: Collective
1210: Input Parameters:
1211: + dm - a `DMSWARM`
1212: . Nf - the number of fields
1213: - fieldnames - the textual names given to the registered fields
1215: Output Parameter:
1216: . vec - the vector
1218: Level: beginner
1220: Notes:
1221: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromFields()`.
1223: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateGlobalVectorFromField()`
1225: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromFields()`
1226: @*/
1227: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1228: {
1229: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1231: PetscFunctionBegin;
1233: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1234: PetscFunctionReturn(PETSC_SUCCESS);
1235: }
1237: /*@
1238: DMSwarmDestroyGlobalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1240: Collective
1242: Input Parameters:
1243: + dm - a `DMSWARM`
1244: . Nf - the number of fields
1245: - fieldnames - the textual names given to the registered fields
1247: Output Parameter:
1248: . vec - the vector
1250: Level: beginner
1252: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1253: @*/
1254: PetscErrorCode DMSwarmDestroyGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1255: {
1256: PetscFunctionBegin;
1258: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1259: PetscFunctionReturn(PETSC_SUCCESS);
1260: }
1262: /*@
1263: DMSwarmCreateLocalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1265: Collective
1267: Input Parameters:
1268: + dm - a `DMSWARM`
1269: . Nf - the number of fields
1270: - fieldnames - the textual names given to the registered fields
1272: Output Parameter:
1273: . vec - the vector
1275: Level: beginner
1277: Notes:
1278: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1280: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateLocalVectorFromField()`
1282: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1283: @*/
1284: PetscErrorCode DMSwarmCreateLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1285: {
1286: MPI_Comm comm = PETSC_COMM_SELF;
1288: PetscFunctionBegin;
1289: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1290: PetscFunctionReturn(PETSC_SUCCESS);
1291: }
1293: /*@
1294: DMSwarmDestroyLocalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1296: Collective
1298: Input Parameters:
1299: + dm - a `DMSWARM`
1300: . Nf - the number of fields
1301: - fieldnames - the textual names given to the registered fields
1303: Output Parameter:
1304: . vec - the vector
1306: Level: beginner
1308: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromFields()`
1309: @*/
1310: PetscErrorCode DMSwarmDestroyLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1311: {
1312: PetscFunctionBegin;
1314: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: /*@
1319: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
1321: Collective
1323: Input Parameter:
1324: . dm - a `DMSWARM`
1326: Level: beginner
1328: Note:
1329: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
1331: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1332: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1333: @*/
1334: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
1335: {
1336: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1338: PetscFunctionBegin;
1339: if (!swarm->field_registration_initialized) {
1340: swarm->field_registration_initialized = PETSC_TRUE;
1341: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
1342: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
1343: }
1344: PetscFunctionReturn(PETSC_SUCCESS);
1345: }
1347: /*@
1348: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
1350: Collective
1352: Input Parameter:
1353: . dm - a `DMSWARM`
1355: Level: beginner
1357: Note:
1358: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
1360: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1361: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1362: @*/
1363: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
1364: {
1365: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1367: PetscFunctionBegin;
1368: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
1369: swarm->field_registration_finalized = PETSC_TRUE;
1370: PetscFunctionReturn(PETSC_SUCCESS);
1371: }
1373: /*@
1374: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
1376: Not Collective
1378: Input Parameters:
1379: + sw - a `DMSWARM`
1380: . nlocal - the length of each registered field
1381: - buffer - the length of the buffer used to efficient dynamic re-sizing
1383: Level: beginner
1385: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
1386: @*/
1387: PetscErrorCode DMSwarmSetLocalSizes(DM sw, PetscInt nlocal, PetscInt buffer)
1388: {
1389: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1390: PetscMPIInt rank;
1391: PetscInt *rankval;
1393: PetscFunctionBegin;
1394: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
1395: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
1396: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
1398: // Initialize values in pid and rank placeholders
1399: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1400: PetscCall(DMSwarmGetField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1401: for (PetscInt p = 0; p < nlocal; p++) rankval[p] = rank;
1402: PetscCall(DMSwarmRestoreField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1403: /* TODO: [pid - use MPI_Scan] */
1404: PetscFunctionReturn(PETSC_SUCCESS);
1405: }
1407: /*@
1408: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
1410: Collective
1412: Input Parameters:
1413: + sw - a `DMSWARM`
1414: - dm - the `DM` to attach to the `DMSWARM`
1416: Level: beginner
1418: Note:
1419: The attached `DM` (dm) will be queried for point location and
1420: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1422: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
1423: @*/
1424: PetscErrorCode DMSwarmSetCellDM(DM sw, DM dm)
1425: {
1426: DMSwarmCellDM celldm;
1427: const char *name;
1428: char *coordName;
1430: PetscFunctionBegin;
1433: PetscCall(PetscStrallocpy(DMSwarmPICField_coor, &coordName));
1434: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 1, (const char **)&coordName, &celldm));
1435: PetscCall(PetscFree(coordName));
1436: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1437: PetscCall(DMSwarmAddCellDM(sw, celldm));
1438: PetscCall(DMSwarmCellDMDestroy(&celldm));
1439: PetscCall(DMSwarmSetCellDMActive(sw, name));
1440: PetscFunctionReturn(PETSC_SUCCESS);
1441: }
1443: /*@
1444: DMSwarmGetCellDM - Fetches the active cell `DM`
1446: Collective
1448: Input Parameter:
1449: . sw - a `DMSWARM`
1451: Output Parameter:
1452: . dm - the active `DM` for the `DMSWARM`
1454: Level: beginner
1456: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1457: @*/
1458: PetscErrorCode DMSwarmGetCellDM(DM sw, DM *dm)
1459: {
1460: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1461: DMSwarmCellDM celldm;
1463: PetscFunctionBegin;
1465: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)&celldm));
1466: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "There is no cell DM named %s in this Swarm", swarm->activeCellDM);
1467: PetscCall(DMSwarmCellDMGetDM(celldm, dm));
1468: PetscFunctionReturn(PETSC_SUCCESS);
1469: }
1471: /*@C
1472: DMSwarmGetCellDMNames - Get the list of cell `DM` names
1474: Not collective
1476: Input Parameter:
1477: . sw - a `DMSWARM`
1479: Output Parameters:
1480: + Ndm - the number of `DMSwarmCellDM` in the `DMSWARM`
1481: - celldms - the name of each `DMSwarmCellDM`
1483: Level: beginner
1485: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMByName()`
1486: @*/
1487: PetscErrorCode DMSwarmGetCellDMNames(DM sw, PetscInt *Ndm, const char **celldms[])
1488: {
1489: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1490: PetscObjectList next = swarm->cellDMs;
1491: PetscInt n = 0;
1493: PetscFunctionBegin;
1495: PetscAssertPointer(Ndm, 2);
1496: PetscAssertPointer(celldms, 3);
1497: while (next) {
1498: next = next->next;
1499: ++n;
1500: }
1501: PetscCall(PetscMalloc1(n, celldms));
1502: next = swarm->cellDMs;
1503: n = 0;
1504: while (next) {
1505: (*celldms)[n] = (const char *)next->obj->name;
1506: next = next->next;
1507: ++n;
1508: }
1509: *Ndm = n;
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }
1513: /*@
1514: DMSwarmSetCellDMActive - Activates a cell `DM` for a `DMSWARM`
1516: Collective
1518: Input Parameters:
1519: + sw - a `DMSWARM`
1520: - name - name of the cell `DM` to active for the `DMSWARM`
1522: Level: beginner
1524: Note:
1525: The attached `DM` (dmcell) will be queried for point location and
1526: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1528: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1529: @*/
1530: PetscErrorCode DMSwarmSetCellDMActive(DM sw, const char name[])
1531: {
1532: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1533: DMSwarmCellDM celldm;
1535: PetscFunctionBegin;
1537: PetscCall(PetscInfo(sw, "Setting cell DM to %s\n", name));
1538: PetscCall(PetscFree(swarm->activeCellDM));
1539: PetscCall(PetscStrallocpy(name, (char **)&swarm->activeCellDM));
1540: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1541: PetscFunctionReturn(PETSC_SUCCESS);
1542: }
1544: /*@
1545: DMSwarmGetCellDMActive - Returns the active cell `DM` for a `DMSWARM`
1547: Collective
1549: Input Parameter:
1550: . sw - a `DMSWARM`
1552: Output Parameter:
1553: . celldm - the active `DMSwarmCellDM`
1555: Level: beginner
1557: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1558: @*/
1559: PetscErrorCode DMSwarmGetCellDMActive(DM sw, DMSwarmCellDM *celldm)
1560: {
1561: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1563: PetscFunctionBegin;
1565: PetscAssertPointer(celldm, 2);
1566: PetscCheck(swarm->activeCellDM, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no active cell DM");
1567: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)celldm));
1568: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", swarm->activeCellDM);
1569: PetscFunctionReturn(PETSC_SUCCESS);
1570: }
1572: /*@C
1573: DMSwarmGetCellDMByName - Get a `DMSwarmCellDM` from its name
1575: Not collective
1577: Input Parameters:
1578: + sw - a `DMSWARM`
1579: - name - the name
1581: Output Parameter:
1582: . celldm - the `DMSwarmCellDM`
1584: Level: beginner
1586: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMNames()`
1587: @*/
1588: PetscErrorCode DMSwarmGetCellDMByName(DM sw, const char name[], DMSwarmCellDM *celldm)
1589: {
1590: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1592: PetscFunctionBegin;
1594: PetscAssertPointer(name, 2);
1595: PetscAssertPointer(celldm, 3);
1596: PetscCall(PetscObjectListFind(swarm->cellDMs, name, (PetscObject *)celldm));
1597: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", name);
1598: PetscFunctionReturn(PETSC_SUCCESS);
1599: }
1601: /*@
1602: DMSwarmAddCellDM - Adds a cell `DM` to the `DMSWARM`
1604: Collective
1606: Input Parameters:
1607: + sw - a `DMSWARM`
1608: - celldm - the `DMSwarmCellDM`
1610: Level: beginner
1612: Note:
1613: Cell DMs with the same name will share the cellid field
1615: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPushCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1616: @*/
1617: PetscErrorCode DMSwarmAddCellDM(DM sw, DMSwarmCellDM celldm)
1618: {
1619: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1620: const char *name;
1621: PetscInt dim;
1622: PetscBool flg;
1623: MPI_Comm comm;
1625: PetscFunctionBegin;
1627: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1629: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1630: PetscCall(PetscObjectListAdd(&swarm->cellDMs, name, (PetscObject)celldm));
1631: PetscCall(DMGetDimension(sw, &dim));
1632: for (PetscInt f = 0; f < celldm->Nfc; ++f) {
1633: PetscCall(DMSwarmDataFieldStringInList(celldm->coordFields[f], swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1634: if (!flg) {
1635: PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->coordFields[f], dim, PETSC_DOUBLE));
1636: } else {
1637: PetscDataType dt;
1638: PetscInt bs;
1640: PetscCall(DMSwarmGetFieldInfo(sw, celldm->coordFields[f], &bs, &dt));
1641: PetscCheck(bs == dim, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has blocksize %" PetscInt_FMT " != %" PetscInt_FMT " spatial dimension", celldm->coordFields[f], bs, dim);
1642: PetscCheck(dt == PETSC_DOUBLE, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has datatype %s != PETSC_DOUBLE", celldm->coordFields[f], PetscDataTypes[dt]);
1643: }
1644: }
1645: // Assume that DMs with the same name share the cellid field
1646: PetscCall(DMSwarmDataFieldStringInList(celldm->cellid, swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1647: if (!flg) {
1648: PetscBool isShell, isDummy;
1649: const char *name;
1651: // Allow dummy DMSHELL (I don't think we should support this mode)
1652: PetscCall(PetscObjectTypeCompare((PetscObject)celldm->dm, DMSHELL, &isShell));
1653: PetscCall(PetscObjectGetName((PetscObject)celldm->dm, &name));
1654: PetscCall(PetscStrcmp(name, "dummy", &isDummy));
1655: if (!isShell || !isDummy) PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->cellid, 1, PETSC_INT));
1656: }
1657: PetscCall(DMSwarmSetCellDMActive(sw, name));
1658: PetscFunctionReturn(PETSC_SUCCESS);
1659: }
1661: /*@
1662: DMSwarmGetLocalSize - Retrieves the local length of fields registered
1664: Not Collective
1666: Input Parameter:
1667: . dm - a `DMSWARM`
1669: Output Parameter:
1670: . nlocal - the length of each registered field
1672: Level: beginner
1674: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
1675: @*/
1676: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
1677: {
1678: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1680: PetscFunctionBegin;
1681: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
1682: PetscFunctionReturn(PETSC_SUCCESS);
1683: }
1685: /*@
1686: DMSwarmGetSize - Retrieves the total length of fields registered
1688: Collective
1690: Input Parameter:
1691: . dm - a `DMSWARM`
1693: Output Parameter:
1694: . n - the total length of each registered field
1696: Level: beginner
1698: Note:
1699: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
1701: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
1702: @*/
1703: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
1704: {
1705: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1706: PetscInt nlocal;
1708: PetscFunctionBegin;
1709: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1710: PetscCallMPI(MPIU_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1711: PetscFunctionReturn(PETSC_SUCCESS);
1712: }
1714: /*@C
1715: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
1717: Collective
1719: Input Parameters:
1720: + dm - a `DMSWARM`
1721: . fieldname - the textual name to identify this field
1722: . blocksize - the number of each data type
1723: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
1725: Level: beginner
1727: Notes:
1728: The textual name for each registered field must be unique.
1730: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1731: @*/
1732: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1733: {
1734: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1735: size_t size;
1737: PetscFunctionBegin;
1738: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1739: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1741: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1742: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1743: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1744: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1745: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1747: PetscCall(PetscDataTypeGetSize(type, &size));
1748: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1749: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1750: {
1751: DMSwarmDataField gfield;
1753: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1754: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1755: }
1756: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1757: PetscFunctionReturn(PETSC_SUCCESS);
1758: }
1760: /*@C
1761: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
1763: Collective
1765: Input Parameters:
1766: + dm - a `DMSWARM`
1767: . fieldname - the textual name to identify this field
1768: - size - the size in bytes of the user struct of each data type
1770: Level: beginner
1772: Note:
1773: The textual name for each registered field must be unique.
1775: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1776: @*/
1777: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1778: {
1779: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1781: PetscFunctionBegin;
1782: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1783: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1784: PetscFunctionReturn(PETSC_SUCCESS);
1785: }
1787: /*@C
1788: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1790: Collective
1792: Input Parameters:
1793: + dm - a `DMSWARM`
1794: . fieldname - the textual name to identify this field
1795: . size - the size in bytes of the user data type
1796: - blocksize - the number of each data type
1798: Level: beginner
1800: Note:
1801: The textual name for each registered field must be unique.
1803: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1804: @*/
1805: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1806: {
1807: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1809: PetscFunctionBegin;
1810: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1811: {
1812: DMSwarmDataField gfield;
1814: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1815: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1816: }
1817: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1818: PetscFunctionReturn(PETSC_SUCCESS);
1819: }
1821: /*@C
1822: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1824: Not Collective, No Fortran Support
1826: Input Parameters:
1827: + dm - a `DMSWARM`
1828: - fieldname - the textual name to identify this field
1830: Output Parameters:
1831: + blocksize - the number of each data type
1832: . type - the data type
1833: - data - pointer to raw array
1835: Level: beginner
1837: Notes:
1838: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1840: Fortran Note:
1841: Only works for `type` of `PETSC_SCALAR`
1843: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1844: @*/
1845: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1846: {
1847: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1848: DMSwarmDataField gfield;
1850: PetscFunctionBegin;
1852: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1853: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1854: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1855: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1856: if (blocksize) *blocksize = gfield->bs;
1857: if (type) *type = gfield->petsc_type;
1858: PetscFunctionReturn(PETSC_SUCCESS);
1859: }
1861: /*@C
1862: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1864: Not Collective
1866: Input Parameters:
1867: + dm - a `DMSWARM`
1868: - fieldname - the textual name to identify this field
1870: Output Parameters:
1871: + blocksize - the number of each data type
1872: . type - the data type
1873: - data - pointer to raw array
1875: Level: beginner
1877: Notes:
1878: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1880: Fortran Note:
1881: Only works for `type` of `PETSC_SCALAR`
1883: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1884: @*/
1885: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1886: {
1887: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1888: DMSwarmDataField gfield;
1890: PetscFunctionBegin;
1892: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1893: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1894: if (data) *data = NULL;
1895: PetscFunctionReturn(PETSC_SUCCESS);
1896: }
1898: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
1899: {
1900: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1901: DMSwarmDataField gfield;
1903: PetscFunctionBegin;
1905: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1906: if (blocksize) *blocksize = gfield->bs;
1907: if (type) *type = gfield->petsc_type;
1908: PetscFunctionReturn(PETSC_SUCCESS);
1909: }
1911: /*@
1912: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
1914: Not Collective
1916: Input Parameter:
1917: . dm - a `DMSWARM`
1919: Level: beginner
1921: Notes:
1922: The new point will have all fields initialized to zero.
1924: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
1925: @*/
1926: PetscErrorCode DMSwarmAddPoint(DM dm)
1927: {
1928: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1930: PetscFunctionBegin;
1931: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1932: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1933: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
1934: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1938: /*@
1939: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
1941: Not Collective
1943: Input Parameters:
1944: + dm - a `DMSWARM`
1945: - npoints - the number of new points to add
1947: Level: beginner
1949: Notes:
1950: The new point will have all fields initialized to zero.
1952: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
1953: @*/
1954: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1955: {
1956: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1957: PetscInt nlocal;
1959: PetscFunctionBegin;
1960: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1961: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1962: nlocal = PetscMax(nlocal, 0) + npoints;
1963: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
1964: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1965: PetscFunctionReturn(PETSC_SUCCESS);
1966: }
1968: /*@
1969: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
1971: Not Collective
1973: Input Parameter:
1974: . dm - a `DMSWARM`
1976: Level: beginner
1978: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
1979: @*/
1980: PetscErrorCode DMSwarmRemovePoint(DM dm)
1981: {
1982: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1984: PetscFunctionBegin;
1985: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1986: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
1987: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1988: PetscFunctionReturn(PETSC_SUCCESS);
1989: }
1991: /*@
1992: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
1994: Not Collective
1996: Input Parameters:
1997: + dm - a `DMSWARM`
1998: - idx - index of point to remove
2000: Level: beginner
2002: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2003: @*/
2004: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
2005: {
2006: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2008: PetscFunctionBegin;
2009: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
2010: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
2011: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
2012: PetscFunctionReturn(PETSC_SUCCESS);
2013: }
2015: /*@
2016: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
2018: Not Collective
2020: Input Parameters:
2021: + dm - a `DMSWARM`
2022: . pi - the index of the point to copy
2023: - pj - the point index where the copy should be located
2025: Level: beginner
2027: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2028: @*/
2029: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
2030: {
2031: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2033: PetscFunctionBegin;
2034: if (!swarm->issetup) PetscCall(DMSetUp(dm));
2035: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
2036: PetscFunctionReturn(PETSC_SUCCESS);
2037: }
2039: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
2040: {
2041: PetscFunctionBegin;
2042: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
2043: PetscFunctionReturn(PETSC_SUCCESS);
2044: }
2046: /*@
2047: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
2049: Collective
2051: Input Parameters:
2052: + dm - the `DMSWARM`
2053: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
2055: Level: advanced
2057: Notes:
2058: The `DM` will be modified to accommodate received points.
2059: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
2060: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
2062: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
2063: @*/
2064: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
2065: {
2066: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2068: PetscFunctionBegin;
2069: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
2070: switch (swarm->migrate_type) {
2071: case DMSWARM_MIGRATE_BASIC:
2072: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
2073: break;
2074: case DMSWARM_MIGRATE_DMCELLNSCATTER:
2075: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
2076: break;
2077: case DMSWARM_MIGRATE_DMCELLEXACT:
2078: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
2079: case DMSWARM_MIGRATE_USER:
2080: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
2081: default:
2082: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
2083: }
2084: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
2085: PetscCall(DMClearGlobalVectors(dm));
2086: PetscFunctionReturn(PETSC_SUCCESS);
2087: }
2089: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
2091: /*
2092: DMSwarmCollectViewCreate
2094: * Applies a collection method and gathers point neighbour points into dm
2096: Notes:
2097: Users should call DMSwarmCollectViewDestroy() after
2098: they have finished computations associated with the collected points
2099: */
2101: /*@
2102: DMSwarmCollectViewCreate - Applies a collection method and gathers points
2103: in neighbour ranks into the `DMSWARM`
2105: Collective
2107: Input Parameter:
2108: . dm - the `DMSWARM`
2110: Level: advanced
2112: Notes:
2113: Users should call `DMSwarmCollectViewDestroy()` after
2114: they have finished computations associated with the collected points
2116: Different collect methods are supported. See `DMSwarmSetCollectType()`.
2118: Developer Note:
2119: Create and Destroy routines create new objects that can get destroyed, they do not change the state
2120: of the current object.
2122: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
2123: @*/
2124: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
2125: {
2126: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2127: PetscInt ng;
2129: PetscFunctionBegin;
2130: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
2131: PetscCall(DMSwarmGetLocalSize(dm, &ng));
2132: switch (swarm->collect_type) {
2133: case DMSWARM_COLLECT_BASIC:
2134: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
2135: break;
2136: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
2137: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
2138: case DMSWARM_COLLECT_GENERAL:
2139: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
2140: default:
2141: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
2142: }
2143: swarm->collect_view_active = PETSC_TRUE;
2144: swarm->collect_view_reset_nlocal = ng;
2145: PetscFunctionReturn(PETSC_SUCCESS);
2146: }
2148: /*@
2149: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
2151: Collective
2153: Input Parameters:
2154: . dm - the `DMSWARM`
2156: Notes:
2157: Users should call `DMSwarmCollectViewCreate()` before this function is called.
2159: Level: advanced
2161: Developer Note:
2162: Create and Destroy routines create new objects that can get destroyed, they do not change the state
2163: of the current object.
2165: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
2166: @*/
2167: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
2168: {
2169: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2171: PetscFunctionBegin;
2172: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
2173: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
2174: swarm->collect_view_active = PETSC_FALSE;
2175: PetscFunctionReturn(PETSC_SUCCESS);
2176: }
2178: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
2179: {
2180: PetscInt dim;
2182: PetscFunctionBegin;
2183: PetscCall(DMSwarmSetNumSpecies(dm, 1));
2184: PetscCall(DMGetDimension(dm, &dim));
2185: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2186: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2187: PetscFunctionReturn(PETSC_SUCCESS);
2188: }
2190: /*@
2191: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
2193: Collective
2195: Input Parameters:
2196: + dm - the `DMSWARM`
2197: - Npc - The number of particles per cell in the cell `DM`
2199: Level: intermediate
2201: Notes:
2202: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
2203: one particle is in each cell, it is placed at the centroid.
2205: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
2206: @*/
2207: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
2208: {
2209: DM cdm;
2210: DMSwarmCellDM celldm;
2211: PetscRandom rnd;
2212: DMPolytopeType ct;
2213: PetscBool simplex;
2214: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
2215: PetscInt dim, d, cStart, cEnd, c, p, Nfc;
2216: const char **coordFields;
2218: PetscFunctionBeginUser;
2219: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
2220: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
2221: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
2223: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2224: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2225: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2226: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2227: PetscCall(DMGetDimension(cdm, &dim));
2228: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
2229: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
2230: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
2232: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
2233: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
2234: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2235: for (c = cStart; c < cEnd; ++c) {
2236: if (Npc == 1) {
2237: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
2238: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
2239: } else {
2240: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
2241: for (p = 0; p < Npc; ++p) {
2242: const PetscInt n = c * Npc + p;
2243: PetscReal sum = 0.0, refcoords[3];
2245: for (d = 0; d < dim; ++d) {
2246: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
2247: sum += refcoords[d];
2248: }
2249: if (simplex && sum > 0.0)
2250: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
2251: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
2252: }
2253: }
2254: }
2255: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2256: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
2257: PetscCall(PetscRandomDestroy(&rnd));
2258: PetscFunctionReturn(PETSC_SUCCESS);
2259: }
2261: /*@
2262: DMSwarmGetType - Get particular flavor of `DMSWARM`
2264: Collective
2266: Input Parameter:
2267: . sw - the `DMSWARM`
2269: Output Parameter:
2270: . stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2272: Level: advanced
2274: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2275: @*/
2276: PetscErrorCode DMSwarmGetType(DM sw, DMSwarmType *stype)
2277: {
2278: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2280: PetscFunctionBegin;
2282: PetscAssertPointer(stype, 2);
2283: *stype = swarm->swarm_type;
2284: PetscFunctionReturn(PETSC_SUCCESS);
2285: }
2287: /*@
2288: DMSwarmSetType - Set particular flavor of `DMSWARM`
2290: Collective
2292: Input Parameters:
2293: + sw - the `DMSWARM`
2294: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2296: Level: advanced
2298: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2299: @*/
2300: PetscErrorCode DMSwarmSetType(DM sw, DMSwarmType stype)
2301: {
2302: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2304: PetscFunctionBegin;
2306: swarm->swarm_type = stype;
2307: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(sw));
2308: PetscFunctionReturn(PETSC_SUCCESS);
2309: }
2311: static PetscErrorCode DMSwarmCreateRemapDM_Private(DM sw, DM *rdm)
2312: {
2313: PetscFE fe;
2314: DMPolytopeType ct;
2315: PetscInt dim, cStart;
2316: const char *prefix = "remap_";
2318: PetscFunctionBegin;
2319: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), rdm));
2320: PetscCall(DMSetType(*rdm, DMPLEX));
2321: PetscCall(DMPlexSetOptionsPrefix(*rdm, prefix));
2322: PetscCall(DMSetFromOptions(*rdm));
2323: PetscCall(PetscObjectSetName((PetscObject)*rdm, "remap"));
2324: PetscCall(DMViewFromOptions(*rdm, NULL, "-dm_view"));
2326: PetscCall(DMGetDimension(*rdm, &dim));
2327: PetscCall(DMPlexGetHeightStratum(*rdm, 0, &cStart, NULL));
2328: PetscCall(DMPlexGetCellType(*rdm, cStart, &ct));
2329: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
2330: PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
2331: PetscCall(DMSetField(*rdm, 0, NULL, (PetscObject)fe));
2332: PetscCall(DMCreateDS(*rdm));
2333: PetscCall(PetscFEDestroy(&fe));
2334: PetscFunctionReturn(PETSC_SUCCESS);
2335: }
2337: static PetscErrorCode DMSetup_Swarm(DM sw)
2338: {
2339: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2341: PetscFunctionBegin;
2342: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
2343: swarm->issetup = PETSC_TRUE;
2345: if (swarm->remap_type != DMSWARM_REMAP_NONE) {
2346: DMSwarmCellDM celldm;
2347: DM rdm;
2348: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
2349: const char *vfieldnames[1] = {"w_q"};
2351: PetscCall(DMSwarmCreateRemapDM_Private(sw, &rdm));
2352: PetscCall(DMSwarmCellDMCreate(rdm, 1, vfieldnames, 2, fieldnames, &celldm));
2353: PetscCall(DMSwarmAddCellDM(sw, celldm));
2354: PetscCall(DMSwarmCellDMDestroy(&celldm));
2355: PetscCall(DMDestroy(&rdm));
2356: }
2358: if (swarm->swarm_type == DMSWARM_PIC) {
2359: DMSwarmCellDM celldm;
2361: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2362: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "No active cell DM. DMSWARM_PIC requires you call DMSwarmSetCellDM() or DMSwarmAddCellDM()");
2363: if (celldm->dm->ops->locatepointssubdomain) {
2364: /* check methods exists for exact ownership identificiation */
2365: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
2366: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
2367: } else {
2368: /* check methods exist for point location AND rank neighbor identification */
2369: PetscCheck(celldm->dm->ops->locatepoints, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
2370: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
2372: PetscCheck(celldm->dm->ops->getneighbors, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
2373: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
2375: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
2376: }
2377: }
2379: PetscCall(DMSwarmFinalizeFieldRegister(sw));
2381: /* check some fields were registered */
2382: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
2383: PetscFunctionReturn(PETSC_SUCCESS);
2384: }
2386: static PetscErrorCode DMDestroy_Swarm(DM dm)
2387: {
2388: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2390: PetscFunctionBegin;
2391: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2392: PetscCall(PetscObjectListDestroy(&swarm->cellDMs));
2393: PetscCall(PetscFree(swarm->activeCellDM));
2394: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
2395: PetscCall(PetscFree(swarm));
2396: PetscFunctionReturn(PETSC_SUCCESS);
2397: }
2399: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
2400: {
2401: DM cdm;
2402: DMSwarmCellDM celldm;
2403: PetscDraw draw;
2404: PetscReal *coords, oldPause, radius = 0.01;
2405: PetscInt Np, p, bs, Nfc;
2406: const char **coordFields;
2408: PetscFunctionBegin;
2409: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
2410: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2411: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2412: PetscCall(PetscDrawGetPause(draw, &oldPause));
2413: PetscCall(PetscDrawSetPause(draw, 0.0));
2414: PetscCall(DMView(cdm, viewer));
2415: PetscCall(PetscDrawSetPause(draw, oldPause));
2417: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2418: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2419: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2420: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2421: PetscCall(DMSwarmGetField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2422: for (p = 0; p < Np; ++p) {
2423: const PetscInt i = p * bs;
2425: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
2426: }
2427: PetscCall(DMSwarmRestoreField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2428: PetscCall(PetscDrawFlush(draw));
2429: PetscCall(PetscDrawPause(draw));
2430: PetscCall(PetscDrawSave(draw));
2431: PetscFunctionReturn(PETSC_SUCCESS);
2432: }
2434: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
2435: {
2436: PetscViewerFormat format;
2437: PetscInt *sizes;
2438: PetscInt dim, Np, maxSize = 17;
2439: MPI_Comm comm;
2440: PetscMPIInt rank, size, p;
2441: const char *name, *cellid;
2443: PetscFunctionBegin;
2444: PetscCall(PetscViewerGetFormat(viewer, &format));
2445: PetscCall(DMGetDimension(dm, &dim));
2446: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2447: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2448: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2449: PetscCallMPI(MPI_Comm_size(comm, &size));
2450: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
2451: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
2452: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
2453: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
2454: else PetscCall(PetscCalloc1(3, &sizes));
2455: if (size < maxSize) {
2456: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
2457: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
2458: for (p = 0; p < size; ++p) {
2459: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
2460: }
2461: } else {
2462: PetscInt locMinMax[2] = {Np, Np};
2464: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
2465: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
2466: }
2467: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2468: PetscCall(PetscFree(sizes));
2469: if (format == PETSC_VIEWER_ASCII_INFO) {
2470: DMSwarmCellDM celldm;
2471: PetscInt *cell;
2473: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
2474: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2475: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2476: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
2477: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&cell));
2478: for (p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%d: %" PetscInt_FMT "\n", p, cell[p]));
2479: PetscCall(PetscViewerFlush(viewer));
2480: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2481: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&cell));
2482: }
2483: PetscFunctionReturn(PETSC_SUCCESS);
2484: }
2486: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
2487: {
2488: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2489: PetscBool isascii, ibinary, isvtk, isdraw, ispython;
2490: #if defined(PETSC_HAVE_HDF5)
2491: PetscBool ishdf5;
2492: #endif
2494: PetscFunctionBegin;
2497: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2498: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
2499: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
2500: #if defined(PETSC_HAVE_HDF5)
2501: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2502: #endif
2503: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2504: PetscCall(PetscObjectHasFunction((PetscObject)viewer, "PetscViewerPythonViewObject_C", &ispython));
2505: if (isascii) {
2506: PetscViewerFormat format;
2508: PetscCall(PetscViewerGetFormat(viewer, &format));
2509: switch (format) {
2510: case PETSC_VIEWER_ASCII_INFO_DETAIL:
2511: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
2512: break;
2513: default:
2514: PetscCall(DMView_Swarm_Ascii(dm, viewer));
2515: }
2516: } else {
2517: #if defined(PETSC_HAVE_HDF5)
2518: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
2519: #endif
2520: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
2521: if (ispython) PetscCall(PetscViewerPythonViewObject(viewer, (PetscObject)dm));
2522: }
2523: PetscFunctionReturn(PETSC_SUCCESS);
2524: }
2526: /*@
2527: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
2528: 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.
2530: Noncollective
2532: Input Parameters:
2533: + sw - the `DMSWARM`
2534: . cellID - the integer id of the cell to be extracted and filtered
2535: - cellswarm - The `DMSWARM` to receive the cell
2537: Level: beginner
2539: Notes:
2540: This presently only supports `DMSWARM_PIC` type
2542: Should be restored with `DMSwarmRestoreCellSwarm()`
2544: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
2546: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
2547: @*/
2548: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2549: {
2550: DM_Swarm *original = (DM_Swarm *)sw->data;
2551: DMLabel label;
2552: DM dmc, subdmc;
2553: PetscInt *pids, particles, dim;
2554: const char *name;
2556: PetscFunctionBegin;
2557: /* Configure new swarm */
2558: PetscCall(DMSetType(cellswarm, DMSWARM));
2559: PetscCall(DMGetDimension(sw, &dim));
2560: PetscCall(DMSetDimension(cellswarm, dim));
2561: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
2562: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
2563: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
2564: PetscCall(DMSwarmSortGetAccess(sw));
2565: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
2566: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2567: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
2568: PetscCall(DMSwarmSortRestoreAccess(sw));
2569: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2570: PetscCall(DMSwarmGetCellDM(sw, &dmc));
2571: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
2572: PetscCall(DMAddLabel(dmc, label));
2573: PetscCall(DMLabelSetValue(label, cellID, 1));
2574: PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdmc));
2575: PetscCall(PetscObjectGetName((PetscObject)dmc, &name));
2576: PetscCall(PetscObjectSetName((PetscObject)subdmc, name));
2577: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
2578: PetscCall(DMLabelDestroy(&label));
2579: PetscFunctionReturn(PETSC_SUCCESS);
2580: }
2582: /*@
2583: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
2585: Noncollective
2587: Input Parameters:
2588: + sw - the parent `DMSWARM`
2589: . cellID - the integer id of the cell to be copied back into the parent swarm
2590: - cellswarm - the cell swarm object
2592: Level: beginner
2594: Note:
2595: This only supports `DMSWARM_PIC` types of `DMSWARM`s
2597: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
2598: @*/
2599: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2600: {
2601: DM dmc;
2602: PetscInt *pids, particles, p;
2604: PetscFunctionBegin;
2605: PetscCall(DMSwarmSortGetAccess(sw));
2606: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2607: PetscCall(DMSwarmSortRestoreAccess(sw));
2608: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
2609: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
2610: /* Free memory, destroy cell dm */
2611: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
2612: PetscCall(DMDestroy(&dmc));
2613: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2614: PetscFunctionReturn(PETSC_SUCCESS);
2615: }
2617: /*@
2618: DMSwarmComputeMoments - Compute the first three particle moments for a given field
2620: Noncollective
2622: Input Parameters:
2623: + sw - the `DMSWARM`
2624: . coordinate - the coordinate field name
2625: - weight - the weight field name
2627: Output Parameter:
2628: . moments - the field moments
2630: Level: intermediate
2632: Notes:
2633: The `moments` array should be of length bs + 2, where bs is the block size of the coordinate field.
2635: The weight field must be a scalar, having blocksize 1.
2637: .seealso: `DM`, `DMSWARM`, `DMPlexComputeMoments()`
2638: @*/
2639: PetscErrorCode DMSwarmComputeMoments(DM sw, const char coordinate[], const char weight[], PetscReal moments[])
2640: {
2641: const PetscReal *coords;
2642: const PetscReal *w;
2643: PetscReal *mom;
2644: PetscDataType dtc, dtw;
2645: PetscInt bsc, bsw, Np;
2646: MPI_Comm comm;
2648: PetscFunctionBegin;
2650: PetscAssertPointer(coordinate, 2);
2651: PetscAssertPointer(weight, 3);
2652: PetscAssertPointer(moments, 4);
2653: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
2654: PetscCall(DMSwarmGetField(sw, coordinate, &bsc, &dtc, (void **)&coords));
2655: PetscCall(DMSwarmGetField(sw, weight, &bsw, &dtw, (void **)&w));
2656: PetscCheck(dtc == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s must be real, not %s", coordinate, PetscDataTypes[dtc]);
2657: PetscCheck(dtw == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be real, not %s", weight, PetscDataTypes[dtw]);
2658: PetscCheck(bsw == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be a scalar, not blocksize %" PetscInt_FMT, weight, bsw);
2659: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2660: PetscCall(DMGetWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2661: PetscCall(PetscArrayzero(mom, bsc + 2));
2662: for (PetscInt p = 0; p < Np; ++p) {
2663: const PetscReal *c = &coords[p * bsc];
2664: const PetscReal wp = w[p];
2666: mom[0] += wp;
2667: for (PetscInt d = 0; d < bsc; ++d) {
2668: mom[d + 1] += wp * c[d];
2669: mom[d + bsc + 1] += wp * PetscSqr(c[d]);
2670: }
2671: }
2672: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&coords));
2673: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
2674: PetscCallMPI(MPIU_Allreduce(mom, moments, bsc + 2, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
2675: PetscCall(DMRestoreWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2676: PetscFunctionReturn(PETSC_SUCCESS);
2677: }
2679: static PetscErrorCode DMSetFromOptions_Swarm(DM dm, PetscOptionItems PetscOptionsObject)
2680: {
2681: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2683: PetscFunctionBegin;
2684: PetscOptionsHeadBegin(PetscOptionsObject, "DMSwarm Options");
2685: PetscCall(PetscOptionsEnum("-dm_swarm_remap_type", "Remap algorithm", "DMSwarmSetRemapType", DMSwarmRemapTypeNames, (PetscEnum)swarm->remap_type, (PetscEnum *)&swarm->remap_type, NULL));
2686: PetscOptionsHeadEnd();
2687: PetscFunctionReturn(PETSC_SUCCESS);
2688: }
2690: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
2692: static PetscErrorCode DMInitialize_Swarm(DM sw)
2693: {
2694: PetscFunctionBegin;
2695: sw->ops->view = DMView_Swarm;
2696: sw->ops->load = NULL;
2697: sw->ops->setfromoptions = DMSetFromOptions_Swarm;
2698: sw->ops->clone = DMClone_Swarm;
2699: sw->ops->setup = DMSetup_Swarm;
2700: sw->ops->createlocalsection = NULL;
2701: sw->ops->createsectionpermutation = NULL;
2702: sw->ops->createdefaultconstraints = NULL;
2703: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
2704: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
2705: sw->ops->getlocaltoglobalmapping = NULL;
2706: sw->ops->createfieldis = NULL;
2707: sw->ops->createcoordinatedm = NULL;
2708: sw->ops->createcellcoordinatedm = NULL;
2709: sw->ops->getcoloring = NULL;
2710: sw->ops->creatematrix = DMCreateMatrix_Swarm;
2711: sw->ops->createinterpolation = NULL;
2712: sw->ops->createinjection = NULL;
2713: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
2714: sw->ops->creategradientmatrix = DMCreateGradientMatrix_Swarm;
2715: sw->ops->refine = NULL;
2716: sw->ops->coarsen = NULL;
2717: sw->ops->refinehierarchy = NULL;
2718: sw->ops->coarsenhierarchy = NULL;
2719: sw->ops->globaltolocalbegin = DMGlobalToLocalBegin_Swarm;
2720: sw->ops->globaltolocalend = DMGlobalToLocalEnd_Swarm;
2721: sw->ops->localtoglobalbegin = DMLocalToGlobalBegin_Swarm;
2722: sw->ops->localtoglobalend = DMLocalToGlobalEnd_Swarm;
2723: sw->ops->destroy = DMDestroy_Swarm;
2724: sw->ops->createsubdm = NULL;
2725: sw->ops->getdimpoints = NULL;
2726: sw->ops->locatepoints = NULL;
2727: sw->ops->projectfieldlocal = DMProjectFieldLocal_Swarm;
2728: PetscFunctionReturn(PETSC_SUCCESS);
2729: }
2731: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
2732: {
2733: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2735: PetscFunctionBegin;
2736: swarm->refct++;
2737: (*newdm)->data = swarm;
2738: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
2739: PetscCall(DMInitialize_Swarm(*newdm));
2740: (*newdm)->dim = dm->dim;
2741: PetscFunctionReturn(PETSC_SUCCESS);
2742: }
2744: /*MC
2745: DMSWARM = "swarm" - A `DM` object for particle methods, such as particle-in-cell (PIC), in which the underlying
2746: data is both (i) dynamic in length, (ii) and of arbitrary data type.
2748: Level: intermediate
2750: Notes:
2751: User data can be represented by `DMSWARM` through a registering "fields" which are to be stored on particles.
2752: To register a field, the user must provide;
2753: (a) a unique name;
2754: (b) the data type (or size in bytes);
2755: (c) the block size of the data.
2757: For example, suppose the application requires a unique id, energy, momentum and density to be stored
2758: on a set of particles. Then the following code could be used
2759: .vb
2760: DMSwarmInitializeFieldRegister(dm)
2761: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
2762: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
2763: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
2764: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
2765: DMSwarmFinalizeFieldRegister(dm)
2766: .ve
2768: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
2769: The only restriction imposed by `DMSWARM` is that all fields contain the same number of particles.
2771: To support particle methods, "migration" techniques are provided. These methods migrate data
2772: between ranks.
2774: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
2775: As a `DMSWARM` may internally define and store values of different data types,
2776: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
2777: fields should be used to define a `Vec` object via `DMSwarmVectorDefineField()`
2778: The specified field can be changed at any time - thereby permitting vectors
2779: compatible with different fields to be created.
2781: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via `DMSwarmCreateGlobalVectorFromField()`
2782: Here the data defining the field in the `DMSWARM` is shared with a `Vec`.
2783: This is inherently unsafe if you alter the size of the field at any time between
2784: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
2785: If the local size of the `DMSWARM` does not match the local size of the global vector
2786: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
2788: Additional high-level support is provided for Particle-In-Cell methods. Refer to `DMSwarmSetType()`.
2790: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`, `DMSwarmSetType()`, `DMSwarmType`, `DMSwarmCreateGlobalVectorFromField()`,
2791: `DMCreateGlobalVector()`, `DMCreateLocalVector()`
2792: M*/
2794: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
2795: {
2796: DM_Swarm *swarm;
2798: PetscFunctionBegin;
2800: PetscCall(PetscNew(&swarm));
2801: dm->data = swarm;
2802: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
2803: PetscCall(DMSwarmInitializeFieldRegister(dm));
2804: dm->dim = 0;
2805: swarm->refct = 1;
2806: swarm->issetup = PETSC_FALSE;
2807: swarm->swarm_type = DMSWARM_BASIC;
2808: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
2809: swarm->collect_type = DMSWARM_COLLECT_BASIC;
2810: swarm->migrate_error_on_missing_point = PETSC_FALSE;
2811: swarm->collect_view_active = PETSC_FALSE;
2812: swarm->collect_view_reset_nlocal = -1;
2813: PetscCall(DMInitialize_Swarm(dm));
2814: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
2815: PetscFunctionReturn(PETSC_SUCCESS);
2816: }
2818: /* Replace dm with the contents of ndm, and then destroy ndm
2819: - Share the DM_Swarm structure
2820: */
2821: PetscErrorCode DMSwarmReplace(DM dm, DM *ndm)
2822: {
2823: DM dmNew = *ndm;
2824: const PetscReal *maxCell, *Lstart, *L;
2825: PetscInt dim;
2827: PetscFunctionBegin;
2828: if (dm == dmNew) {
2829: PetscCall(DMDestroy(ndm));
2830: PetscFunctionReturn(PETSC_SUCCESS);
2831: }
2832: dm->setupcalled = dmNew->setupcalled;
2833: if (!dm->hdr.name) {
2834: const char *name;
2836: PetscCall(PetscObjectGetName((PetscObject)*ndm, &name));
2837: PetscCall(PetscObjectSetName((PetscObject)dm, name));
2838: }
2839: PetscCall(DMGetDimension(dmNew, &dim));
2840: PetscCall(DMSetDimension(dm, dim));
2841: PetscCall(DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L));
2842: PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
2843: PetscCall(DMDestroy_Swarm(dm));
2844: PetscCall(DMInitialize_Swarm(dm));
2845: dm->data = dmNew->data;
2846: ((DM_Swarm *)dmNew->data)->refct++;
2847: PetscCall(DMDestroy(ndm));
2848: PetscFunctionReturn(PETSC_SUCCESS);
2849: }
2851: /*@
2852: DMSwarmDuplicate - Creates a new `DMSWARM` with the same fields and cell `DM`s but no particles
2854: Collective
2856: Input Parameter:
2857: . sw - the `DMSWARM`
2859: Output Parameter:
2860: . nsw - the new `DMSWARM`
2862: Level: beginner
2864: .seealso: `DM`, `DMSWARM`, `DMSwarmCreate()`, `DMClone()`
2865: @*/
2866: PetscErrorCode DMSwarmDuplicate(DM sw, DM *nsw)
2867: {
2868: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2869: DMSwarmDataField *fields;
2870: DMSwarmCellDM celldm, ncelldm;
2871: DMSwarmType stype;
2872: const char *name, **celldmnames;
2873: void *ctx;
2874: PetscInt dim, Nf, Ndm;
2875: PetscBool flg;
2877: PetscFunctionBegin;
2878: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), nsw));
2879: PetscCall(DMSetType(*nsw, DMSWARM));
2880: PetscCall(PetscObjectGetName((PetscObject)sw, &name));
2881: PetscCall(PetscObjectSetName((PetscObject)*nsw, name));
2882: PetscCall(DMGetDimension(sw, &dim));
2883: PetscCall(DMSetDimension(*nsw, dim));
2884: PetscCall(DMSwarmGetType(sw, &stype));
2885: PetscCall(DMSwarmSetType(*nsw, stype));
2886: PetscCall(DMGetApplicationContext(sw, &ctx));
2887: PetscCall(DMSetApplicationContext(*nsw, ctx));
2889: PetscCall(DMSwarmDataBucketGetDMSwarmDataFields(swarm->db, &Nf, &fields));
2890: for (PetscInt f = 0; f < Nf; ++f) {
2891: PetscCall(DMSwarmDataFieldStringInList(fields[f]->name, ((DM_Swarm *)(*nsw)->data)->db->nfields, (const DMSwarmDataField *)((DM_Swarm *)(*nsw)->data)->db->field, &flg));
2892: if (!flg) PetscCall(DMSwarmRegisterPetscDatatypeField(*nsw, fields[f]->name, fields[f]->bs, fields[f]->petsc_type));
2893: }
2895: PetscCall(DMSwarmGetCellDMNames(sw, &Ndm, &celldmnames));
2896: for (PetscInt c = 0; c < Ndm; ++c) {
2897: DM dm;
2898: PetscInt Ncf;
2899: const char **coordfields, **fields;
2901: PetscCall(DMSwarmGetCellDMByName(sw, celldmnames[c], &celldm));
2902: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
2903: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncf, &coordfields));
2904: PetscCall(DMSwarmCellDMGetFields(celldm, &Nf, &fields));
2905: PetscCall(DMSwarmCellDMCreate(dm, Nf, fields, Ncf, coordfields, &ncelldm));
2906: PetscCall(DMSwarmAddCellDM(*nsw, ncelldm));
2907: PetscCall(DMSwarmCellDMDestroy(&ncelldm));
2908: }
2909: PetscCall(PetscFree(celldmnames));
2911: PetscCall(DMSetFromOptions(*nsw));
2912: PetscCall(DMSetUp(*nsw));
2913: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2914: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
2915: PetscCall(DMSwarmSetCellDMActive(*nsw, name));
2916: PetscFunctionReturn(PETSC_SUCCESS);
2917: }
2919: PetscErrorCode DMLocalToGlobalBegin_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
2920: {
2921: PetscFunctionBegin;
2922: PetscFunctionReturn(PETSC_SUCCESS);
2923: }
2925: PetscErrorCode DMLocalToGlobalEnd_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
2926: {
2927: PetscFunctionBegin;
2928: switch (mode) {
2929: case INSERT_VALUES:
2930: PetscCall(VecCopy(l, g));
2931: break;
2932: case ADD_VALUES:
2933: PetscCall(VecAXPY(g, 1., l));
2934: break;
2935: default:
2936: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mode not supported: %d", mode);
2937: }
2938: PetscFunctionReturn(PETSC_SUCCESS);
2939: }
2941: PetscErrorCode DMGlobalToLocalBegin_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
2942: {
2943: PetscFunctionBegin;
2944: PetscFunctionReturn(PETSC_SUCCESS);
2945: }
2947: PetscErrorCode DMGlobalToLocalEnd_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
2948: {
2949: PetscFunctionBegin;
2950: PetscCall(DMLocalToGlobalEnd_Swarm(dm, g, mode, l));
2951: PetscFunctionReturn(PETSC_SUCCESS);
2952: }