Actual source code: swarm.c
1: #include "petscdmswarm.h"
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petsc/private/hashsetij.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petscviewer.h>
6: #include <petscdraw.h>
7: #include <petscdmplex.h>
8: #include <petscblaslapack.h>
9: #include "../src/dm/impls/swarm/data_bucket.h"
10: #include <petscdmlabel.h>
11: #include <petscsection.h>
13: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
14: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
15: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
17: const char *DMSwarmTypeNames[] = {"basic", "pic", NULL};
18: const char *DMSwarmMigrateTypeNames[] = {"basic", "dmcellnscatter", "dmcellexact", "user", NULL};
19: const char *DMSwarmCollectTypeNames[] = {"basic", "boundingbox", "general", "user", NULL};
20: const char *DMSwarmRemapTypeNames[] = {"none", "pfak", "colella", "DMSwarmRemapType", "DMSWARM_REMAP_", NULL};
21: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
23: const char DMSwarmField_pid[] = "DMSwarm_pid";
24: const char DMSwarmField_rank[] = "DMSwarm_rank";
25: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
27: PetscInt SwarmDataFieldId = -1;
29: #if defined(PETSC_HAVE_HDF5)
30: #include <petscviewerhdf5.h>
32: static PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
33: {
34: DM dm;
35: PetscReal seqval;
36: PetscInt seqnum, bs;
37: PetscBool isseq, ists;
39: PetscFunctionBegin;
40: PetscCall(VecGetDM(v, &dm));
41: PetscCall(VecGetBlockSize(v, &bs));
42: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particle_fields"));
43: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
44: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
45: PetscCall(PetscViewerHDF5IsTimestepping(viewer, &ists));
46: if (ists) PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
47: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
48: PetscCall(VecViewNative(v, viewer));
49: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
50: PetscCall(PetscViewerHDF5PopGroup(viewer));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
55: {
56: DMSwarmCellDM celldm;
57: Vec coordinates;
58: PetscInt Np, Nfc;
59: PetscBool isseq;
60: const char **coordFields;
62: PetscFunctionBegin;
63: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
64: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
65: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
66: PetscCall(DMSwarmGetSize(dm, &Np));
67: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, coordFields[0], &coordinates));
68: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
69: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
70: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
71: PetscCall(VecViewNative(coordinates, viewer));
72: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
73: PetscCall(PetscViewerHDF5PopGroup(viewer));
74: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, coordFields[0], &coordinates));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
77: #endif
79: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
80: {
81: DM dm;
82: #if defined(PETSC_HAVE_HDF5)
83: PetscBool ishdf5;
84: #endif
86: PetscFunctionBegin;
87: PetscCall(VecGetDM(v, &dm));
88: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
89: #if defined(PETSC_HAVE_HDF5)
90: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
91: if (ishdf5) {
92: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
95: #endif
96: PetscCall(VecViewNative(v, viewer));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: /*@C
101: DMSwarmVectorGetField - Gets the fields from which to define a `Vec` object
102: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
104: Not collective
106: Input Parameter:
107: . sw - a `DMSWARM`
109: Output Parameters:
110: + Nf - the number of fields
111: - fieldnames - the textual name given to each registered field, or NULL if it has not been set
113: Level: beginner
115: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
116: @*/
117: PetscErrorCode DMSwarmVectorGetField(DM sw, PetscInt *Nf, const char **fieldnames[])
118: {
119: DMSwarmCellDM celldm;
121: PetscFunctionBegin;
123: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
124: PetscCall(DMSwarmCellDMGetFields(celldm, Nf, fieldnames));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*@
129: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
130: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
132: Collective
134: Input Parameters:
135: + dm - a `DMSWARM`
136: - fieldname - the textual name given to each registered field
138: Level: beginner
140: Notes:
141: The field with name `fieldname` must be defined as having a data type of `PetscScalar`.
143: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
144: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
146: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineFields()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
147: @*/
148: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
149: {
150: PetscFunctionBegin;
151: PetscCall(DMSwarmVectorDefineFields(dm, 1, &fieldname));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: /*@C
156: DMSwarmVectorDefineFields - Sets the fields from which to define a `Vec` object
157: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
159: Collective, No Fortran support
161: Input Parameters:
162: + sw - a `DMSWARM`
163: . Nf - the number of fields
164: - fieldnames - the textual name given to each registered field
166: Level: beginner
168: Notes:
169: Each field with name in `fieldnames` must be defined as having a data type of `PetscScalar`.
171: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
172: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
174: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
175: @*/
176: PetscErrorCode DMSwarmVectorDefineFields(DM sw, PetscInt Nf, const char *fieldnames[])
177: {
178: DM_Swarm *swarm = (DM_Swarm *)sw->data;
179: DMSwarmCellDM celldm;
181: PetscFunctionBegin;
183: if (fieldnames) PetscAssertPointer(fieldnames, 3);
184: if (!swarm->issetup) PetscCall(DMSetUp(sw));
185: PetscCheck(Nf >= 0, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields must be non-negative, not %" PetscInt_FMT, Nf);
186: // Create a dummy cell DM if none has been specified (I think we should not support this mode)
187: if (!swarm->activeCellDM) {
188: DM dm;
189: DMSwarmCellDM celldm;
191: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), &dm));
192: PetscCall(DMSetType(dm, DMSHELL));
193: PetscCall(PetscObjectSetName((PetscObject)dm, "dummy"));
194: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 0, NULL, &celldm));
195: PetscCall(DMDestroy(&dm));
196: PetscCall(DMSwarmAddCellDM(sw, celldm));
197: PetscCall(DMSwarmCellDMDestroy(&celldm));
198: PetscCall(DMSwarmSetCellDMActive(sw, "dummy"));
199: }
200: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
201: for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscFree(celldm->dmFields[f]));
202: PetscCall(PetscFree(celldm->dmFields));
204: celldm->Nf = Nf;
205: PetscCall(PetscMalloc1(Nf, &celldm->dmFields));
206: for (PetscInt f = 0; f < Nf; ++f) {
207: PetscDataType type;
209: // Check all fields are of type PETSC_REAL or PETSC_SCALAR
210: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], NULL, &type));
211: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
212: PetscCall(PetscStrallocpy(fieldnames[f], (char **)&celldm->dmFields[f]));
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /* requires DMSwarmDefineFieldVector has been called */
218: static PetscErrorCode DMCreateGlobalVector_Swarm(DM sw, Vec *vec)
219: {
220: DM_Swarm *swarm = (DM_Swarm *)sw->data;
221: DMSwarmCellDM celldm;
222: Vec x;
223: char name[PETSC_MAX_PATH_LEN];
224: PetscInt bs = 0, n;
226: PetscFunctionBegin;
227: if (!swarm->issetup) PetscCall(DMSetUp(sw));
228: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
229: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
230: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
232: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
233: for (PetscInt f = 0; f < celldm->Nf; ++f) {
234: PetscInt fbs;
235: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
236: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
237: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
238: bs += fbs;
239: }
240: PetscCall(VecCreate(PetscObjectComm((PetscObject)sw), &x));
241: PetscCall(PetscObjectSetName((PetscObject)x, name));
242: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
243: PetscCall(VecSetBlockSize(x, bs));
244: PetscCall(VecSetDM(x, sw));
245: PetscCall(VecSetFromOptions(x));
246: PetscCall(VecSetOperation(x, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
247: *vec = x;
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: /* requires DMSwarmDefineFieldVector has been called */
252: static PetscErrorCode DMCreateLocalVector_Swarm(DM sw, Vec *vec)
253: {
254: DM_Swarm *swarm = (DM_Swarm *)sw->data;
255: DMSwarmCellDM celldm;
256: Vec x;
257: char name[PETSC_MAX_PATH_LEN];
258: PetscInt bs = 0, n;
260: PetscFunctionBegin;
261: if (!swarm->issetup) PetscCall(DMSetUp(sw));
262: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
263: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
264: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
266: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
267: for (PetscInt f = 0; f < celldm->Nf; ++f) {
268: PetscInt fbs;
269: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
270: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
271: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
272: bs += fbs;
273: }
274: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
275: PetscCall(PetscObjectSetName((PetscObject)x, name));
276: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
277: PetscCall(VecSetBlockSize(x, bs));
278: PetscCall(VecSetDM(x, sw));
279: PetscCall(VecSetFromOptions(x));
280: *vec = x;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
285: {
286: DM_Swarm *swarm = (DM_Swarm *)dm->data;
287: DMSwarmDataField gfield;
288: PetscInt bs, nlocal, fid = -1, cfid = -2;
289: PetscBool flg;
291: PetscFunctionBegin;
292: /* check vector is an inplace array */
293: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
294: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
295: (void)flg; /* avoid compiler warning */
296: 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);
297: PetscCall(VecGetLocalSize(*vec, &nlocal));
298: PetscCall(VecGetBlockSize(*vec, &bs));
299: 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");
300: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
301: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
302: PetscCall(VecResetArray(*vec));
303: PetscCall(VecDestroy(vec));
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
308: {
309: DM_Swarm *swarm = (DM_Swarm *)dm->data;
310: PetscDataType type;
311: PetscScalar *array;
312: PetscInt bs, n, fid;
313: char name[PETSC_MAX_PATH_LEN];
314: PetscMPIInt size;
315: PetscBool iscuda, iskokkos, iship;
317: PetscFunctionBegin;
318: if (!swarm->issetup) PetscCall(DMSetUp(dm));
319: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
320: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
321: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
323: PetscCallMPI(MPI_Comm_size(comm, &size));
324: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
325: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
326: PetscCall(PetscStrcmp(dm->vectype, VECHIP, &iship));
327: PetscCall(VecCreate(comm, vec));
328: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
329: PetscCall(VecSetBlockSize(*vec, bs));
330: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
331: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
332: else if (iship) PetscCall(VecSetType(*vec, VECHIP));
333: else PetscCall(VecSetType(*vec, VECSTANDARD));
334: PetscCall(VecPlaceArray(*vec, array));
336: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
337: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
339: /* Set guard */
340: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
341: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
343: PetscCall(VecSetDM(*vec, dm));
344: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: static PetscErrorCode DMSwarmDestroyVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], Vec *vec)
349: {
350: DM_Swarm *swarm = (DM_Swarm *)sw->data;
351: const PetscScalar *array;
352: PetscInt bs, n, id = 0, cid = -2;
353: PetscBool flg;
355: PetscFunctionBegin;
356: for (PetscInt f = 0; f < Nf; ++f) {
357: PetscInt fid;
359: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
360: id += fid;
361: }
362: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cid, flg));
363: (void)flg; /* avoid compiler warning */
364: 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);
365: PetscCall(VecGetLocalSize(*vec, &n));
366: PetscCall(VecGetBlockSize(*vec, &bs));
367: n /= bs;
368: PetscCheck(n == swarm->db->L, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
369: PetscCall(VecGetArrayRead(*vec, &array));
370: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
371: PetscScalar *farray;
372: PetscDataType ftype;
373: PetscInt fbs;
375: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
376: PetscCheck(off + fbs <= bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid blocksize %" PetscInt_FMT " < %" PetscInt_FMT, bs, off + fbs);
377: for (PetscInt i = 0; i < n; ++i) {
378: for (PetscInt b = 0; b < fbs; ++b) farray[i * fbs + b] = array[i * bs + off + b];
379: }
380: off += fbs;
381: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
382: }
383: PetscCall(VecRestoreArrayRead(*vec, &array));
384: PetscCall(VecDestroy(vec));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: static PetscErrorCode DMSwarmCreateVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], MPI_Comm comm, Vec *vec)
389: {
390: DM_Swarm *swarm = (DM_Swarm *)sw->data;
391: PetscScalar *array;
392: PetscInt n, bs = 0, id = 0;
393: char name[PETSC_MAX_PATH_LEN];
395: PetscFunctionBegin;
396: if (!swarm->issetup) PetscCall(DMSetUp(sw));
397: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
398: for (PetscInt f = 0; f < Nf; ++f) {
399: PetscDataType ftype;
400: PetscInt fbs;
402: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
403: PetscCheck(ftype == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
404: bs += fbs;
405: }
407: PetscCall(VecCreate(comm, vec));
408: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
409: PetscCall(VecSetBlockSize(*vec, bs));
410: PetscCall(VecSetType(*vec, sw->vectype));
412: PetscCall(VecGetArrayWrite(*vec, &array));
413: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
414: PetscScalar *farray;
415: PetscDataType ftype;
416: PetscInt fbs;
418: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
419: for (PetscInt i = 0; i < n; ++i) {
420: for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = farray[i * fbs + b];
421: }
422: off += fbs;
423: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
424: }
425: PetscCall(VecRestoreArrayWrite(*vec, &array));
427: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
428: for (PetscInt f = 0; f < Nf; ++f) {
429: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
430: PetscCall(PetscStrlcat(name, fieldnames[f], PETSC_MAX_PATH_LEN));
431: }
432: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
434: for (PetscInt f = 0; f < Nf; ++f) {
435: PetscInt fid;
437: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
438: id += fid;
439: }
440: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, id));
442: PetscCall(VecSetDM(*vec, sw));
443: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
449: \hat f = \sum_i f_i \phi_i
451: and a particle function is given by
453: f = \sum_i w_i \delta(x - x_i)
455: then we want to require that
457: M \hat f = M_p f
459: where the particle mass matrix is given by
461: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
463: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
464: his integral. We allow this with the boolean flag.
465: */
466: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
467: {
468: const char *name = "Mass Matrix";
469: MPI_Comm comm;
470: DMSwarmCellDM celldm;
471: PetscDS prob;
472: PetscSection fsection, globalFSection;
473: PetscHSetIJ ht;
474: PetscLayout rLayout, colLayout;
475: PetscInt *dnz, *onz;
476: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
477: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
478: PetscScalar *elemMat;
479: PetscInt dim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0, totNc = 0;
480: const char **coordFields;
481: PetscReal **coordVals;
482: PetscInt *bs;
484: PetscFunctionBegin;
485: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
486: PetscCall(DMGetCoordinateDim(dmf, &dim));
487: PetscCall(DMGetDS(dmf, &prob));
488: PetscCall(PetscDSGetNumFields(prob, &Nf));
489: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
490: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
491: PetscCall(DMGetLocalSection(dmf, &fsection));
492: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
493: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
494: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
496: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
497: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
498: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
500: PetscCall(PetscLayoutCreate(comm, &colLayout));
501: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
502: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
503: PetscCall(PetscLayoutSetUp(colLayout));
504: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
505: PetscCall(PetscLayoutDestroy(&colLayout));
507: PetscCall(PetscLayoutCreate(comm, &rLayout));
508: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
509: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
510: PetscCall(PetscLayoutSetUp(rLayout));
511: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
512: PetscCall(PetscLayoutDestroy(&rLayout));
514: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
515: PetscCall(PetscHSetIJCreate(&ht));
517: PetscCall(PetscSynchronizedFlush(comm, NULL));
518: for (PetscInt field = 0; field < Nf; ++field) {
519: PetscObject obj;
520: PetscClassId id;
521: PetscInt Nc;
523: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
524: PetscCall(PetscObjectGetClassId(obj, &id));
525: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
526: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
527: totNc += Nc;
528: }
529: /* count non-zeros */
530: PetscCall(DMSwarmSortGetAccess(dmc));
531: for (PetscInt field = 0; field < Nf; ++field) {
532: PetscObject obj;
533: PetscClassId id;
534: PetscInt Nc;
536: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
537: PetscCall(PetscObjectGetClassId(obj, &id));
538: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
539: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
541: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
542: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
543: PetscInt numFIndices, numCIndices;
545: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
546: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
547: maxC = PetscMax(maxC, numCIndices);
548: {
549: PetscHashIJKey key;
550: PetscBool missing;
551: for (PetscInt i = 0; i < numFIndices; ++i) {
552: key.j = findices[i]; /* global column (from Plex) */
553: if (key.j >= 0) {
554: /* Get indices for coarse elements */
555: for (PetscInt j = 0; j < numCIndices; ++j) {
556: for (PetscInt c = 0; c < Nc; ++c) {
557: // TODO Need field offset on particle here
558: key.i = cindices[j] * totNc + c + rStart; /* global cols (from Swarm) */
559: if (key.i < 0) continue;
560: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
561: PetscCheck(missing, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
562: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
563: else ++onz[key.i - rStart];
564: }
565: }
566: }
567: }
568: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
569: }
570: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
571: }
572: }
573: PetscCall(PetscHSetIJDestroy(&ht));
574: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
575: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
576: PetscCall(PetscFree2(dnz, onz));
577: PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
578: for (PetscInt field = 0; field < Nf; ++field) {
579: PetscTabulation Tcoarse;
580: PetscObject obj;
581: PetscClassId id;
582: PetscInt Nc;
584: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
585: PetscCall(PetscObjectGetClassId(obj, &id));
586: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
587: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
589: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
590: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
591: PetscInt *findices, *cindices;
592: PetscInt numFIndices, numCIndices;
594: /* TODO: Use DMField instead of assuming affine */
595: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
596: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
597: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
598: for (PetscInt j = 0; j < numCIndices; ++j) {
599: PetscReal xr[8];
600: PetscInt off = 0;
602: for (PetscInt i = 0; i < Nfc; ++i) {
603: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[j] * bs[i] + b];
604: }
605: 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);
606: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, xr, &xi[j * dim]);
607: }
608: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
609: else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
610: /* Get elemMat entries by multiplying by weight */
611: PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
612: for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
613: for (PetscInt j = 0; j < numCIndices; ++j) {
614: for (PetscInt c = 0; c < Nc; ++c) {
615: // TODO Need field offset on particle and field here
616: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
617: elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
618: }
619: }
620: }
621: for (PetscInt j = 0; j < numCIndices; ++j)
622: // TODO Need field offset on particle here
623: for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
624: if (0) PetscCall(DMPrintCellMatrix(cell, name, numCIndices * Nc, numFIndices, elemMat));
625: PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
626: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
627: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
628: PetscCall(PetscTabulationDestroy(&Tcoarse));
629: }
630: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
631: }
632: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
633: PetscCall(DMSwarmSortRestoreAccess(dmc));
634: PetscCall(PetscFree3(v0, J, invJ));
635: PetscCall(PetscFree2(coordVals, bs));
636: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
637: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: /* Returns empty matrix for use with SNES FD */
642: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
643: {
644: Vec field;
645: PetscInt size;
647: PetscFunctionBegin;
648: PetscCall(DMGetGlobalVector(sw, &field));
649: PetscCall(VecGetLocalSize(field, &size));
650: PetscCall(DMRestoreGlobalVector(sw, &field));
651: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
652: PetscCall(MatSetFromOptions(*m));
653: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
654: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
655: PetscCall(MatZeroEntries(*m));
656: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
657: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
658: PetscCall(MatShift(*m, 1.0));
659: PetscCall(MatSetDM(*m, sw));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /* FEM cols, Particle rows */
664: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
665: {
666: DMSwarmCellDM celldm;
667: PetscSection gsf;
668: PetscInt m, n, Np, bs;
669: void *ctx;
671: PetscFunctionBegin;
672: PetscCall(DMSwarmGetCellDMActive(dmCoarse, &celldm));
673: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)dmCoarse), PETSC_ERR_USER, "Active cell DM does not define any fields");
674: PetscCall(DMGetGlobalSection(dmFine, &gsf));
675: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
676: PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
677: PetscCall(DMSwarmCellDMGetBlockSize(celldm, dmCoarse, &bs));
678: n = Np * bs;
679: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
680: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
681: PetscCall(MatSetType(*mass, dmCoarse->mattype));
682: PetscCall(DMGetApplicationContext(dmFine, &ctx));
684: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
685: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }
689: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
690: {
691: const char *name = "Mass Matrix Square";
692: MPI_Comm comm;
693: DMSwarmCellDM celldm;
694: PetscDS prob;
695: PetscSection fsection, globalFSection;
696: PetscHSetIJ ht;
697: PetscLayout rLayout, colLayout;
698: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
699: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
700: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
701: PetscScalar *elemMat, *elemMatSq;
702: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0;
703: const char **coordFields;
704: PetscReal **coordVals;
705: PetscInt *bs;
707: PetscFunctionBegin;
708: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
709: PetscCall(DMGetCoordinateDim(dmf, &cdim));
710: PetscCall(DMGetDS(dmf, &prob));
711: PetscCall(PetscDSGetNumFields(prob, &Nf));
712: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
713: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
714: PetscCall(DMGetLocalSection(dmf, &fsection));
715: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
716: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
717: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
719: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
720: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
721: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
723: PetscCall(PetscLayoutCreate(comm, &colLayout));
724: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
725: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
726: PetscCall(PetscLayoutSetUp(colLayout));
727: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
728: PetscCall(PetscLayoutDestroy(&colLayout));
730: PetscCall(PetscLayoutCreate(comm, &rLayout));
731: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
732: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
733: PetscCall(PetscLayoutSetUp(rLayout));
734: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
735: PetscCall(PetscLayoutDestroy(&rLayout));
737: PetscCall(DMPlexGetDepth(dmf, &depth));
738: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
739: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
740: PetscCall(PetscMalloc1(maxAdjSize, &adj));
742: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
743: PetscCall(PetscHSetIJCreate(&ht));
744: /* Count nonzeros
745: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
746: */
747: PetscCall(DMSwarmSortGetAccess(dmc));
748: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
749: PetscInt *cindices;
750: PetscInt numCIndices;
751: #if 0
752: PetscInt adjSize = maxAdjSize, a, j;
753: #endif
755: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
756: maxC = PetscMax(maxC, numCIndices);
757: /* Diagonal block */
758: for (PetscInt i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
759: #if 0
760: /* Off-diagonal blocks */
761: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
762: for (a = 0; a < adjSize; ++a) {
763: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
764: const PetscInt ncell = adj[a];
765: PetscInt *ncindices;
766: PetscInt numNCIndices;
768: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
769: {
770: PetscHashIJKey key;
771: PetscBool missing;
773: for (i = 0; i < numCIndices; ++i) {
774: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
775: if (key.i < 0) continue;
776: for (j = 0; j < numNCIndices; ++j) {
777: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
778: if (key.j < 0) continue;
779: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
780: if (missing) {
781: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
782: else ++onz[key.i - rStart];
783: }
784: }
785: }
786: }
787: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
788: }
789: }
790: #endif
791: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
792: }
793: PetscCall(PetscHSetIJDestroy(&ht));
794: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
795: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
796: PetscCall(PetscFree2(dnz, onz));
797: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
798: /* Fill in values
799: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
800: Start just by producing block diagonal
801: Could loop over adjacent cells
802: Produce neighboring element matrix
803: TODO Determine which columns and rows correspond to shared dual vector
804: Do MatMatMult with rectangular matrices
805: Insert block
806: */
807: for (PetscInt field = 0; field < Nf; ++field) {
808: PetscTabulation Tcoarse;
809: PetscObject obj;
810: PetscInt Nc;
812: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
813: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
814: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
815: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
816: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
817: PetscInt *findices, *cindices;
818: PetscInt numFIndices, numCIndices;
820: /* TODO: Use DMField instead of assuming affine */
821: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
822: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
823: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
824: for (PetscInt p = 0; p < numCIndices; ++p) {
825: PetscReal xr[8];
826: PetscInt off = 0;
828: for (PetscInt i = 0; i < Nfc; ++i) {
829: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[p] * bs[i] + b];
830: }
831: 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);
832: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[p * cdim]);
833: }
834: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
835: /* Get elemMat entries by multiplying by weight */
836: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
837: for (PetscInt i = 0; i < numFIndices; ++i) {
838: for (PetscInt p = 0; p < numCIndices; ++p) {
839: for (PetscInt c = 0; c < Nc; ++c) {
840: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
841: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
842: }
843: }
844: }
845: PetscCall(PetscTabulationDestroy(&Tcoarse));
846: for (PetscInt p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
847: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
848: /* Block diagonal */
849: if (numCIndices) {
850: PetscBLASInt blasn, blask;
851: PetscScalar one = 1.0, zero = 0.0;
853: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
854: PetscCall(PetscBLASIntCast(numFIndices, &blask));
855: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
856: }
857: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
858: /* TODO off-diagonal */
859: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
860: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
861: }
862: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
863: }
864: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
865: PetscCall(PetscFree(adj));
866: PetscCall(DMSwarmSortRestoreAccess(dmc));
867: PetscCall(PetscFree3(v0, J, invJ));
868: PetscCall(PetscFree2(coordVals, bs));
869: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
870: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: /*@
875: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
877: Collective
879: Input Parameters:
880: + dmCoarse - a `DMSWARM`
881: - dmFine - a `DMPLEX`
883: Output Parameter:
884: . mass - the square of the particle mass matrix
886: Level: advanced
888: Note:
889: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
890: future to compute the full normal equations.
892: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
893: @*/
894: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
895: {
896: PetscInt n;
897: void *ctx;
899: PetscFunctionBegin;
900: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
901: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
902: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
903: PetscCall(MatSetType(*mass, dmCoarse->mattype));
904: PetscCall(DMGetApplicationContext(dmFine, &ctx));
906: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
907: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
911: /* 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
913: \int_X \psi_i \nabla \cdot \hat f = \int_X \psi_i \nabla \cdot f
915: and then integrate by parts
917: \int_X \nabla \psi_i \cdot \hat f = \int_X \nabla \psi_i \cdot f
919: where \psi is from a scalar FE space. If a finite element interpolant is given by
921: \hat f^c = \sum_i f_i \phi^c_i
923: and a particle function is given by
925: f^c = \sum_p f^c_p \delta(x - x_p)
927: then we want to require that
929: D_f \hat f = D_p f
931: where the gradient matrices are given by
933: (D_f)_{i(jc)} = \int \partial_c \psi_i \phi_j
934: (D_p)_{i(jc)} = \int \partial_c \psi_i \delta(x - x_j)
936: Thus we need two finite element spaces, a scalar and a vector. The vector space holds the representer for the
937: vector particle field. The scalar space holds the output of D_p or D_f, which is the weak divergence of the field.
939: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
940: his integral. We allow this with the boolean flag.
941: */
942: static PetscErrorCode DMSwarmComputeGradientMatrix_Private(DM sw, DM dm, Mat derv, PetscBool useDeltaFunction, void *ctx)
943: {
944: const char *name = "Derivative Matrix";
945: MPI_Comm comm;
946: DMSwarmCellDM celldm;
947: PetscDS ds;
948: PetscSection fsection, globalFSection;
949: PetscLayout rLayout;
950: PetscInt locRows, rStart, *rowIDXs;
951: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
952: PetscScalar *elemMat;
953: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxNpc = 0, totNc = 0;
954: const char **coordFields;
955: PetscReal **coordVals;
956: PetscInt *bs;
958: PetscFunctionBegin;
959: PetscCall(PetscObjectGetComm((PetscObject)derv, &comm));
960: PetscCall(DMGetCoordinateDim(dm, &cdim));
961: PetscCall(DMGetDS(dm, &ds));
962: PetscCall(PetscDSGetNumFields(ds, &Nf));
963: PetscCheck(Nf == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
964: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
965: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
966: PetscCall(DMGetLocalSection(dm, &fsection));
967: PetscCall(DMGetGlobalSection(dm, &globalFSection));
968: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
969: PetscCall(MatGetLocalSize(derv, &locRows, NULL));
971: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
972: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
973: PetscCheck(Nfc == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
974: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
976: PetscCall(PetscLayoutCreate(comm, &rLayout));
977: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
978: PetscCall(PetscLayoutSetBlockSize(rLayout, cdim));
979: PetscCall(PetscLayoutSetUp(rLayout));
980: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
981: PetscCall(PetscLayoutDestroy(&rLayout));
983: for (PetscInt field = 0; field < Nf; ++field) {
984: PetscObject obj;
985: PetscInt Nc;
987: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
988: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
989: totNc += Nc;
990: }
991: PetscCheck(totNc == 1, comm, PETSC_ERR_ARG_WRONG, "The number of field components %" PetscInt_FMT " != 1", totNc);
992: /* count non-zeros */
993: PetscCall(DMSwarmSortGetAccess(sw));
994: for (PetscInt field = 0; field < Nf; ++field) {
995: PetscObject obj;
996: PetscInt Nc;
998: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
999: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
1000: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1001: PetscInt *pind;
1002: PetscInt Npc;
1004: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1005: maxNpc = PetscMax(maxNpc, Npc);
1006: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1007: }
1008: }
1009: PetscCall(PetscMalloc3(maxNpc * cdim * totDim, &elemMat, maxNpc * cdim, &rowIDXs, maxNpc * cdim, &xi));
1010: for (PetscInt field = 0; field < Nf; ++field) {
1011: PetscTabulation Tcoarse;
1012: PetscFE fe;
1014: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1015: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1016: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1017: PetscInt *findices, *pind;
1018: PetscInt numFIndices, Npc;
1020: /* TODO: Use DMField instead of assuming affine */
1021: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
1022: PetscCall(DMPlexGetClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1023: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1024: for (PetscInt j = 0; j < Npc; ++j) {
1025: PetscReal xr[8];
1026: PetscInt off = 0;
1028: for (PetscInt i = 0; i < Nfc; ++i) {
1029: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][pind[j] * bs[i] + b];
1030: }
1031: 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);
1032: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[j * cdim]);
1033: }
1034: PetscCall(PetscFECreateTabulation(fe, 1, Npc, xi, 1, &Tcoarse));
1035: /* Get elemMat entries by multiplying by weight */
1036: PetscCall(PetscArrayzero(elemMat, Npc * cdim * totDim));
1037: for (PetscInt i = 0; i < numFIndices; ++i) {
1038: for (PetscInt j = 0; j < Npc; ++j) {
1039: /* D[((p*pdim + i)*Nc + c)*cdim + d] is the value at point p for basis function i, component c, derivative d */
1040: for (PetscInt d = 0; d < cdim; ++d) {
1041: xi[d] = 0.;
1042: for (PetscInt e = 0; e < cdim; ++e) xi[d] += invJ[e * cdim + d] * Tcoarse->T[1][(j * numFIndices + i) * cdim + e];
1043: elemMat[(j * cdim + d) * numFIndices + i] += xi[d] * (useDeltaFunction ? 1.0 : detJ);
1044: }
1045: }
1046: }
1047: for (PetscInt j = 0; j < Npc; ++j)
1048: for (PetscInt d = 0; d < cdim; ++d) rowIDXs[j * cdim + d] = pind[j] * cdim + d + rStart;
1049: if (0) PetscCall(DMPrintCellMatrix(cell, name, Npc * cdim, numFIndices, elemMat));
1050: PetscCall(MatSetValues(derv, Npc * cdim, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
1051: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1052: PetscCall(DMPlexRestoreClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1053: PetscCall(PetscTabulationDestroy(&Tcoarse));
1054: }
1055: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1056: }
1057: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
1058: PetscCall(DMSwarmSortRestoreAccess(sw));
1059: PetscCall(PetscFree3(v0, J, invJ));
1060: PetscCall(PetscFree2(coordVals, bs));
1061: PetscCall(MatAssemblyBegin(derv, MAT_FINAL_ASSEMBLY));
1062: PetscCall(MatAssemblyEnd(derv, MAT_FINAL_ASSEMBLY));
1063: PetscFunctionReturn(PETSC_SUCCESS);
1064: }
1066: /* FEM cols: this is a scalar space
1067: Particle rows: this is a vector space that contracts with the derivative
1068: */
1069: static PetscErrorCode DMCreateGradientMatrix_Swarm(DM sw, DM dm, Mat *derv)
1070: {
1071: DMSwarmCellDM celldm;
1072: PetscSection gs;
1073: PetscInt cdim, m, n, Np, bs;
1074: void *ctx;
1075: MPI_Comm comm;
1077: PetscFunctionBegin;
1078: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1079: PetscCall(DMGetCoordinateDim(dm, &cdim));
1080: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1081: PetscCheck(celldm->Nf, comm, PETSC_ERR_USER, "Active cell DM does not define any fields");
1082: PetscCall(DMGetGlobalSection(dm, &gs));
1083: PetscCall(PetscSectionGetConstrainedStorageSize(gs, &n));
1084: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1085: PetscCall(DMSwarmCellDMGetBlockSize(celldm, sw, &bs));
1086: PetscCheck(cdim == bs, comm, PETSC_ERR_ARG_WRONG, "Coordinate dimension %" PetscInt_FMT " != %" PetscInt_FMT " swarm field block size", cdim, bs);
1087: m = Np * bs;
1088: PetscCall(MatCreate(PetscObjectComm((PetscObject)sw), derv));
1089: PetscCall(PetscObjectSetName((PetscObject)*derv, "Swarm Derivative Matrix"));
1090: PetscCall(MatSetSizes(*derv, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
1091: PetscCall(MatSetType(*derv, sw->mattype));
1092: PetscCall(DMGetApplicationContext(dm, &ctx));
1094: PetscCall(DMSwarmComputeGradientMatrix_Private(sw, dm, *derv, PETSC_TRUE, ctx));
1095: PetscCall(MatViewFromOptions(*derv, NULL, "-gradient_mat_view"));
1096: PetscFunctionReturn(PETSC_SUCCESS);
1097: }
1099: /*@
1100: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
1102: Collective
1104: Input Parameters:
1105: + dm - a `DMSWARM`
1106: - fieldname - the textual name given to a registered field
1108: Output Parameter:
1109: . vec - the vector
1111: Level: beginner
1113: Note:
1114: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
1116: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
1117: @*/
1118: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1119: {
1120: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1122: PetscFunctionBegin;
1124: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: /*@
1129: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
1131: Collective
1133: Input Parameters:
1134: + dm - a `DMSWARM`
1135: - fieldname - the textual name given to a registered field
1137: Output Parameter:
1138: . vec - the vector
1140: Level: beginner
1142: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1143: @*/
1144: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1145: {
1146: PetscFunctionBegin;
1148: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1149: PetscFunctionReturn(PETSC_SUCCESS);
1150: }
1152: /*@
1153: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
1155: Collective
1157: Input Parameters:
1158: + dm - a `DMSWARM`
1159: - fieldname - the textual name given to a registered field
1161: Output Parameter:
1162: . vec - the vector
1164: Level: beginner
1166: Note:
1167: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1169: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1170: @*/
1171: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1172: {
1173: MPI_Comm comm = PETSC_COMM_SELF;
1175: PetscFunctionBegin;
1176: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
1177: PetscFunctionReturn(PETSC_SUCCESS);
1178: }
1180: /*@
1181: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
1183: Collective
1185: Input Parameters:
1186: + dm - a `DMSWARM`
1187: - fieldname - the textual name given to a registered field
1189: Output Parameter:
1190: . vec - the vector
1192: Level: beginner
1194: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
1195: @*/
1196: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1197: {
1198: PetscFunctionBegin;
1200: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1201: PetscFunctionReturn(PETSC_SUCCESS);
1202: }
1204: /*@
1205: DMSwarmCreateGlobalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1207: Collective
1209: Input Parameters:
1210: + dm - a `DMSWARM`
1211: . Nf - the number of fields
1212: - fieldnames - the textual names given to the registered fields
1214: Output Parameter:
1215: . vec - the vector
1217: Level: beginner
1219: Notes:
1220: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromFields()`.
1222: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateGlobalVectorFromField()`
1224: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromFields()`
1225: @*/
1226: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1227: {
1228: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1230: PetscFunctionBegin;
1232: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1233: PetscFunctionReturn(PETSC_SUCCESS);
1234: }
1236: /*@
1237: DMSwarmDestroyGlobalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1239: Collective
1241: Input Parameters:
1242: + dm - a `DMSWARM`
1243: . Nf - the number of fields
1244: - fieldnames - the textual names given to the registered fields
1246: Output Parameter:
1247: . vec - the vector
1249: Level: beginner
1251: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1252: @*/
1253: PetscErrorCode DMSwarmDestroyGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1254: {
1255: PetscFunctionBegin;
1257: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1261: /*@
1262: DMSwarmCreateLocalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1264: Collective
1266: Input Parameters:
1267: + dm - a `DMSWARM`
1268: . Nf - the number of fields
1269: - fieldnames - the textual names given to the registered fields
1271: Output Parameter:
1272: . vec - the vector
1274: Level: beginner
1276: Notes:
1277: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1279: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateLocalVectorFromField()`
1281: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1282: @*/
1283: PetscErrorCode DMSwarmCreateLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1284: {
1285: MPI_Comm comm = PETSC_COMM_SELF;
1287: PetscFunctionBegin;
1288: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1289: PetscFunctionReturn(PETSC_SUCCESS);
1290: }
1292: /*@
1293: DMSwarmDestroyLocalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1295: Collective
1297: Input Parameters:
1298: + dm - a `DMSWARM`
1299: . Nf - the number of fields
1300: - fieldnames - the textual names given to the registered fields
1302: Output Parameter:
1303: . vec - the vector
1305: Level: beginner
1307: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromFields()`
1308: @*/
1309: PetscErrorCode DMSwarmDestroyLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1310: {
1311: PetscFunctionBegin;
1313: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1314: PetscFunctionReturn(PETSC_SUCCESS);
1315: }
1317: /*@
1318: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
1320: Collective
1322: Input Parameter:
1323: . dm - a `DMSWARM`
1325: Level: beginner
1327: Note:
1328: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
1330: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1331: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1332: @*/
1333: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
1334: {
1335: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1337: PetscFunctionBegin;
1338: if (!swarm->field_registration_initialized) {
1339: swarm->field_registration_initialized = PETSC_TRUE;
1340: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
1341: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
1342: }
1343: PetscFunctionReturn(PETSC_SUCCESS);
1344: }
1346: /*@
1347: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
1349: Collective
1351: Input Parameter:
1352: . dm - a `DMSWARM`
1354: Level: beginner
1356: Note:
1357: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
1359: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1360: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1361: @*/
1362: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
1363: {
1364: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1366: PetscFunctionBegin;
1367: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
1368: swarm->field_registration_finalized = PETSC_TRUE;
1369: PetscFunctionReturn(PETSC_SUCCESS);
1370: }
1372: /*@
1373: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
1375: Not Collective
1377: Input Parameters:
1378: + sw - a `DMSWARM`
1379: . nlocal - the length of each registered field
1380: - buffer - the length of the buffer used to efficient dynamic re-sizing
1382: Level: beginner
1384: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
1385: @*/
1386: PetscErrorCode DMSwarmSetLocalSizes(DM sw, PetscInt nlocal, PetscInt buffer)
1387: {
1388: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1389: PetscMPIInt rank;
1390: PetscInt *rankval;
1392: PetscFunctionBegin;
1393: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
1394: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
1395: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
1397: // Initialize values in pid and rank placeholders
1398: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1399: PetscCall(DMSwarmGetField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1400: for (PetscInt p = 0; p < nlocal; p++) rankval[p] = rank;
1401: PetscCall(DMSwarmRestoreField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1402: /* TODO: [pid - use MPI_Scan] */
1403: PetscFunctionReturn(PETSC_SUCCESS);
1404: }
1406: /*@
1407: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
1409: Collective
1411: Input Parameters:
1412: + sw - a `DMSWARM`
1413: - dm - the `DM` to attach to the `DMSWARM`
1415: Level: beginner
1417: Note:
1418: The attached `DM` (dm) will be queried for point location and
1419: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1421: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
1422: @*/
1423: PetscErrorCode DMSwarmSetCellDM(DM sw, DM dm)
1424: {
1425: DMSwarmCellDM celldm;
1426: const char *name;
1427: char *coordName;
1429: PetscFunctionBegin;
1432: PetscCall(PetscStrallocpy(DMSwarmPICField_coor, &coordName));
1433: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 1, (const char **)&coordName, &celldm));
1434: PetscCall(PetscFree(coordName));
1435: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1436: PetscCall(DMSwarmAddCellDM(sw, celldm));
1437: PetscCall(DMSwarmCellDMDestroy(&celldm));
1438: PetscCall(DMSwarmSetCellDMActive(sw, name));
1439: PetscFunctionReturn(PETSC_SUCCESS);
1440: }
1442: /*@
1443: DMSwarmGetCellDM - Fetches the active cell `DM`
1445: Collective
1447: Input Parameter:
1448: . sw - a `DMSWARM`
1450: Output Parameter:
1451: . dm - the active `DM` for the `DMSWARM`
1453: Level: beginner
1455: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1456: @*/
1457: PetscErrorCode DMSwarmGetCellDM(DM sw, DM *dm)
1458: {
1459: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1460: DMSwarmCellDM celldm;
1462: PetscFunctionBegin;
1464: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)&celldm));
1465: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "There is no cell DM named %s in this Swarm", swarm->activeCellDM);
1466: PetscCall(DMSwarmCellDMGetDM(celldm, dm));
1467: PetscFunctionReturn(PETSC_SUCCESS);
1468: }
1470: /*@C
1471: DMSwarmGetCellDMNames - Get the list of cell `DM` names
1473: Not collective
1475: Input Parameter:
1476: . sw - a `DMSWARM`
1478: Output Parameters:
1479: + Ndm - the number of `DMSwarmCellDM` in the `DMSWARM`
1480: - celldms - the name of each `DMSwarmCellDM`
1482: Level: beginner
1484: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMByName()`
1485: @*/
1486: PetscErrorCode DMSwarmGetCellDMNames(DM sw, PetscInt *Ndm, const char **celldms[])
1487: {
1488: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1489: PetscObjectList next = swarm->cellDMs;
1490: PetscInt n = 0;
1492: PetscFunctionBegin;
1494: PetscAssertPointer(Ndm, 2);
1495: PetscAssertPointer(celldms, 3);
1496: while (next) {
1497: next = next->next;
1498: ++n;
1499: }
1500: PetscCall(PetscMalloc1(n, celldms));
1501: next = swarm->cellDMs;
1502: n = 0;
1503: while (next) {
1504: (*celldms)[n] = (const char *)next->obj->name;
1505: next = next->next;
1506: ++n;
1507: }
1508: *Ndm = n;
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1512: /*@
1513: DMSwarmSetCellDMActive - Activates a cell `DM` for a `DMSWARM`
1515: Collective
1517: Input Parameters:
1518: + sw - a `DMSWARM`
1519: - name - name of the cell `DM` to active for the `DMSWARM`
1521: Level: beginner
1523: Note:
1524: The attached `DM` (dmcell) will be queried for point location and
1525: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1527: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1528: @*/
1529: PetscErrorCode DMSwarmSetCellDMActive(DM sw, const char name[])
1530: {
1531: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1532: DMSwarmCellDM celldm;
1534: PetscFunctionBegin;
1536: PetscCall(PetscInfo(sw, "Setting cell DM to %s\n", name));
1537: PetscCall(PetscFree(swarm->activeCellDM));
1538: PetscCall(PetscStrallocpy(name, (char **)&swarm->activeCellDM));
1539: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1540: PetscFunctionReturn(PETSC_SUCCESS);
1541: }
1543: /*@
1544: DMSwarmGetCellDMActive - Returns the active cell `DM` for a `DMSWARM`
1546: Collective
1548: Input Parameter:
1549: . sw - a `DMSWARM`
1551: Output Parameter:
1552: . celldm - the active `DMSwarmCellDM`
1554: Level: beginner
1556: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1557: @*/
1558: PetscErrorCode DMSwarmGetCellDMActive(DM sw, DMSwarmCellDM *celldm)
1559: {
1560: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1562: PetscFunctionBegin;
1564: PetscAssertPointer(celldm, 2);
1565: PetscCheck(swarm->activeCellDM, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no active cell DM");
1566: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)celldm));
1567: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", swarm->activeCellDM);
1568: PetscFunctionReturn(PETSC_SUCCESS);
1569: }
1571: /*@C
1572: DMSwarmGetCellDMByName - Get a `DMSwarmCellDM` from its name
1574: Not collective
1576: Input Parameters:
1577: + sw - a `DMSWARM`
1578: - name - the name
1580: Output Parameter:
1581: . celldm - the `DMSwarmCellDM`
1583: Level: beginner
1585: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMNames()`
1586: @*/
1587: PetscErrorCode DMSwarmGetCellDMByName(DM sw, const char name[], DMSwarmCellDM *celldm)
1588: {
1589: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1591: PetscFunctionBegin;
1593: PetscAssertPointer(name, 2);
1594: PetscAssertPointer(celldm, 3);
1595: PetscCall(PetscObjectListFind(swarm->cellDMs, name, (PetscObject *)celldm));
1596: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", name);
1597: PetscFunctionReturn(PETSC_SUCCESS);
1598: }
1600: /*@
1601: DMSwarmAddCellDM - Adds a cell `DM` to the `DMSWARM`
1603: Collective
1605: Input Parameters:
1606: + sw - a `DMSWARM`
1607: - celldm - the `DMSwarmCellDM`
1609: Level: beginner
1611: Note:
1612: Cell DMs with the same name will share the cellid field
1614: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPushCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1615: @*/
1616: PetscErrorCode DMSwarmAddCellDM(DM sw, DMSwarmCellDM celldm)
1617: {
1618: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1619: const char *name;
1620: PetscInt dim;
1621: PetscBool flg;
1622: MPI_Comm comm;
1624: PetscFunctionBegin;
1626: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1628: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1629: PetscCall(PetscObjectListAdd(&swarm->cellDMs, name, (PetscObject)celldm));
1630: PetscCall(DMGetDimension(sw, &dim));
1631: for (PetscInt f = 0; f < celldm->Nfc; ++f) {
1632: PetscCall(DMSwarmDataFieldStringInList(celldm->coordFields[f], swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1633: if (!flg) {
1634: PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->coordFields[f], dim, PETSC_DOUBLE));
1635: } else {
1636: PetscDataType dt;
1637: PetscInt bs;
1639: PetscCall(DMSwarmGetFieldInfo(sw, celldm->coordFields[f], &bs, &dt));
1640: PetscCheck(bs == dim, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has blocksize %" PetscInt_FMT " != %" PetscInt_FMT " spatial dimension", celldm->coordFields[f], bs, dim);
1641: PetscCheck(dt == PETSC_DOUBLE, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has datatype %s != PETSC_DOUBLE", celldm->coordFields[f], PetscDataTypes[dt]);
1642: }
1643: }
1644: // Assume that DMs with the same name share the cellid field
1645: PetscCall(DMSwarmDataFieldStringInList(celldm->cellid, swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1646: if (!flg) {
1647: PetscBool isShell, isDummy;
1648: const char *name;
1650: // Allow dummy DMSHELL (I don't think we should support this mode)
1651: PetscCall(PetscObjectTypeCompare((PetscObject)celldm->dm, DMSHELL, &isShell));
1652: PetscCall(PetscObjectGetName((PetscObject)celldm->dm, &name));
1653: PetscCall(PetscStrcmp(name, "dummy", &isDummy));
1654: if (!isShell || !isDummy) PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->cellid, 1, PETSC_INT));
1655: }
1656: PetscCall(DMSwarmSetCellDMActive(sw, name));
1657: PetscFunctionReturn(PETSC_SUCCESS);
1658: }
1660: /*@
1661: DMSwarmGetLocalSize - Retrieves the local length of fields registered
1663: Not Collective
1665: Input Parameter:
1666: . dm - a `DMSWARM`
1668: Output Parameter:
1669: . nlocal - the length of each registered field
1671: Level: beginner
1673: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
1674: @*/
1675: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
1676: {
1677: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1679: PetscFunctionBegin;
1680: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: /*@
1685: DMSwarmGetSize - Retrieves the total length of fields registered
1687: Collective
1689: Input Parameter:
1690: . dm - a `DMSWARM`
1692: Output Parameter:
1693: . n - the total length of each registered field
1695: Level: beginner
1697: Note:
1698: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
1700: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
1701: @*/
1702: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
1703: {
1704: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1705: PetscInt nlocal;
1707: PetscFunctionBegin;
1708: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1709: PetscCallMPI(MPIU_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1710: PetscFunctionReturn(PETSC_SUCCESS);
1711: }
1713: /*@C
1714: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
1716: Collective
1718: Input Parameters:
1719: + dm - a `DMSWARM`
1720: . fieldname - the textual name to identify this field
1721: . blocksize - the number of each data type
1722: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
1724: Level: beginner
1726: Notes:
1727: The textual name for each registered field must be unique.
1729: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1730: @*/
1731: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1732: {
1733: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1734: size_t size;
1736: PetscFunctionBegin;
1737: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1738: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1740: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1741: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1742: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1743: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1744: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1746: PetscCall(PetscDataTypeGetSize(type, &size));
1747: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1748: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1749: {
1750: DMSwarmDataField gfield;
1752: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1753: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1754: }
1755: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1756: PetscFunctionReturn(PETSC_SUCCESS);
1757: }
1759: /*@C
1760: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
1762: Collective
1764: Input Parameters:
1765: + dm - a `DMSWARM`
1766: . fieldname - the textual name to identify this field
1767: - size - the size in bytes of the user struct of each data type
1769: Level: beginner
1771: Note:
1772: The textual name for each registered field must be unique.
1774: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1775: @*/
1776: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1777: {
1778: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1780: PetscFunctionBegin;
1781: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1782: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1783: PetscFunctionReturn(PETSC_SUCCESS);
1784: }
1786: /*@C
1787: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1789: Collective
1791: Input Parameters:
1792: + dm - a `DMSWARM`
1793: . fieldname - the textual name to identify this field
1794: . size - the size in bytes of the user data type
1795: - blocksize - the number of each data type
1797: Level: beginner
1799: Note:
1800: The textual name for each registered field must be unique.
1802: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1803: @*/
1804: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1805: {
1806: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1808: PetscFunctionBegin;
1809: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1810: {
1811: DMSwarmDataField gfield;
1813: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1814: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1815: }
1816: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1817: PetscFunctionReturn(PETSC_SUCCESS);
1818: }
1820: /*@C
1821: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1823: Not Collective, No Fortran Support
1825: Input Parameters:
1826: + dm - a `DMSWARM`
1827: - fieldname - the textual name to identify this field
1829: Output Parameters:
1830: + blocksize - the number of each data type
1831: . type - the data type
1832: - data - pointer to raw array
1834: Level: beginner
1836: Notes:
1837: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1839: Fortran Note:
1840: Only works for `type` of `PETSC_SCALAR`
1842: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1843: @*/
1844: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1845: {
1846: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1847: DMSwarmDataField gfield;
1849: PetscFunctionBegin;
1851: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1852: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1853: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1854: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1855: if (blocksize) *blocksize = gfield->bs;
1856: if (type) *type = gfield->petsc_type;
1857: PetscFunctionReturn(PETSC_SUCCESS);
1858: }
1860: /*@C
1861: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1863: Not Collective
1865: Input Parameters:
1866: + dm - a `DMSWARM`
1867: - fieldname - the textual name to identify this field
1869: Output Parameters:
1870: + blocksize - the number of each data type
1871: . type - the data type
1872: - data - pointer to raw array
1874: Level: beginner
1876: Notes:
1877: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1879: Fortran Note:
1880: Only works for `type` of `PETSC_SCALAR`
1882: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1883: @*/
1884: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1885: {
1886: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1887: DMSwarmDataField gfield;
1889: PetscFunctionBegin;
1891: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1892: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1893: if (data) *data = NULL;
1894: PetscFunctionReturn(PETSC_SUCCESS);
1895: }
1897: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
1898: {
1899: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1900: DMSwarmDataField gfield;
1902: PetscFunctionBegin;
1904: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1905: if (blocksize) *blocksize = gfield->bs;
1906: if (type) *type = gfield->petsc_type;
1907: PetscFunctionReturn(PETSC_SUCCESS);
1908: }
1910: /*@
1911: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
1913: Not Collective
1915: Input Parameter:
1916: . dm - a `DMSWARM`
1918: Level: beginner
1920: Notes:
1921: The new point will have all fields initialized to zero.
1923: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
1924: @*/
1925: PetscErrorCode DMSwarmAddPoint(DM dm)
1926: {
1927: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1929: PetscFunctionBegin;
1930: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1931: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1932: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
1933: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1934: PetscFunctionReturn(PETSC_SUCCESS);
1935: }
1937: /*@
1938: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
1940: Not Collective
1942: Input Parameters:
1943: + dm - a `DMSWARM`
1944: - npoints - the number of new points to add
1946: Level: beginner
1948: Notes:
1949: The new point will have all fields initialized to zero.
1951: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
1952: @*/
1953: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1954: {
1955: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1956: PetscInt nlocal;
1958: PetscFunctionBegin;
1959: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1960: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1961: nlocal = PetscMax(nlocal, 0) + npoints;
1962: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
1963: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1964: PetscFunctionReturn(PETSC_SUCCESS);
1965: }
1967: /*@
1968: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
1970: Not Collective
1972: Input Parameter:
1973: . dm - a `DMSWARM`
1975: Level: beginner
1977: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
1978: @*/
1979: PetscErrorCode DMSwarmRemovePoint(DM dm)
1980: {
1981: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1983: PetscFunctionBegin;
1984: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1985: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
1986: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1987: PetscFunctionReturn(PETSC_SUCCESS);
1988: }
1990: /*@
1991: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
1993: Not Collective
1995: Input Parameters:
1996: + dm - a `DMSWARM`
1997: - idx - index of point to remove
1999: Level: beginner
2001: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2002: @*/
2003: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
2004: {
2005: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2007: PetscFunctionBegin;
2008: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
2009: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
2010: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
2011: PetscFunctionReturn(PETSC_SUCCESS);
2012: }
2014: /*@
2015: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
2017: Not Collective
2019: Input Parameters:
2020: + dm - a `DMSWARM`
2021: . pi - the index of the point to copy
2022: - pj - the point index where the copy should be located
2024: Level: beginner
2026: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2027: @*/
2028: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
2029: {
2030: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2032: PetscFunctionBegin;
2033: if (!swarm->issetup) PetscCall(DMSetUp(dm));
2034: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
2035: PetscFunctionReturn(PETSC_SUCCESS);
2036: }
2038: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
2039: {
2040: PetscFunctionBegin;
2041: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
2042: PetscFunctionReturn(PETSC_SUCCESS);
2043: }
2045: /*@
2046: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
2048: Collective
2050: Input Parameters:
2051: + dm - the `DMSWARM`
2052: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
2054: Level: advanced
2056: Notes:
2057: The `DM` will be modified to accommodate received points.
2058: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
2059: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
2061: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
2062: @*/
2063: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
2064: {
2065: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2067: PetscFunctionBegin;
2068: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
2069: switch (swarm->migrate_type) {
2070: case DMSWARM_MIGRATE_BASIC:
2071: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
2072: break;
2073: case DMSWARM_MIGRATE_DMCELLNSCATTER:
2074: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
2075: break;
2076: case DMSWARM_MIGRATE_DMCELLEXACT:
2077: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
2078: case DMSWARM_MIGRATE_USER:
2079: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
2080: default:
2081: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
2082: }
2083: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
2084: PetscCall(DMClearGlobalVectors(dm));
2085: PetscFunctionReturn(PETSC_SUCCESS);
2086: }
2088: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
2090: /*
2091: DMSwarmCollectViewCreate
2093: * Applies a collection method and gathers point neighbour points into dm
2095: Notes:
2096: Users should call DMSwarmCollectViewDestroy() after
2097: they have finished computations associated with the collected points
2098: */
2100: /*@
2101: DMSwarmCollectViewCreate - Applies a collection method and gathers points
2102: in neighbour ranks into the `DMSWARM`
2104: Collective
2106: Input Parameter:
2107: . dm - the `DMSWARM`
2109: Level: advanced
2111: Notes:
2112: Users should call `DMSwarmCollectViewDestroy()` after
2113: they have finished computations associated with the collected points
2115: Different collect methods are supported. See `DMSwarmSetCollectType()`.
2117: Developer Note:
2118: Create and Destroy routines create new objects that can get destroyed, they do not change the state
2119: of the current object.
2121: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
2122: @*/
2123: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
2124: {
2125: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2126: PetscInt ng;
2128: PetscFunctionBegin;
2129: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
2130: PetscCall(DMSwarmGetLocalSize(dm, &ng));
2131: switch (swarm->collect_type) {
2132: case DMSWARM_COLLECT_BASIC:
2133: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
2134: break;
2135: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
2136: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
2137: case DMSWARM_COLLECT_GENERAL:
2138: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
2139: default:
2140: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
2141: }
2142: swarm->collect_view_active = PETSC_TRUE;
2143: swarm->collect_view_reset_nlocal = ng;
2144: PetscFunctionReturn(PETSC_SUCCESS);
2145: }
2147: /*@
2148: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
2150: Collective
2152: Input Parameters:
2153: . dm - the `DMSWARM`
2155: Notes:
2156: Users should call `DMSwarmCollectViewCreate()` before this function is called.
2158: Level: advanced
2160: Developer Note:
2161: Create and Destroy routines create new objects that can get destroyed, they do not change the state
2162: of the current object.
2164: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
2165: @*/
2166: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
2167: {
2168: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2170: PetscFunctionBegin;
2171: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
2172: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
2173: swarm->collect_view_active = PETSC_FALSE;
2174: PetscFunctionReturn(PETSC_SUCCESS);
2175: }
2177: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
2178: {
2179: PetscInt dim;
2181: PetscFunctionBegin;
2182: PetscCall(DMSwarmSetNumSpecies(dm, 1));
2183: PetscCall(DMGetDimension(dm, &dim));
2184: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2185: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2186: PetscFunctionReturn(PETSC_SUCCESS);
2187: }
2189: /*@
2190: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
2192: Collective
2194: Input Parameters:
2195: + dm - the `DMSWARM`
2196: - Npc - The number of particles per cell in the cell `DM`
2198: Level: intermediate
2200: Notes:
2201: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
2202: one particle is in each cell, it is placed at the centroid.
2204: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
2205: @*/
2206: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
2207: {
2208: DM cdm;
2209: DMSwarmCellDM celldm;
2210: PetscRandom rnd;
2211: DMPolytopeType ct;
2212: PetscBool simplex;
2213: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
2214: PetscInt dim, d, cStart, cEnd, c, p, Nfc;
2215: const char **coordFields;
2217: PetscFunctionBeginUser;
2218: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
2219: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
2220: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
2222: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2223: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2224: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2225: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2226: PetscCall(DMGetDimension(cdm, &dim));
2227: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
2228: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
2229: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
2231: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
2232: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
2233: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2234: for (c = cStart; c < cEnd; ++c) {
2235: if (Npc == 1) {
2236: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
2237: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
2238: } else {
2239: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
2240: for (p = 0; p < Npc; ++p) {
2241: const PetscInt n = c * Npc + p;
2242: PetscReal sum = 0.0, refcoords[3];
2244: for (d = 0; d < dim; ++d) {
2245: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
2246: sum += refcoords[d];
2247: }
2248: if (simplex && sum > 0.0)
2249: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
2250: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
2251: }
2252: }
2253: }
2254: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2255: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
2256: PetscCall(PetscRandomDestroy(&rnd));
2257: PetscFunctionReturn(PETSC_SUCCESS);
2258: }
2260: /*@
2261: DMSwarmGetType - Get particular flavor of `DMSWARM`
2263: Collective
2265: Input Parameter:
2266: . sw - the `DMSWARM`
2268: Output Parameter:
2269: . stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2271: Level: advanced
2273: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2274: @*/
2275: PetscErrorCode DMSwarmGetType(DM sw, DMSwarmType *stype)
2276: {
2277: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2279: PetscFunctionBegin;
2281: PetscAssertPointer(stype, 2);
2282: *stype = swarm->swarm_type;
2283: PetscFunctionReturn(PETSC_SUCCESS);
2284: }
2286: /*@
2287: DMSwarmSetType - Set particular flavor of `DMSWARM`
2289: Collective
2291: Input Parameters:
2292: + sw - the `DMSWARM`
2293: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2295: Level: advanced
2297: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2298: @*/
2299: PetscErrorCode DMSwarmSetType(DM sw, DMSwarmType stype)
2300: {
2301: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2303: PetscFunctionBegin;
2305: swarm->swarm_type = stype;
2306: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(sw));
2307: PetscFunctionReturn(PETSC_SUCCESS);
2308: }
2310: static PetscErrorCode DMSwarmCreateRemapDM_Private(DM sw, DM *rdm)
2311: {
2312: PetscFE fe;
2313: DMPolytopeType ct;
2314: PetscInt dim, cStart;
2315: const char *prefix = "remap_";
2317: PetscFunctionBegin;
2318: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), rdm));
2319: PetscCall(DMSetType(*rdm, DMPLEX));
2320: PetscCall(DMPlexSetOptionsPrefix(*rdm, prefix));
2321: PetscCall(DMSetFromOptions(*rdm));
2322: PetscCall(PetscObjectSetName((PetscObject)*rdm, "remap"));
2323: PetscCall(DMViewFromOptions(*rdm, NULL, "-dm_view"));
2325: PetscCall(DMGetDimension(*rdm, &dim));
2326: PetscCall(DMPlexGetHeightStratum(*rdm, 0, &cStart, NULL));
2327: PetscCall(DMPlexGetCellType(*rdm, cStart, &ct));
2328: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
2329: PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
2330: PetscCall(DMSetField(*rdm, 0, NULL, (PetscObject)fe));
2331: PetscCall(DMCreateDS(*rdm));
2332: PetscCall(PetscFEDestroy(&fe));
2333: PetscFunctionReturn(PETSC_SUCCESS);
2334: }
2336: static PetscErrorCode DMSetup_Swarm(DM sw)
2337: {
2338: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2340: PetscFunctionBegin;
2341: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
2342: swarm->issetup = PETSC_TRUE;
2344: if (swarm->remap_type != DMSWARM_REMAP_NONE) {
2345: DMSwarmCellDM celldm;
2346: DM rdm;
2347: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
2348: const char *vfieldnames[1] = {"w_q"};
2350: PetscCall(DMSwarmCreateRemapDM_Private(sw, &rdm));
2351: PetscCall(DMSwarmCellDMCreate(rdm, 1, vfieldnames, 2, fieldnames, &celldm));
2352: PetscCall(DMSwarmAddCellDM(sw, celldm));
2353: PetscCall(DMSwarmCellDMDestroy(&celldm));
2354: PetscCall(DMDestroy(&rdm));
2355: }
2357: if (swarm->swarm_type == DMSWARM_PIC) {
2358: DMSwarmCellDM celldm;
2360: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2361: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "No active cell DM. DMSWARM_PIC requires you call DMSwarmSetCellDM() or DMSwarmAddCellDM()");
2362: if (celldm->dm->ops->locatepointssubdomain) {
2363: /* check methods exists for exact ownership identificiation */
2364: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
2365: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
2366: } else {
2367: /* check methods exist for point location AND rank neighbor identification */
2368: PetscCheck(celldm->dm->ops->locatepoints, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
2369: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
2371: PetscCheck(celldm->dm->ops->getneighbors, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
2372: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
2374: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
2375: }
2376: }
2378: PetscCall(DMSwarmFinalizeFieldRegister(sw));
2380: /* check some fields were registered */
2381: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
2382: PetscFunctionReturn(PETSC_SUCCESS);
2383: }
2385: static PetscErrorCode DMDestroy_Swarm(DM dm)
2386: {
2387: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2389: PetscFunctionBegin;
2390: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2391: PetscCall(PetscObjectListDestroy(&swarm->cellDMs));
2392: PetscCall(PetscFree(swarm->activeCellDM));
2393: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
2394: PetscCall(PetscFree(swarm));
2395: PetscFunctionReturn(PETSC_SUCCESS);
2396: }
2398: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
2399: {
2400: DM cdm;
2401: DMSwarmCellDM celldm;
2402: PetscDraw draw;
2403: PetscReal *coords, oldPause, radius = 0.01;
2404: PetscInt Np, p, bs, Nfc;
2405: const char **coordFields;
2407: PetscFunctionBegin;
2408: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
2409: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2410: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2411: PetscCall(PetscDrawGetPause(draw, &oldPause));
2412: PetscCall(PetscDrawSetPause(draw, 0.0));
2413: PetscCall(DMView(cdm, viewer));
2414: PetscCall(PetscDrawSetPause(draw, oldPause));
2416: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2417: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2418: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2419: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2420: PetscCall(DMSwarmGetField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2421: for (p = 0; p < Np; ++p) {
2422: const PetscInt i = p * bs;
2424: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
2425: }
2426: PetscCall(DMSwarmRestoreField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2427: PetscCall(PetscDrawFlush(draw));
2428: PetscCall(PetscDrawPause(draw));
2429: PetscCall(PetscDrawSave(draw));
2430: PetscFunctionReturn(PETSC_SUCCESS);
2431: }
2433: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
2434: {
2435: PetscViewerFormat format;
2436: PetscInt *sizes;
2437: PetscInt dim, Np, maxSize = 17;
2438: MPI_Comm comm;
2439: PetscMPIInt rank, size, p;
2440: const char *name, *cellid;
2442: PetscFunctionBegin;
2443: PetscCall(PetscViewerGetFormat(viewer, &format));
2444: PetscCall(DMGetDimension(dm, &dim));
2445: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2446: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2447: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2448: PetscCallMPI(MPI_Comm_size(comm, &size));
2449: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
2450: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
2451: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
2452: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
2453: else PetscCall(PetscCalloc1(3, &sizes));
2454: if (size < maxSize) {
2455: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
2456: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
2457: for (p = 0; p < size; ++p) {
2458: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
2459: }
2460: } else {
2461: PetscInt locMinMax[2] = {Np, Np};
2463: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
2464: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
2465: }
2466: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2467: PetscCall(PetscFree(sizes));
2468: if (format == PETSC_VIEWER_ASCII_INFO) {
2469: DMSwarmCellDM celldm;
2470: PetscInt *cell;
2472: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
2473: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2474: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2475: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
2476: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&cell));
2477: for (p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%d: %" PetscInt_FMT "\n", p, cell[p]));
2478: PetscCall(PetscViewerFlush(viewer));
2479: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2480: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&cell));
2481: }
2482: PetscFunctionReturn(PETSC_SUCCESS);
2483: }
2485: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
2486: {
2487: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2488: PetscBool isascii, ibinary, isvtk, isdraw, ispython;
2489: #if defined(PETSC_HAVE_HDF5)
2490: PetscBool ishdf5;
2491: #endif
2493: PetscFunctionBegin;
2496: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2497: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
2498: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
2499: #if defined(PETSC_HAVE_HDF5)
2500: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2501: #endif
2502: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2503: PetscCall(PetscObjectHasFunction((PetscObject)viewer, "PetscViewerPythonViewObject_C", &ispython));
2504: if (isascii) {
2505: PetscViewerFormat format;
2507: PetscCall(PetscViewerGetFormat(viewer, &format));
2508: switch (format) {
2509: case PETSC_VIEWER_ASCII_INFO_DETAIL:
2510: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
2511: break;
2512: default:
2513: PetscCall(DMView_Swarm_Ascii(dm, viewer));
2514: }
2515: } else {
2516: #if defined(PETSC_HAVE_HDF5)
2517: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
2518: #endif
2519: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
2520: if (ispython) PetscCall(PetscViewerPythonViewObject(viewer, (PetscObject)dm));
2521: }
2522: PetscFunctionReturn(PETSC_SUCCESS);
2523: }
2525: /*@
2526: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
2527: 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.
2529: Noncollective
2531: Input Parameters:
2532: + sw - the `DMSWARM`
2533: . cellID - the integer id of the cell to be extracted and filtered
2534: - cellswarm - The `DMSWARM` to receive the cell
2536: Level: beginner
2538: Notes:
2539: This presently only supports `DMSWARM_PIC` type
2541: Should be restored with `DMSwarmRestoreCellSwarm()`
2543: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
2545: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
2546: @*/
2547: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2548: {
2549: DM_Swarm *original = (DM_Swarm *)sw->data;
2550: DMLabel label;
2551: DM dmc, subdmc;
2552: PetscInt *pids, particles, dim;
2553: const char *name;
2555: PetscFunctionBegin;
2556: /* Configure new swarm */
2557: PetscCall(DMSetType(cellswarm, DMSWARM));
2558: PetscCall(DMGetDimension(sw, &dim));
2559: PetscCall(DMSetDimension(cellswarm, dim));
2560: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
2561: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
2562: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
2563: PetscCall(DMSwarmSortGetAccess(sw));
2564: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
2565: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2566: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
2567: PetscCall(DMSwarmSortRestoreAccess(sw));
2568: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2569: PetscCall(DMSwarmGetCellDM(sw, &dmc));
2570: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
2571: PetscCall(DMAddLabel(dmc, label));
2572: PetscCall(DMLabelSetValue(label, cellID, 1));
2573: PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdmc));
2574: PetscCall(PetscObjectGetName((PetscObject)dmc, &name));
2575: PetscCall(PetscObjectSetName((PetscObject)subdmc, name));
2576: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
2577: PetscCall(DMLabelDestroy(&label));
2578: PetscFunctionReturn(PETSC_SUCCESS);
2579: }
2581: /*@
2582: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
2584: Noncollective
2586: Input Parameters:
2587: + sw - the parent `DMSWARM`
2588: . cellID - the integer id of the cell to be copied back into the parent swarm
2589: - cellswarm - the cell swarm object
2591: Level: beginner
2593: Note:
2594: This only supports `DMSWARM_PIC` types of `DMSWARM`s
2596: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
2597: @*/
2598: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2599: {
2600: DM dmc;
2601: PetscInt *pids, particles, p;
2603: PetscFunctionBegin;
2604: PetscCall(DMSwarmSortGetAccess(sw));
2605: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2606: PetscCall(DMSwarmSortRestoreAccess(sw));
2607: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
2608: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
2609: /* Free memory, destroy cell dm */
2610: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
2611: PetscCall(DMDestroy(&dmc));
2612: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2613: PetscFunctionReturn(PETSC_SUCCESS);
2614: }
2616: /*@
2617: DMSwarmComputeMoments - Compute the first three particle moments for a given field
2619: Noncollective
2621: Input Parameters:
2622: + sw - the `DMSWARM`
2623: . coordinate - the coordinate field name
2624: - weight - the weight field name
2626: Output Parameter:
2627: . moments - the field moments
2629: Level: intermediate
2631: Notes:
2632: The `moments` array should be of length bs + 2, where bs is the block size of the coordinate field.
2634: The weight field must be a scalar, having blocksize 1.
2636: .seealso: `DM`, `DMSWARM`, `DMPlexComputeMoments()`
2637: @*/
2638: PetscErrorCode DMSwarmComputeMoments(DM sw, const char coordinate[], const char weight[], PetscReal moments[])
2639: {
2640: const PetscReal *coords;
2641: const PetscReal *w;
2642: PetscReal *mom;
2643: PetscDataType dtc, dtw;
2644: PetscInt bsc, bsw, Np;
2645: MPI_Comm comm;
2647: PetscFunctionBegin;
2649: PetscAssertPointer(coordinate, 2);
2650: PetscAssertPointer(weight, 3);
2651: PetscAssertPointer(moments, 4);
2652: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
2653: PetscCall(DMSwarmGetField(sw, coordinate, &bsc, &dtc, (void **)&coords));
2654: PetscCall(DMSwarmGetField(sw, weight, &bsw, &dtw, (void **)&w));
2655: PetscCheck(dtc == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s must be real, not %s", coordinate, PetscDataTypes[dtc]);
2656: PetscCheck(dtw == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be real, not %s", weight, PetscDataTypes[dtw]);
2657: PetscCheck(bsw == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be a scalar, not blocksize %" PetscInt_FMT, weight, bsw);
2658: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2659: PetscCall(DMGetWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2660: PetscCall(PetscArrayzero(mom, bsc + 2));
2661: for (PetscInt p = 0; p < Np; ++p) {
2662: const PetscReal *c = &coords[p * bsc];
2663: const PetscReal wp = w[p];
2665: mom[0] += wp;
2666: for (PetscInt d = 0; d < bsc; ++d) {
2667: mom[d + 1] += wp * c[d];
2668: mom[d + bsc + 1] += wp * PetscSqr(c[d]);
2669: }
2670: }
2671: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&coords));
2672: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
2673: PetscCallMPI(MPIU_Allreduce(mom, moments, bsc + 2, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
2674: PetscCall(DMRestoreWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2675: PetscFunctionReturn(PETSC_SUCCESS);
2676: }
2678: static PetscErrorCode DMSetFromOptions_Swarm(DM dm, PetscOptionItems PetscOptionsObject)
2679: {
2680: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2682: PetscFunctionBegin;
2683: PetscOptionsHeadBegin(PetscOptionsObject, "DMSwarm Options");
2684: PetscCall(PetscOptionsEnum("-dm_swarm_remap_type", "Remap algorithm", "DMSwarmSetRemapType", DMSwarmRemapTypeNames, (PetscEnum)swarm->remap_type, (PetscEnum *)&swarm->remap_type, NULL));
2685: PetscOptionsHeadEnd();
2686: PetscFunctionReturn(PETSC_SUCCESS);
2687: }
2689: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
2691: static PetscErrorCode DMInitialize_Swarm(DM sw)
2692: {
2693: PetscFunctionBegin;
2694: sw->ops->view = DMView_Swarm;
2695: sw->ops->load = NULL;
2696: sw->ops->setfromoptions = DMSetFromOptions_Swarm;
2697: sw->ops->clone = DMClone_Swarm;
2698: sw->ops->setup = DMSetup_Swarm;
2699: sw->ops->createlocalsection = NULL;
2700: sw->ops->createsectionpermutation = NULL;
2701: sw->ops->createdefaultconstraints = NULL;
2702: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
2703: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
2704: sw->ops->getlocaltoglobalmapping = NULL;
2705: sw->ops->createfieldis = NULL;
2706: sw->ops->createcoordinatedm = NULL;
2707: sw->ops->createcellcoordinatedm = NULL;
2708: sw->ops->getcoloring = NULL;
2709: sw->ops->creatematrix = DMCreateMatrix_Swarm;
2710: sw->ops->createinterpolation = NULL;
2711: sw->ops->createinjection = NULL;
2712: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
2713: sw->ops->creategradientmatrix = DMCreateGradientMatrix_Swarm;
2714: sw->ops->refine = NULL;
2715: sw->ops->coarsen = NULL;
2716: sw->ops->refinehierarchy = NULL;
2717: sw->ops->coarsenhierarchy = NULL;
2718: sw->ops->globaltolocalbegin = DMGlobalToLocalBegin_Swarm;
2719: sw->ops->globaltolocalend = DMGlobalToLocalEnd_Swarm;
2720: sw->ops->localtoglobalbegin = DMLocalToGlobalBegin_Swarm;
2721: sw->ops->localtoglobalend = DMLocalToGlobalEnd_Swarm;
2722: sw->ops->destroy = DMDestroy_Swarm;
2723: sw->ops->createsubdm = NULL;
2724: sw->ops->getdimpoints = NULL;
2725: sw->ops->locatepoints = NULL;
2726: sw->ops->projectfieldlocal = DMProjectFieldLocal_Swarm;
2727: PetscFunctionReturn(PETSC_SUCCESS);
2728: }
2730: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
2731: {
2732: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2734: PetscFunctionBegin;
2735: swarm->refct++;
2736: (*newdm)->data = swarm;
2737: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
2738: PetscCall(DMInitialize_Swarm(*newdm));
2739: (*newdm)->dim = dm->dim;
2740: PetscFunctionReturn(PETSC_SUCCESS);
2741: }
2743: /*MC
2744: DMSWARM = "swarm" - A `DM` object for particle methods, such as particle-in-cell (PIC), in which the underlying
2745: data is both (i) dynamic in length, (ii) and of arbitrary data type.
2747: Level: intermediate
2749: Notes:
2750: User data can be represented by `DMSWARM` through a registering "fields" which are to be stored on particles.
2751: To register a field, the user must provide;
2752: (a) a unique name;
2753: (b) the data type (or size in bytes);
2754: (c) the block size of the data.
2756: For example, suppose the application requires a unique id, energy, momentum and density to be stored
2757: on a set of particles. Then the following code could be used
2758: .vb
2759: DMSwarmInitializeFieldRegister(dm)
2760: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
2761: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
2762: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
2763: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
2764: DMSwarmFinalizeFieldRegister(dm)
2765: .ve
2767: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
2768: The only restriction imposed by `DMSWARM` is that all fields contain the same number of particles.
2770: To support particle methods, "migration" techniques are provided. These methods migrate data
2771: between ranks.
2773: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
2774: As a `DMSWARM` may internally define and store values of different data types,
2775: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
2776: fields should be used to define a `Vec` object via `DMSwarmVectorDefineField()`
2777: The specified field can be changed at any time - thereby permitting vectors
2778: compatible with different fields to be created.
2780: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via `DMSwarmCreateGlobalVectorFromField()`
2781: Here the data defining the field in the `DMSWARM` is shared with a `Vec`.
2782: This is inherently unsafe if you alter the size of the field at any time between
2783: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
2784: If the local size of the `DMSWARM` does not match the local size of the global vector
2785: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
2787: Additional high-level support is provided for Particle-In-Cell methods. Refer to `DMSwarmSetType()`.
2789: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`, `DMSwarmSetType()`, `DMSwarmType`, `DMSwarmCreateGlobalVectorFromField()`,
2790: `DMCreateGlobalVector()`, `DMCreateLocalVector()`
2791: M*/
2793: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
2794: {
2795: DM_Swarm *swarm;
2797: PetscFunctionBegin;
2799: PetscCall(PetscNew(&swarm));
2800: dm->data = swarm;
2801: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
2802: PetscCall(DMSwarmInitializeFieldRegister(dm));
2803: dm->dim = 0;
2804: swarm->refct = 1;
2805: swarm->issetup = PETSC_FALSE;
2806: swarm->swarm_type = DMSWARM_BASIC;
2807: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
2808: swarm->collect_type = DMSWARM_COLLECT_BASIC;
2809: swarm->migrate_error_on_missing_point = PETSC_FALSE;
2810: swarm->collect_view_active = PETSC_FALSE;
2811: swarm->collect_view_reset_nlocal = -1;
2812: PetscCall(DMInitialize_Swarm(dm));
2813: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
2814: PetscFunctionReturn(PETSC_SUCCESS);
2815: }
2817: /* Replace dm with the contents of ndm, and then destroy ndm
2818: - Share the DM_Swarm structure
2819: */
2820: PetscErrorCode DMSwarmReplace(DM dm, DM *ndm)
2821: {
2822: DM dmNew = *ndm;
2823: const PetscReal *maxCell, *Lstart, *L;
2824: PetscInt dim;
2826: PetscFunctionBegin;
2827: if (dm == dmNew) {
2828: PetscCall(DMDestroy(ndm));
2829: PetscFunctionReturn(PETSC_SUCCESS);
2830: }
2831: dm->setupcalled = dmNew->setupcalled;
2832: if (!dm->hdr.name) {
2833: const char *name;
2835: PetscCall(PetscObjectGetName((PetscObject)*ndm, &name));
2836: PetscCall(PetscObjectSetName((PetscObject)dm, name));
2837: }
2838: PetscCall(DMGetDimension(dmNew, &dim));
2839: PetscCall(DMSetDimension(dm, dim));
2840: PetscCall(DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L));
2841: PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
2842: PetscCall(DMDestroy_Swarm(dm));
2843: PetscCall(DMInitialize_Swarm(dm));
2844: dm->data = dmNew->data;
2845: ((DM_Swarm *)dmNew->data)->refct++;
2846: PetscCall(DMDestroy(ndm));
2847: PetscFunctionReturn(PETSC_SUCCESS);
2848: }
2850: /*@
2851: DMSwarmDuplicate - Creates a new `DMSWARM` with the same fields and cell `DM`s but no particles
2853: Collective
2855: Input Parameter:
2856: . sw - the `DMSWARM`
2858: Output Parameter:
2859: . nsw - the new `DMSWARM`
2861: Level: beginner
2863: .seealso: `DM`, `DMSWARM`, `DMSwarmCreate()`, `DMClone()`
2864: @*/
2865: PetscErrorCode DMSwarmDuplicate(DM sw, DM *nsw)
2866: {
2867: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2868: DMSwarmDataField *fields;
2869: DMSwarmCellDM celldm, ncelldm;
2870: DMSwarmType stype;
2871: const char *name, **celldmnames;
2872: void *ctx;
2873: PetscInt dim, Nf, Ndm;
2874: PetscBool flg;
2876: PetscFunctionBegin;
2877: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), nsw));
2878: PetscCall(DMSetType(*nsw, DMSWARM));
2879: PetscCall(PetscObjectGetName((PetscObject)sw, &name));
2880: PetscCall(PetscObjectSetName((PetscObject)*nsw, name));
2881: PetscCall(DMGetDimension(sw, &dim));
2882: PetscCall(DMSetDimension(*nsw, dim));
2883: PetscCall(DMSwarmGetType(sw, &stype));
2884: PetscCall(DMSwarmSetType(*nsw, stype));
2885: PetscCall(DMGetApplicationContext(sw, &ctx));
2886: PetscCall(DMSetApplicationContext(*nsw, ctx));
2888: PetscCall(DMSwarmDataBucketGetDMSwarmDataFields(swarm->db, &Nf, &fields));
2889: for (PetscInt f = 0; f < Nf; ++f) {
2890: PetscCall(DMSwarmDataFieldStringInList(fields[f]->name, ((DM_Swarm *)(*nsw)->data)->db->nfields, (const DMSwarmDataField *)((DM_Swarm *)(*nsw)->data)->db->field, &flg));
2891: if (!flg) PetscCall(DMSwarmRegisterPetscDatatypeField(*nsw, fields[f]->name, fields[f]->bs, fields[f]->petsc_type));
2892: }
2894: PetscCall(DMSwarmGetCellDMNames(sw, &Ndm, &celldmnames));
2895: for (PetscInt c = 0; c < Ndm; ++c) {
2896: DM dm;
2897: PetscInt Ncf;
2898: const char **coordfields, **fields;
2900: PetscCall(DMSwarmGetCellDMByName(sw, celldmnames[c], &celldm));
2901: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
2902: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncf, &coordfields));
2903: PetscCall(DMSwarmCellDMGetFields(celldm, &Nf, &fields));
2904: PetscCall(DMSwarmCellDMCreate(dm, Nf, fields, Ncf, coordfields, &ncelldm));
2905: PetscCall(DMSwarmAddCellDM(*nsw, ncelldm));
2906: PetscCall(DMSwarmCellDMDestroy(&ncelldm));
2907: }
2908: PetscCall(PetscFree(celldmnames));
2910: PetscCall(DMSetFromOptions(*nsw));
2911: PetscCall(DMSetUp(*nsw));
2912: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2913: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
2914: PetscCall(DMSwarmSetCellDMActive(*nsw, name));
2915: PetscFunctionReturn(PETSC_SUCCESS);
2916: }
2918: PetscErrorCode DMLocalToGlobalBegin_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
2919: {
2920: PetscFunctionBegin;
2921: PetscFunctionReturn(PETSC_SUCCESS);
2922: }
2924: PetscErrorCode DMLocalToGlobalEnd_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
2925: {
2926: PetscFunctionBegin;
2927: switch (mode) {
2928: case INSERT_VALUES:
2929: PetscCall(VecCopy(l, g));
2930: break;
2931: case ADD_VALUES:
2932: PetscCall(VecAXPY(g, 1., l));
2933: break;
2934: default:
2935: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mode not supported: %d", mode);
2936: }
2937: PetscFunctionReturn(PETSC_SUCCESS);
2938: }
2940: PetscErrorCode DMGlobalToLocalBegin_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
2941: {
2942: PetscFunctionBegin;
2943: PetscFunctionReturn(PETSC_SUCCESS);
2944: }
2946: PetscErrorCode DMGlobalToLocalEnd_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
2947: {
2948: PetscFunctionBegin;
2949: PetscCall(DMLocalToGlobalEnd_Swarm(dm, g, mode, l));
2950: PetscFunctionReturn(PETSC_SUCCESS);
2951: }