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 <petscdraw.h>
6: #include <petscdmplex.h>
7: #include <petscblaslapack.h>
8: #include "../src/dm/impls/swarm/data_bucket.h"
9: #include <petscdmlabel.h>
10: #include <petscsection.h>
12: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
13: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
14: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
16: const char *DMSwarmTypeNames[] = {"basic", "pic", NULL};
17: const char *DMSwarmMigrateTypeNames[] = {"basic", "dmcellnscatter", "dmcellexact", "user", NULL};
18: const char *DMSwarmCollectTypeNames[] = {"basic", "boundingbox", "general", "user", NULL};
19: const char *DMSwarmRemapTypeNames[] = {"none", "pfak", "colella", "DMSwarmRemapType", "DMSWARM_REMAP_", NULL};
20: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
22: const char DMSwarmField_pid[] = "DMSwarm_pid";
23: const char DMSwarmField_rank[] = "DMSwarm_rank";
24: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
26: PetscInt SwarmDataFieldId = -1;
28: #if PetscDefined(HAVE_HDF5)
29: #include <petscviewerhdf5.h>
31: static PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
32: {
33: DM dm;
34: PetscReal seqval;
35: PetscInt seqnum, bs;
36: PetscBool isseq, ists;
38: PetscFunctionBegin;
39: PetscCall(VecGetDM(v, &dm));
40: PetscCall(VecGetBlockSize(v, &bs));
41: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particle_fields"));
42: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
43: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
44: PetscCall(PetscViewerHDF5IsTimestepping(viewer, &ists));
45: if (ists) PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
46: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
47: PetscCall(VecViewNative(v, viewer));
48: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
49: PetscCall(PetscViewerHDF5PopGroup(viewer));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
54: {
55: DMSwarmCellDM celldm;
56: Vec coordinates;
57: PetscInt Np, Nfc;
58: PetscBool isseq;
59: const char **coordFields;
61: PetscFunctionBegin;
62: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
63: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
64: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
65: PetscCall(DMSwarmGetSize(dm, &Np));
66: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, coordFields[0], &coordinates));
67: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
68: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
69: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
70: PetscCall(VecViewNative(coordinates, viewer));
71: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
72: PetscCall(PetscViewerHDF5PopGroup(viewer));
73: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, coordFields[0], &coordinates));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
76: #endif
78: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
79: {
80: DM dm;
81: #if PetscDefined(HAVE_HDF5)
82: PetscBool ishdf5;
83: #endif
85: PetscFunctionBegin;
86: PetscCall(VecGetDM(v, &dm));
87: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
88: #if PetscDefined(HAVE_HDF5)
89: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
90: if (ishdf5) {
91: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
94: #endif
95: PetscCall(VecViewNative(v, viewer));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: /*@
100: DMSwarmVectorGetField - Gets the fields from which to define a `Vec` object
101: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
103: Not collective
105: Input Parameter:
106: . sw - a `DMSWARM`
108: Output Parameters:
109: + Nf - the number of fields
110: - fieldnames - the textual name given to each registered field, or NULL if it has not been set
112: Level: beginner
114: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
115: @*/
116: PetscErrorCode DMSwarmVectorGetField(DM sw, PetscInt *Nf, const char **fieldnames[])
117: {
118: DMSwarmCellDM celldm;
120: PetscFunctionBegin;
122: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
123: PetscCall(DMSwarmCellDMGetFields(celldm, Nf, fieldnames));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: /*@
128: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
129: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
131: Collective
133: Input Parameters:
134: + dm - a `DMSWARM`
135: - fieldname - the textual name given to each registered field
137: Level: beginner
139: Notes:
140: The field with name `fieldname` must have data type `PETSC_REAL` or `PETSC_SCALAR`.
142: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
143: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
145: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineFields()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
146: @*/
147: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
148: {
149: PetscFunctionBegin;
150: PetscCall(DMSwarmVectorDefineFields(dm, 1, &fieldname));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /*@
155: DMSwarmVectorDefineFields - Sets the fields from which to define a `Vec` object
156: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
158: Collective, No Fortran support
160: Input Parameters:
161: + sw - a `DMSWARM`
162: . Nf - the number of fields
163: - fieldnames - the textual name given to each registered field
165: Level: beginner
167: Notes:
168: Each field with name in `fieldnames` must have data type `PETSC_REAL` or `PETSC_SCALAR`.
170: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
171: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
173: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
174: @*/
175: PetscErrorCode DMSwarmVectorDefineFields(DM sw, PetscInt Nf, const char *fieldnames[])
176: {
177: DM_Swarm *swarm = (DM_Swarm *)sw->data;
178: DMSwarmCellDM celldm;
180: PetscFunctionBegin;
182: if (fieldnames) PetscAssertPointer(fieldnames, 3);
183: if (!swarm->issetup) PetscCall(DMSetUp(sw));
184: PetscCheck(Nf >= 0, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields must be non-negative, not %" PetscInt_FMT, Nf);
185: for (PetscInt f = 0; f < Nf; ++f) {
186: PetscDataType type;
188: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], NULL, &type));
189: PetscCheck(type == PETSC_REAL || type == PETSC_SCALAR, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Field %s must have type PETSC_REAL or PETSC_SCALAR, not %s", fieldnames[f], PetscDataTypes[type]);
190: }
191: // Create a dummy cell DM if none has been specified (I think we should not support this mode)
192: if (!swarm->activeCellDM) {
193: DM dm;
194: DMSwarmCellDM celldm;
196: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), &dm));
197: PetscCall(DMSetType(dm, DMSHELL));
198: PetscCall(PetscObjectSetName((PetscObject)dm, "dummy"));
199: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 0, NULL, &celldm));
200: PetscCall(DMDestroy(&dm));
201: PetscCall(DMSwarmAddCellDM(sw, celldm));
202: PetscCall(DMSwarmCellDMDestroy(&celldm));
203: PetscCall(DMSwarmSetCellDMActive(sw, "dummy"));
204: }
205: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
206: for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscFree(celldm->dmFields[f]));
207: PetscCall(PetscFree(celldm->dmFields));
209: celldm->Nf = Nf;
210: PetscCall(PetscMalloc1(Nf, &celldm->dmFields));
211: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscStrallocpy(fieldnames[f], (char **)&celldm->dmFields[f]));
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /* requires DMSwarmDefineFieldVector has been called */
216: static PetscErrorCode DMCreateGlobalVector_Swarm(DM sw, Vec *vec)
217: {
218: DM_Swarm *swarm = (DM_Swarm *)sw->data;
219: DMSwarmCellDM celldm;
220: Vec x;
221: char name[PETSC_MAX_PATH_LEN];
222: PetscInt bs = 0, n;
224: PetscFunctionBegin;
225: if (!swarm->issetup) PetscCall(DMSetUp(sw));
226: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
227: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
228: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
230: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
231: for (PetscInt f = 0; f < celldm->Nf; ++f) {
232: PetscInt fbs;
233: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
234: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
235: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
236: bs += fbs;
237: }
238: PetscCall(VecCreate(PetscObjectComm((PetscObject)sw), &x));
239: PetscCall(PetscObjectSetName((PetscObject)x, name));
240: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
241: PetscCall(VecSetBlockSize(x, bs));
242: PetscCall(VecSetDM(x, sw));
243: PetscCall(VecSetFromOptions(x));
244: PetscCall(VecSetOperation(x, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
245: *vec = x;
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: /* requires DMSwarmDefineFieldVector has been called */
250: static PetscErrorCode DMCreateLocalVector_Swarm(DM sw, Vec *vec)
251: {
252: DM_Swarm *swarm = (DM_Swarm *)sw->data;
253: DMSwarmCellDM celldm;
254: Vec x;
255: char name[PETSC_MAX_PATH_LEN];
256: PetscInt bs = 0, n;
258: PetscFunctionBegin;
259: if (!swarm->issetup) PetscCall(DMSetUp(sw));
260: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
261: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
262: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
264: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
265: for (PetscInt f = 0; f < celldm->Nf; ++f) {
266: PetscInt fbs;
267: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
268: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
269: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
270: bs += fbs;
271: }
272: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
273: PetscCall(PetscObjectSetName((PetscObject)x, name));
274: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
275: PetscCall(VecSetBlockSize(x, bs));
276: PetscCall(VecSetDM(x, sw));
277: PetscCall(VecSetFromOptions(x));
278: *vec = x;
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
283: {
284: DM_Swarm *swarm = (DM_Swarm *)dm->data;
285: DMSwarmDataField gfield;
286: PetscInt bs, nlocal, fid = -1, cfid = -2;
287: PetscBool flg;
289: PetscFunctionBegin;
290: /* check vector is an inplace array */
291: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
292: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
293: (void)flg; /* avoid compiler warning */
294: 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);
295: PetscCall(VecGetLocalSize(*vec, &nlocal));
296: PetscCall(VecGetBlockSize(*vec, &bs));
297: 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");
298: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
299: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
300: PetscCall(VecResetArray(*vec));
301: PetscCall(VecDestroy(vec));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
306: {
307: DM_Swarm *swarm = (DM_Swarm *)dm->data;
308: PetscDataType type;
309: PetscScalar *array;
310: PetscInt bs, n, fid;
311: char name[PETSC_MAX_PATH_LEN];
312: PetscMPIInt size;
313: PetscBool iscuda, iskokkos, iship;
315: PetscFunctionBegin;
316: if (!swarm->issetup) PetscCall(DMSetUp(dm));
317: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
318: PetscCall(DMSwarmGetFieldInfo(dm, fieldname, &bs, &type));
319: PetscCheck(type == PETSC_SCALAR, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Field %s must have type PETSC_SCALAR, not %s", fieldname, PetscDataTypes[type]);
320: PetscCall(DMSwarmGetField(dm, fieldname, NULL, NULL, (void **)&array));
322: PetscCallMPI(MPI_Comm_size(comm, &size));
323: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
324: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
325: PetscCall(PetscStrcmp(dm->vectype, VECHIP, &iship));
326: PetscCall(VecCreate(comm, vec));
327: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
328: PetscCall(VecSetBlockSize(*vec, bs));
329: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
330: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
331: else if (iship) PetscCall(VecSetType(*vec, VECHIP));
332: else PetscCall(VecSetType(*vec, VECSTANDARD));
333: PetscCall(VecPlaceArray(*vec, array));
335: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
336: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
338: /* Set guard */
339: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
340: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
342: PetscCall(VecSetDM(*vec, dm));
343: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: static PetscErrorCode DMSwarmDestroyVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], Vec *vec)
348: {
349: DM_Swarm *swarm = (DM_Swarm *)sw->data;
350: const PetscScalar *array;
351: PetscInt bs, n, id = 0, cid = -2;
352: PetscBool flg;
354: PetscFunctionBegin;
355: for (PetscInt f = 0; f < Nf; ++f) {
356: PetscInt fid;
358: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
359: id += fid;
360: }
361: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cid, flg));
362: (void)flg; /* avoid compiler warning */
363: 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);
364: PetscCall(VecGetLocalSize(*vec, &n));
365: PetscCall(VecGetBlockSize(*vec, &bs));
366: n /= bs;
367: PetscCheck(n == swarm->db->L, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
368: PetscCall(VecGetArrayRead(*vec, &array));
369: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
370: void *farray;
371: PetscDataType ftype;
372: PetscInt fbs;
374: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
375: PetscCheck(ftype == PETSC_REAL || ftype == PETSC_SCALAR, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Field %s must have type PETSC_REAL or PETSC_SCALAR, not %s", fieldnames[f], PetscDataTypes[ftype]);
376: PetscCall(DMSwarmGetField(sw, fieldnames[f], NULL, NULL, &farray));
377: PetscCheck(off + fbs <= bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid blocksize %" PetscInt_FMT " < %" PetscInt_FMT, bs, off + fbs);
378: if (ftype == PETSC_REAL) {
379: PetscReal *rarray = (PetscReal *)farray;
381: for (PetscInt i = 0; i < n; ++i) {
382: for (PetscInt b = 0; b < fbs; ++b) rarray[i * fbs + b] = PetscRealPart(array[i * bs + off + b]);
383: }
384: } else {
385: PetscScalar *sarray = (PetscScalar *)farray;
387: for (PetscInt i = 0; i < n; ++i) {
388: for (PetscInt b = 0; b < fbs; ++b) sarray[i * fbs + b] = array[i * bs + off + b];
389: }
390: }
391: off += fbs;
392: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], NULL, NULL, &farray));
393: }
394: PetscCall(VecRestoreArrayRead(*vec, &array));
395: PetscCall(VecDestroy(vec));
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: static PetscErrorCode DMSwarmCreateVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], MPI_Comm comm, Vec *vec)
400: {
401: DM_Swarm *swarm = (DM_Swarm *)sw->data;
402: PetscScalar *array;
403: PetscInt n, bs = 0, id = 0;
404: char name[PETSC_MAX_PATH_LEN];
406: PetscFunctionBegin;
407: if (!swarm->issetup) PetscCall(DMSetUp(sw));
408: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
409: for (PetscInt f = 0; f < Nf; ++f) {
410: PetscDataType ftype;
411: PetscInt fbs;
413: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
414: PetscCheck(ftype == PETSC_REAL || ftype == PETSC_SCALAR, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Field %s must have type PETSC_REAL or PETSC_SCALAR, not %s", fieldnames[f], PetscDataTypes[ftype]);
415: bs += fbs;
416: }
418: PetscCall(VecCreate(comm, vec));
419: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
420: PetscCall(VecSetBlockSize(*vec, bs));
421: PetscCall(VecSetType(*vec, sw->vectype));
423: PetscCall(VecGetArrayWrite(*vec, &array));
424: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
425: void *farray;
426: PetscDataType ftype;
427: PetscInt fbs;
429: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, &farray));
430: if (ftype == PETSC_REAL) {
431: const PetscReal *rarray = (const PetscReal *)farray;
433: for (PetscInt i = 0; i < n; ++i) {
434: for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = rarray[i * fbs + b];
435: }
436: } else {
437: const PetscScalar *sarray = (const PetscScalar *)farray;
439: for (PetscInt i = 0; i < n; ++i) {
440: for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = sarray[i * fbs + b];
441: }
442: }
443: off += fbs;
444: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, &farray));
445: }
446: PetscCall(VecRestoreArrayWrite(*vec, &array));
448: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
449: for (PetscInt f = 0; f < Nf; ++f) {
450: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
451: PetscCall(PetscStrlcat(name, fieldnames[f], PETSC_MAX_PATH_LEN));
452: }
453: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
455: for (PetscInt f = 0; f < Nf; ++f) {
456: PetscInt fid;
458: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
459: id += fid;
460: }
461: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, id));
463: PetscCall(VecSetDM(*vec, sw));
464: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: DMSwarmPreallocateMassMatrix - Preallocate particle mass matrix between a `DMSWARM` and a `DMPLEX`
471: Collective
473: Input Parameters:
474: + dmc - The `DMSWARM` object
475: . dmf - The `DMPLEX` object
476: . mass - The mass matrix to preallocate
477: . rStart - The starting row index for this process
478: . maxC - The maximum number of columns per row
479: - ctx - The user context
481: Level: developer
483: .seealso: [](ch_unstructured), `DM`, `DMSWARM`, `DMPLEX`, `DMCreateMassMatrix()`, `DMSwarmFillMassMatrix()`
484: @*/
485: PetscErrorCode DMSwarmPreallocateMassMatrix(DM dmc, DM dmf, Mat mass, PetscInt *rStart, PetscInt *maxC, PetscCtx ctx)
486: {
487: MPI_Comm comm;
488: PetscHSetIJ ht;
489: PetscDS ds;
490: PetscSection fsection, globalFSection;
491: PetscLayout rLayout, colLayout;
492: PetscInt *dnz, *onz;
493: PetscInt cStart, cEnd, locRows, locCols, colStart, colEnd, Nf, totNc = 0;
495: PetscFunctionBegin;
499: PetscAssertPointer(rStart, 4);
500: PetscAssertPointer(maxC, 5);
501: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
502: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
503: PetscCall(DMGetLocalSection(dmf, &fsection));
504: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
505: PetscCall(DMGetDS(dmf, &ds));
506: PetscCall(PetscDSGetNumFields(ds, &Nf));
507: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
508: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
509: PetscCall(PetscHSetIJCreate(&ht));
511: PetscCall(PetscLayoutCreate(comm, &colLayout));
512: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
513: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
514: PetscCall(PetscLayoutSetUp(colLayout));
515: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
516: PetscCall(PetscLayoutDestroy(&colLayout));
518: PetscCall(PetscLayoutCreate(comm, &rLayout));
519: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
520: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
521: PetscCall(PetscLayoutSetUp(rLayout));
522: PetscCall(PetscLayoutGetRange(rLayout, rStart, NULL));
523: PetscCall(PetscLayoutDestroy(&rLayout));
525: for (PetscInt field = 0; field < Nf; ++field) {
526: PetscObject obj;
527: PetscClassId id;
528: PetscInt Nc;
530: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
531: PetscCall(PetscObjectGetClassId(obj, &id));
532: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
533: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
534: totNc += Nc;
535: }
536: PetscCall(DMSwarmSortGetAccess(dmc));
537: for (PetscInt field = 0; field < Nf; ++field) {
538: PetscObject obj;
539: PetscClassId id;
540: PetscInt Nc;
542: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
543: PetscCall(PetscObjectGetClassId(obj, &id));
544: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
545: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
547: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
548: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
549: PetscInt numFIndices, numCIndices;
551: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
552: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
553: *maxC = PetscMax(*maxC, numCIndices);
554: {
555: PetscHashIJKey key;
556: PetscBool missing;
557: for (PetscInt i = 0; i < numFIndices; ++i) {
558: key.j = findices[i]; /* global column (from Plex) */
559: if (key.j >= 0) {
560: /* Get indices for coarse elements */
561: for (PetscInt j = 0; j < numCIndices; ++j) {
562: for (PetscInt c = 0; c < Nc; ++c) {
563: // TODO Need field offset on particle here
564: key.i = cindices[j] * totNc + c + *rStart; /* global cols (from Swarm) */
565: if (key.i < 0) continue;
566: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
567: PetscCheck(missing, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
568: if (key.j >= colStart && key.j < colEnd) ++dnz[key.i - *rStart];
569: else ++onz[key.i - *rStart];
570: }
571: }
572: }
573: }
574: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
575: }
576: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
577: }
578: }
579: PetscCall(DMSwarmSortRestoreAccess(dmc));
580: PetscCall(PetscHSetIJDestroy(&ht));
581: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
582: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
583: PetscCall(PetscFree2(dnz, onz));
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: /*@
588: DMSwarmFillMassMatrix - Assemble the particle mass matrix between a `DMSWARM` and a `DMPLEX`
590: Collective
592: Input Parameters:
593: + dmc - The `DMSWARM` object
594: . dmf - The `DMPLEX` object
595: . mass - The mass matrix to fill up
596: . rStart - The starting row index for this process
597: . maxC - The maximum number of columns per row
598: . useDeltaFunction - Flag to use a delta function for the particle shape function
599: . Nfc - The number of swarm coordinate fields
600: . bs - The block size for each swarm coordinate field
601: . coordVals - The values for each particle for each swarm coordinate field
602: - ctx - The user context
604: Level: developer
606: .seealso: [](ch_unstructured), `DM`, `DMSWARM`, `DMPLEX`, `DMCreateMassMatrix()`, `DMSwarmPreallocateMassMatrix()`
607: @*/
608: PetscErrorCode DMSwarmFillMassMatrix(DM dmc, DM dmf, Mat mass, PetscInt rStart, PetscInt maxC, PetscBool useDeltaFunction, PetscInt Nfc, const PetscInt bs[], PetscReal *coordVals[], PetscCtx ctx)
609: {
610: PetscDS ds;
611: PetscSection fsection, globalFSection;
612: PetscInt dim, cStart, cEnd, Nf, totDim, totNc = 0, *rowIDXs;
613: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
614: PetscScalar *elemMat;
616: PetscFunctionBegin;
620: PetscCall(DMGetCoordinateDim(dmf, &dim));
621: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
622: PetscCall(DMGetLocalSection(dmf, &fsection));
623: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
624: PetscCall(DMGetDS(dmf, &ds));
625: PetscCall(PetscDSGetNumFields(ds, &Nf));
626: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
627: for (PetscInt field = 0; field < Nf; ++field) {
628: PetscObject obj;
629: PetscClassId id;
630: PetscInt Nc;
632: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
633: PetscCall(PetscObjectGetClassId(obj, &id));
634: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
635: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
636: totNc += Nc;
637: }
639: PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
640: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
641: PetscCall(DMSwarmSortGetAccess(dmc));
642: for (PetscInt field = 0; field < Nf; ++field) {
643: PetscTabulation Tcoarse;
644: PetscObject obj;
645: PetscClassId id;
646: PetscInt Nc;
648: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
649: PetscCall(PetscObjectGetClassId(obj, &id));
650: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
651: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
653: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
654: PetscInt *findices, *cindices;
655: PetscInt numFIndices, numCIndices;
657: /* TODO: Use DMField instead of assuming affine */
658: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
659: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
660: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
661: for (PetscInt j = 0; j < numCIndices; ++j) {
662: PetscReal xr[8];
663: PetscInt off = 0;
665: for (PetscInt i = 0; i < Nfc; ++i) {
666: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[j] * bs[i] + b];
667: }
668: 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);
669: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, xr, &xi[j * dim]);
670: }
671: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
672: else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
673: /* Get elemMat entries by multiplying by weight */
674: PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
675: for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
676: for (PetscInt j = 0; j < numCIndices; ++j) {
677: for (PetscInt c = 0; c < Nc; ++c) {
678: // TODO Need field offset on particle and field here
679: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
680: elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
681: }
682: }
683: }
684: for (PetscInt j = 0; j < numCIndices; ++j)
685: // TODO Need field offset on particle here
686: for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
687: if (0) PetscCall(DMPrintCellMatrix(cell, "Mass Matrix", numCIndices * Nc, numFIndices, elemMat));
688: PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
689: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
690: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
691: PetscCall(PetscTabulationDestroy(&Tcoarse));
692: }
693: }
694: PetscCall(DMSwarmSortRestoreAccess(dmc));
695: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
696: PetscCall(PetscFree3(v0, J, invJ));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
702: \hat f = \sum_i f_i \phi_i
704: and a particle function is given by
706: f = \sum_i w_i \delta(x - x_i)
708: then we want to require that
710: M \hat f = M_p f
712: where the particle mass matrix is given by
714: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
716: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
717: his integral. We allow this with the boolean flag.
718: */
719: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, PetscCtx ctx)
720: {
721: DMSwarmCellDM celldm;
722: PetscInt rStart, maxC = 0;
723: PetscInt Nfc;
724: const char **coordFields;
725: PetscReal **coordVals;
726: PetscInt *bs;
728: PetscFunctionBegin;
729: PetscCall(DMSwarmPreallocateMassMatrix(dmc, dmf, mass, &rStart, &maxC, ctx));
731: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
732: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
733: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
734: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
735: PetscCall(DMSwarmFillMassMatrix(dmc, dmf, mass, rStart, maxC, useDeltaFunction, Nfc, bs, coordVals, ctx));
736: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
737: PetscCall(PetscFree2(coordVals, bs));
739: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
740: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: /* Returns empty matrix for use with SNES FD */
745: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
746: {
747: Vec field;
748: PetscInt size;
750: PetscFunctionBegin;
751: PetscCall(DMGetGlobalVector(sw, &field));
752: PetscCall(VecGetLocalSize(field, &size));
753: PetscCall(DMRestoreGlobalVector(sw, &field));
754: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
755: PetscCall(MatSetFromOptions(*m));
756: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
757: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
758: PetscCall(MatZeroEntries(*m));
759: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
760: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
761: PetscCall(MatShift(*m, 1.0));
762: PetscCall(MatSetDM(*m, sw));
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: /* FEM cols, Particle rows */
767: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
768: {
769: DMSwarmCellDM celldm;
770: PetscSection gsf;
771: PetscInt m, n, Np, bs;
772: void *ctx;
774: PetscFunctionBegin;
775: PetscCall(DMSwarmGetCellDMActive(dmCoarse, &celldm));
776: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)dmCoarse), PETSC_ERR_USER, "Active cell DM does not define any fields");
777: PetscCall(DMGetGlobalSection(dmFine, &gsf));
778: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
779: PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
780: PetscCall(DMSwarmCellDMGetBlockSize(celldm, dmCoarse, &bs));
781: n = Np * bs;
782: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
783: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
784: PetscCall(MatSetType(*mass, dmCoarse->mattype));
785: PetscCall(DMGetApplicationContext(dmFine, &ctx));
787: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
788: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
789: PetscFunctionReturn(PETSC_SUCCESS);
790: }
792: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, PetscCtx ctx)
793: {
794: MPI_Comm comm;
795: DMSwarmCellDM celldm;
796: PetscDS prob;
797: PetscSection fsection, globalFSection;
798: PetscHSetIJ ht;
799: PetscLayout rLayout, colLayout;
800: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
801: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
802: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
803: PetscScalar *elemMat, *elemMatSq;
804: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0;
805: const char **coordFields;
806: PetscReal **coordVals;
807: PetscInt *bs;
809: PetscFunctionBegin;
810: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
811: PetscCall(DMGetCoordinateDim(dmf, &cdim));
812: PetscCall(DMGetDS(dmf, &prob));
813: PetscCall(PetscDSGetNumFields(prob, &Nf));
814: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
815: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
816: PetscCall(DMGetLocalSection(dmf, &fsection));
817: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
818: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
819: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
821: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
822: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
823: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
825: PetscCall(PetscLayoutCreate(comm, &colLayout));
826: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
827: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
828: PetscCall(PetscLayoutSetUp(colLayout));
829: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
830: PetscCall(PetscLayoutDestroy(&colLayout));
832: PetscCall(PetscLayoutCreate(comm, &rLayout));
833: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
834: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
835: PetscCall(PetscLayoutSetUp(rLayout));
836: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
837: PetscCall(PetscLayoutDestroy(&rLayout));
839: PetscCall(DMPlexGetDepth(dmf, &depth));
840: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
841: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
842: PetscCall(PetscMalloc1(maxAdjSize, &adj));
844: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
845: PetscCall(PetscHSetIJCreate(&ht));
846: /* Count nonzeros
847: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
848: */
849: PetscCall(DMSwarmSortGetAccess(dmc));
850: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
851: PetscInt *cindices;
852: PetscInt numCIndices;
853: #if 0
854: PetscInt adjSize = maxAdjSize, a, j;
855: #endif
857: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
858: maxC = PetscMax(maxC, numCIndices);
859: /* Diagonal block */
860: for (PetscInt i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
861: #if 0
862: /* Off-diagonal blocks */
863: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
864: for (a = 0; a < adjSize; ++a) {
865: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
866: const PetscInt ncell = adj[a];
867: PetscInt *ncindices;
868: PetscInt numNCIndices;
870: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
871: {
872: PetscHashIJKey key;
873: PetscBool missing;
875: for (i = 0; i < numCIndices; ++i) {
876: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
877: if (key.i < 0) continue;
878: for (j = 0; j < numNCIndices; ++j) {
879: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
880: if (key.j < 0) continue;
881: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
882: if (missing) {
883: if (key.j >= colStart && key.j < colEnd) ++dnz[key.i - rStart];
884: else ++onz[key.i - rStart];
885: }
886: }
887: }
888: }
889: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
890: }
891: }
892: #endif
893: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
894: }
895: PetscCall(PetscHSetIJDestroy(&ht));
896: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
897: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
898: PetscCall(PetscFree2(dnz, onz));
899: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
900: /* Fill in values
901: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
902: Start just by producing block diagonal
903: Could loop over adjacent cells
904: Produce neighboring element matrix
905: TODO Determine which columns and rows correspond to shared dual vector
906: Do MatMatMult with rectangular matrices
907: Insert block
908: */
909: for (PetscInt field = 0; field < Nf; ++field) {
910: PetscTabulation Tcoarse;
911: PetscObject obj;
912: PetscInt Nc;
914: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
915: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
916: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
917: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
918: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
919: PetscInt *findices, *cindices;
920: PetscInt numFIndices, numCIndices;
922: /* TODO: Use DMField instead of assuming affine */
923: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
924: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
925: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
926: for (PetscInt p = 0; p < numCIndices; ++p) {
927: PetscReal xr[8];
928: PetscInt off = 0;
930: for (PetscInt i = 0; i < Nfc; ++i) {
931: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[p] * bs[i] + b];
932: }
933: 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);
934: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[p * cdim]);
935: }
936: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
937: /* Get elemMat entries by multiplying by weight */
938: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
939: for (PetscInt i = 0; i < numFIndices; ++i) {
940: for (PetscInt p = 0; p < numCIndices; ++p) {
941: for (PetscInt c = 0; c < Nc; ++c) {
942: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
943: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
944: }
945: }
946: }
947: PetscCall(PetscTabulationDestroy(&Tcoarse));
948: for (PetscInt p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
949: if (0) PetscCall(DMPrintCellMatrix(cell, "Mass Matrix Square", 1, numCIndices, elemMat));
950: /* Block diagonal */
951: if (numCIndices) {
952: PetscBLASInt blasn, blask;
953: PetscScalar one = 1.0, zero = 0.0;
955: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
956: PetscCall(PetscBLASIntCast(numFIndices, &blask));
957: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
958: }
959: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
960: /* TODO off-diagonal */
961: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
962: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
963: }
964: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
965: }
966: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
967: PetscCall(PetscFree(adj));
968: PetscCall(DMSwarmSortRestoreAccess(dmc));
969: PetscCall(PetscFree3(v0, J, invJ));
970: PetscCall(PetscFree2(coordVals, bs));
971: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
972: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: /*@
977: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
979: Collective
981: Input Parameters:
982: + dmCoarse - a `DMSWARM`
983: - dmFine - a `DMPLEX`
985: Output Parameter:
986: . mass - the square of the particle mass matrix
988: Level: advanced
990: Note:
991: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
992: future to compute the full normal equations.
994: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
995: @*/
996: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
997: {
998: PetscInt n;
999: void *ctx;
1001: PetscFunctionBegin;
1002: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
1003: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
1004: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
1005: PetscCall(MatSetType(*mass, dmCoarse->mattype));
1006: PetscCall(DMGetApplicationContext(dmFine, &ctx));
1008: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
1009: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: /* 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
1015: \int_X \psi_i \nabla \cdot \hat f = \int_X \psi_i \nabla \cdot f
1017: and then integrate by parts
1019: \int_X \nabla \psi_i \cdot \hat f = \int_X \nabla \psi_i \cdot f
1021: where \psi is from a scalar FE space. If a finite element interpolant is given by
1023: \hat f^c = \sum_i f_i \phi^c_i
1025: and a particle function is given by
1027: f^c = \sum_p f^c_p \delta(x - x_p)
1029: then we want to require that
1031: D_f \hat f = D_p f
1033: where the gradient matrices are given by
1035: (D_f)_{i(jc)} = \int \partial_c \psi_i \phi_j
1036: (D_p)_{i(jc)} = \int \partial_c \psi_i \delta(x - x_j)
1038: Thus we need two finite element spaces, a scalar and a vector. The vector space holds the representer for the
1039: vector particle field. The scalar space holds the output of D_p or D_f, which is the weak divergence of the field.
1041: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
1042: his integral. We allow this with the boolean flag.
1043: */
1044: static PetscErrorCode DMSwarmComputeGradientMatrix_Private(DM sw, DM dm, Mat derv, PetscBool useDeltaFunction, PetscCtx ctx)
1045: {
1046: MPI_Comm comm;
1047: DMSwarmCellDM celldm;
1048: PetscDS ds;
1049: PetscSection fsection, globalFSection;
1050: PetscLayout rLayout;
1051: PetscInt locRows, rStart, *rowIDXs;
1052: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
1053: PetscScalar *elemMat;
1054: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxNpc = 0, totNc = 0;
1055: const char **coordFields;
1056: PetscReal **coordVals;
1057: PetscInt *bs;
1059: PetscFunctionBegin;
1060: PetscCall(PetscObjectGetComm((PetscObject)derv, &comm));
1061: PetscCall(DMGetCoordinateDim(dm, &cdim));
1062: PetscCall(DMGetDS(dm, &ds));
1063: PetscCall(PetscDSGetNumFields(ds, &Nf));
1064: PetscCheck(Nf == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
1065: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
1066: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
1067: PetscCall(DMGetLocalSection(dm, &fsection));
1068: PetscCall(DMGetGlobalSection(dm, &globalFSection));
1069: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1070: PetscCall(MatGetLocalSize(derv, &locRows, NULL));
1072: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1073: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1074: PetscCheck(Nfc == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
1075: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
1077: PetscCall(PetscLayoutCreate(comm, &rLayout));
1078: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
1079: PetscCall(PetscLayoutSetBlockSize(rLayout, cdim));
1080: PetscCall(PetscLayoutSetUp(rLayout));
1081: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
1082: PetscCall(PetscLayoutDestroy(&rLayout));
1084: for (PetscInt field = 0; field < Nf; ++field) {
1085: PetscObject obj;
1086: PetscInt Nc;
1088: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
1089: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
1090: totNc += Nc;
1091: }
1092: PetscCheck(totNc == 1, comm, PETSC_ERR_ARG_WRONG, "The number of field components %" PetscInt_FMT " != 1", totNc);
1093: /* count non-zeros */
1094: PetscCall(DMSwarmSortGetAccess(sw));
1095: for (PetscInt field = 0; field < Nf; ++field) {
1096: PetscObject obj;
1097: PetscInt Nc;
1099: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
1100: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
1101: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1102: PetscInt *pind;
1103: PetscInt Npc;
1105: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1106: maxNpc = PetscMax(maxNpc, Npc);
1107: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1108: }
1109: }
1110: PetscCall(PetscMalloc3(maxNpc * cdim * totDim, &elemMat, maxNpc * cdim, &rowIDXs, maxNpc * cdim, &xi));
1111: for (PetscInt field = 0; field < Nf; ++field) {
1112: PetscTabulation Tcoarse;
1113: PetscFE fe;
1115: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1116: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1117: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1118: PetscInt *findices, *pind;
1119: PetscInt numFIndices, Npc;
1121: /* TODO: Use DMField instead of assuming affine */
1122: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
1123: PetscCall(DMPlexGetClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1124: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1125: for (PetscInt j = 0; j < Npc; ++j) {
1126: PetscReal xr[8];
1127: PetscInt off = 0;
1129: for (PetscInt i = 0; i < Nfc; ++i) {
1130: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][pind[j] * bs[i] + b];
1131: }
1132: 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);
1133: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[j * cdim]);
1134: }
1135: PetscCall(PetscFECreateTabulation(fe, 1, Npc, xi, 1, &Tcoarse));
1136: /* Get elemMat entries by multiplying by weight */
1137: PetscCall(PetscArrayzero(elemMat, Npc * cdim * totDim));
1138: for (PetscInt i = 0; i < numFIndices; ++i) {
1139: for (PetscInt j = 0; j < Npc; ++j) {
1140: /* D[((p*pdim + i)*Nc + c)*cdim + d] is the value at point p for basis function i, component c, derivative d */
1141: for (PetscInt d = 0; d < cdim; ++d) {
1142: xi[d] = 0.;
1143: for (PetscInt e = 0; e < cdim; ++e) xi[d] += invJ[e * cdim + d] * Tcoarse->T[1][(j * numFIndices + i) * cdim + e];
1144: elemMat[(j * cdim + d) * numFIndices + i] += xi[d] * (useDeltaFunction ? 1.0 : detJ);
1145: }
1146: }
1147: }
1148: for (PetscInt j = 0; j < Npc; ++j)
1149: for (PetscInt d = 0; d < cdim; ++d) rowIDXs[j * cdim + d] = pind[j] * cdim + d + rStart;
1150: if (0) PetscCall(DMPrintCellMatrix(cell, "Derivative Matrix", Npc * cdim, numFIndices, elemMat));
1151: PetscCall(MatSetValues(derv, Npc * cdim, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
1152: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1153: PetscCall(DMPlexRestoreClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1154: PetscCall(PetscTabulationDestroy(&Tcoarse));
1155: }
1156: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1157: }
1158: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
1159: PetscCall(DMSwarmSortRestoreAccess(sw));
1160: PetscCall(PetscFree3(v0, J, invJ));
1161: PetscCall(PetscFree2(coordVals, bs));
1162: PetscCall(MatAssemblyBegin(derv, MAT_FINAL_ASSEMBLY));
1163: PetscCall(MatAssemblyEnd(derv, MAT_FINAL_ASSEMBLY));
1164: PetscFunctionReturn(PETSC_SUCCESS);
1165: }
1167: /* FEM cols: this is a scalar space
1168: Particle rows: this is a vector space that contracts with the derivative
1169: */
1170: static PetscErrorCode DMCreateGradientMatrix_Swarm(DM sw, DM dm, Mat *derv)
1171: {
1172: DMSwarmCellDM celldm;
1173: PetscSection gs;
1174: PetscInt cdim, m, n, Np, bs;
1175: void *ctx;
1176: MPI_Comm comm;
1178: PetscFunctionBegin;
1179: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1180: PetscCall(DMGetCoordinateDim(dm, &cdim));
1181: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1182: PetscCheck(celldm->Nf, comm, PETSC_ERR_USER, "Active cell DM does not define any fields");
1183: PetscCall(DMGetGlobalSection(dm, &gs));
1184: PetscCall(PetscSectionGetConstrainedStorageSize(gs, &n));
1185: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1186: PetscCall(DMSwarmCellDMGetBlockSize(celldm, sw, &bs));
1187: PetscCheck(cdim == bs, comm, PETSC_ERR_ARG_WRONG, "Coordinate dimension %" PetscInt_FMT " != %" PetscInt_FMT " swarm field block size", cdim, bs);
1188: m = Np * bs;
1189: PetscCall(MatCreate(PetscObjectComm((PetscObject)sw), derv));
1190: PetscCall(PetscObjectSetName((PetscObject)*derv, "Swarm Derivative Matrix"));
1191: PetscCall(MatSetSizes(*derv, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
1192: PetscCall(MatSetType(*derv, sw->mattype));
1193: PetscCall(DMGetApplicationContext(dm, &ctx));
1195: PetscCall(DMSwarmComputeGradientMatrix_Private(sw, dm, *derv, PETSC_TRUE, ctx));
1196: PetscCall(MatViewFromOptions(*derv, NULL, "-gradient_mat_view"));
1197: PetscFunctionReturn(PETSC_SUCCESS);
1198: }
1200: /*@
1201: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
1203: Collective
1205: Input Parameters:
1206: + dm - a `DMSWARM`
1207: - fieldname - the textual name given to a registered field
1209: Output Parameter:
1210: . vec - the vector
1212: Level: beginner
1214: Note:
1215: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
1217: The field must have data type `PETSC_SCALAR`.
1219: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
1220: @*/
1221: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1222: {
1223: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1225: PetscFunctionBegin;
1227: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: /*@
1232: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
1234: Collective
1236: Input Parameters:
1237: + dm - a `DMSWARM`
1238: - fieldname - the textual name given to a registered field
1240: Output Parameter:
1241: . vec - the vector
1243: Level: beginner
1245: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1246: @*/
1247: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1248: {
1249: PetscFunctionBegin;
1251: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: /*@
1256: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
1258: Collective
1260: Input Parameters:
1261: + dm - a `DMSWARM`
1262: - fieldname - the textual name given to a registered field
1264: Output Parameter:
1265: . vec - the vector
1267: Level: beginner
1269: Note:
1270: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1272: The field must have data type `PETSC_SCALAR`.
1274: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1275: @*/
1276: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1277: {
1278: PetscFunctionBegin;
1279: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, PETSC_COMM_SELF, vec));
1280: PetscFunctionReturn(PETSC_SUCCESS);
1281: }
1283: /*@
1284: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
1286: Collective
1288: Input Parameters:
1289: + dm - a `DMSWARM`
1290: - fieldname - the textual name given to a registered field
1292: Output Parameter:
1293: . vec - the vector
1295: Level: beginner
1297: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
1298: @*/
1299: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1300: {
1301: PetscFunctionBegin;
1303: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1304: PetscFunctionReturn(PETSC_SUCCESS);
1305: }
1307: /*@
1308: DMSwarmCreateGlobalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1310: Collective
1312: Input Parameters:
1313: + dm - a `DMSWARM`
1314: . Nf - the number of fields
1315: - fieldnames - the textual names given to the registered fields
1317: Output Parameter:
1318: . vec - the vector
1320: Level: beginner
1322: Notes:
1323: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromFields()`.
1325: All fields must have data type `PETSC_REAL` or `PETSC_SCALAR`. When PETSc is configured with complex scalars,
1326: only the real part of vector entries is copied back to `PETSC_REAL` fields.
1328: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateGlobalVectorFromField()`
1330: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromFields()`
1331: @*/
1332: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1333: {
1334: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1336: PetscFunctionBegin;
1338: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: /*@
1343: DMSwarmDestroyGlobalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1345: Collective
1347: Input Parameters:
1348: + dm - a `DMSWARM`
1349: . Nf - the number of fields
1350: - fieldnames - the textual names given to the registered fields
1352: Output Parameter:
1353: . vec - the vector
1355: Level: beginner
1357: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1358: @*/
1359: PetscErrorCode DMSwarmDestroyGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1360: {
1361: PetscFunctionBegin;
1363: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1364: PetscFunctionReturn(PETSC_SUCCESS);
1365: }
1367: /*@
1368: DMSwarmCreateLocalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1370: Collective
1372: Input Parameters:
1373: + dm - a `DMSWARM`
1374: . Nf - the number of fields
1375: - fieldnames - the textual names given to the registered fields
1377: Output Parameter:
1378: . vec - the vector
1380: Level: beginner
1382: Notes:
1383: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1385: All fields must have data type `PETSC_REAL` or `PETSC_SCALAR`. When PETSc is configured with complex scalars,
1386: only the real part of vector entries is copied back to `PETSC_REAL` fields.
1388: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateLocalVectorFromField()`
1390: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1391: @*/
1392: PetscErrorCode DMSwarmCreateLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1393: {
1394: PetscFunctionBegin;
1395: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, PETSC_COMM_SELF, vec));
1396: PetscFunctionReturn(PETSC_SUCCESS);
1397: }
1399: /*@
1400: DMSwarmDestroyLocalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1402: Collective
1404: Input Parameters:
1405: + dm - a `DMSWARM`
1406: . Nf - the number of fields
1407: - fieldnames - the textual names given to the registered fields
1409: Output Parameter:
1410: . vec - the vector
1412: Level: beginner
1414: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromFields()`
1415: @*/
1416: PetscErrorCode DMSwarmDestroyLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1417: {
1418: PetscFunctionBegin;
1420: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1421: PetscFunctionReturn(PETSC_SUCCESS);
1422: }
1424: /*@
1425: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
1427: Collective
1429: Input Parameter:
1430: . dm - a `DMSWARM`
1432: Level: beginner
1434: Note:
1435: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
1437: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1438: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1439: @*/
1440: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
1441: {
1442: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1444: PetscFunctionBegin;
1445: if (!swarm->field_registration_initialized) {
1446: swarm->field_registration_initialized = PETSC_TRUE;
1447: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
1448: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
1449: }
1450: PetscFunctionReturn(PETSC_SUCCESS);
1451: }
1453: /*@
1454: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
1456: Collective
1458: Input Parameter:
1459: . dm - a `DMSWARM`
1461: Level: beginner
1463: Note:
1464: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
1466: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1467: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1468: @*/
1469: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
1470: {
1471: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1473: PetscFunctionBegin;
1474: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
1475: swarm->field_registration_finalized = PETSC_TRUE;
1476: PetscFunctionReturn(PETSC_SUCCESS);
1477: }
1479: /*@
1480: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
1482: Not Collective
1484: Input Parameters:
1485: + sw - a `DMSWARM`
1486: . nlocal - the length of each registered field
1487: - buffer - the length of the buffer used to efficient dynamic re-sizing
1489: Level: beginner
1491: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
1492: @*/
1493: PetscErrorCode DMSwarmSetLocalSizes(DM sw, PetscInt nlocal, PetscInt buffer)
1494: {
1495: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1496: PetscMPIInt rank;
1497: PetscInt *rankval;
1499: PetscFunctionBegin;
1500: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
1501: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
1502: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
1504: // Initialize values in pid and rank placeholders
1505: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1506: PetscCall(DMSwarmGetField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1507: for (PetscInt p = 0; p < nlocal; p++) rankval[p] = rank;
1508: PetscCall(DMSwarmRestoreField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1509: /* TODO: [pid - use MPI_Scan] */
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }
1513: /*@
1514: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
1516: Collective
1518: Input Parameters:
1519: + sw - a `DMSWARM`
1520: - dm - the `DM` to attach to the `DMSWARM`
1522: Level: beginner
1524: Note:
1525: The attached `DM` (dm) will be queried for point location and
1526: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1528: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
1529: @*/
1530: PetscErrorCode DMSwarmSetCellDM(DM sw, DM dm)
1531: {
1532: DMSwarmCellDM celldm;
1533: const char *name;
1534: char *coordName;
1536: PetscFunctionBegin;
1539: PetscCall(PetscStrallocpy(DMSwarmPICField_coor, &coordName));
1540: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 1, (const char **)&coordName, &celldm));
1541: PetscCall(PetscFree(coordName));
1542: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1543: PetscCall(DMSwarmAddCellDM(sw, celldm));
1544: PetscCall(DMSwarmCellDMDestroy(&celldm));
1545: PetscCall(DMSwarmSetCellDMActive(sw, name));
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: /*@
1550: DMSwarmGetCellDM - Fetches the active cell `DM`
1552: Collective
1554: Input Parameter:
1555: . sw - a `DMSWARM`
1557: Output Parameter:
1558: . dm - the active `DM` for the `DMSWARM`
1560: Level: beginner
1562: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1563: @*/
1564: PetscErrorCode DMSwarmGetCellDM(DM sw, DM *dm)
1565: {
1566: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1567: DMSwarmCellDM celldm;
1569: PetscFunctionBegin;
1571: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)&celldm));
1572: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "There is no cell DM named %s in this Swarm", swarm->activeCellDM);
1573: PetscCall(DMSwarmCellDMGetDM(celldm, dm));
1574: PetscFunctionReturn(PETSC_SUCCESS);
1575: }
1577: /*@
1578: DMSwarmGetCellDMNames - Get the list of cell `DM` names
1580: Not collective
1582: Input Parameter:
1583: . sw - a `DMSWARM`
1585: Output Parameters:
1586: + Ndm - the number of `DMSwarmCellDM` in the `DMSWARM`
1587: - celldms - the name of each `DMSwarmCellDM`
1589: Level: beginner
1591: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMByName()`
1592: @*/
1593: PetscErrorCode DMSwarmGetCellDMNames(DM sw, PetscInt *Ndm, const char **celldms[])
1594: {
1595: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1596: PetscObjectList next = swarm->cellDMs;
1597: PetscInt n = 0;
1599: PetscFunctionBegin;
1601: PetscAssertPointer(Ndm, 2);
1602: PetscAssertPointer(celldms, 3);
1603: while (next) {
1604: next = next->next;
1605: ++n;
1606: }
1607: PetscCall(PetscMalloc1(n, celldms));
1608: next = swarm->cellDMs;
1609: n = 0;
1610: while (next) {
1611: (*celldms)[n] = (const char *)next->obj->name;
1612: next = next->next;
1613: ++n;
1614: }
1615: *Ndm = n;
1616: PetscFunctionReturn(PETSC_SUCCESS);
1617: }
1619: /*@
1620: DMSwarmSetCellDMActive - Activates a cell `DM` for a `DMSWARM`
1622: Collective
1624: Input Parameters:
1625: + sw - a `DMSWARM`
1626: - name - name of the cell `DM` to active for the `DMSWARM`
1628: Level: beginner
1630: Note:
1631: The attached `DM` (dmcell) will be queried for point location and
1632: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1634: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1635: @*/
1636: PetscErrorCode DMSwarmSetCellDMActive(DM sw, const char name[])
1637: {
1638: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1639: DMSwarmCellDM celldm;
1641: PetscFunctionBegin;
1643: PetscCall(PetscInfo(sw, "Setting cell DM to %s\n", name));
1644: PetscCall(PetscFree(swarm->activeCellDM));
1645: PetscCall(PetscStrallocpy(name, (char **)&swarm->activeCellDM));
1646: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: /*@
1651: DMSwarmGetCellDMActive - Returns the active cell `DM` for a `DMSWARM`
1653: Collective
1655: Input Parameter:
1656: . sw - a `DMSWARM`
1658: Output Parameter:
1659: . celldm - the active `DMSwarmCellDM`
1661: Level: beginner
1663: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1664: @*/
1665: PetscErrorCode DMSwarmGetCellDMActive(DM sw, DMSwarmCellDM *celldm)
1666: {
1667: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1669: PetscFunctionBegin;
1671: PetscAssertPointer(celldm, 2);
1672: PetscCheck(swarm->activeCellDM, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no active cell DM");
1673: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)celldm));
1674: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", swarm->activeCellDM);
1675: PetscFunctionReturn(PETSC_SUCCESS);
1676: }
1678: /*@
1679: DMSwarmGetCellDMByName - Get a `DMSwarmCellDM` from its name
1681: Not collective
1683: Input Parameters:
1684: + sw - a `DMSWARM`
1685: - name - the name
1687: Output Parameter:
1688: . celldm - the `DMSwarmCellDM`, or `NULL` if the name is unknown
1690: Level: beginner
1692: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMNames()`
1693: @*/
1694: PetscErrorCode DMSwarmGetCellDMByName(DM sw, const char name[], DMSwarmCellDM *celldm)
1695: {
1696: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1698: PetscFunctionBegin;
1700: PetscAssertPointer(name, 2);
1701: PetscAssertPointer(celldm, 3);
1702: PetscCall(PetscObjectListFind(swarm->cellDMs, name, (PetscObject *)celldm));
1703: PetscFunctionReturn(PETSC_SUCCESS);
1704: }
1706: /*@
1707: DMSwarmAddCellDM - Adds a cell `DM` to the `DMSWARM`
1709: Collective
1711: Input Parameters:
1712: + sw - a `DMSWARM`
1713: - celldm - the `DMSwarmCellDM`
1715: Level: beginner
1717: Note:
1718: Cell DMs with the same name will share the cellid field
1720: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPushCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1721: @*/
1722: PetscErrorCode DMSwarmAddCellDM(DM sw, DMSwarmCellDM celldm)
1723: {
1724: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1725: const char *name;
1726: PetscInt dim;
1727: PetscBool flg;
1728: MPI_Comm comm;
1730: PetscFunctionBegin;
1732: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1734: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1735: PetscCall(PetscObjectListAdd(&swarm->cellDMs, name, (PetscObject)celldm));
1736: PetscCall(DMGetDimension(sw, &dim));
1737: for (PetscInt f = 0; f < celldm->Nfc; ++f) {
1738: PetscCall(DMSwarmDataFieldStringInList(celldm->coordFields[f], swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1739: if (!flg) {
1740: PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->coordFields[f], dim, PETSC_REAL));
1741: } else {
1742: PetscDataType dt;
1743: PetscInt bs;
1745: PetscCall(DMSwarmGetFieldInfo(sw, celldm->coordFields[f], &bs, &dt));
1746: PetscCheck(bs == dim, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has blocksize %" PetscInt_FMT " != %" PetscInt_FMT " spatial dimension", celldm->coordFields[f], bs, dim);
1747: PetscCheck(dt == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has datatype %s != PETSC_REAL", celldm->coordFields[f], PetscDataTypes[dt]);
1748: }
1749: }
1750: // Assume that DMs with the same name share the cellid field
1751: PetscCall(DMSwarmDataFieldStringInList(celldm->cellid, swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1752: if (!flg) {
1753: PetscBool isShell, isDummy;
1754: const char *name;
1756: // Allow dummy DMSHELL (I don't think we should support this mode)
1757: PetscCall(PetscObjectTypeCompare((PetscObject)celldm->dm, DMSHELL, &isShell));
1758: PetscCall(PetscObjectGetName((PetscObject)celldm->dm, &name));
1759: PetscCall(PetscStrcmp(name, "dummy", &isDummy));
1760: if (!isShell || !isDummy) PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->cellid, 1, PETSC_INT));
1761: }
1762: PetscCall(DMSwarmSetCellDMActive(sw, name));
1763: PetscFunctionReturn(PETSC_SUCCESS);
1764: }
1766: /*@
1767: DMSwarmGetLocalSize - Retrieves the local length of fields registered
1769: Not Collective
1771: Input Parameter:
1772: . dm - a `DMSWARM`
1774: Output Parameter:
1775: . nlocal - the length of each registered field
1777: Level: beginner
1779: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
1780: @*/
1781: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
1782: {
1783: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1785: PetscFunctionBegin;
1786: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
1787: PetscFunctionReturn(PETSC_SUCCESS);
1788: }
1790: /*@
1791: DMSwarmGetSize - Retrieves the total length of fields registered
1793: Collective
1795: Input Parameter:
1796: . dm - a `DMSWARM`
1798: Output Parameter:
1799: . n - the total length of each registered field
1801: Level: beginner
1803: Note:
1804: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
1806: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
1807: @*/
1808: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
1809: {
1810: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1812: PetscFunctionBegin;
1813: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, n, NULL, NULL));
1814: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1815: PetscFunctionReturn(PETSC_SUCCESS);
1816: }
1818: /*@
1819: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
1821: Collective
1823: Input Parameters:
1824: + dm - a `DMSWARM`
1825: . fieldname - the textual name to identify this field
1826: . blocksize - the number of each data type
1827: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
1829: Level: beginner
1831: Notes:
1832: The textual name for each registered field must be unique.
1834: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1835: @*/
1836: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1837: {
1838: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1839: size_t size;
1841: PetscFunctionBegin;
1842: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1843: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1845: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1846: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1847: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1848: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1849: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1851: PetscCall(PetscDataTypeGetSize(type, &size));
1852: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1853: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1854: {
1855: DMSwarmDataField gfield;
1857: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1858: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1859: }
1860: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1861: PetscFunctionReturn(PETSC_SUCCESS);
1862: }
1864: /*@
1865: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
1867: Collective
1869: Input Parameters:
1870: + dm - a `DMSWARM`
1871: . fieldname - the textual name to identify this field
1872: - size - the size in bytes of the user struct of each data type
1874: Level: beginner
1876: Note:
1877: The textual name for each registered field must be unique.
1879: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1880: @*/
1881: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1882: {
1883: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1885: PetscFunctionBegin;
1886: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1887: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1888: PetscFunctionReturn(PETSC_SUCCESS);
1889: }
1891: /*@
1892: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1894: Collective
1896: Input Parameters:
1897: + dm - a `DMSWARM`
1898: . fieldname - the textual name to identify this field
1899: . size - the size in bytes of the user data type
1900: - blocksize - the number of each data type
1902: Level: beginner
1904: Note:
1905: The textual name for each registered field must be unique.
1907: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1908: @*/
1909: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1910: {
1911: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1913: PetscFunctionBegin;
1914: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1915: {
1916: DMSwarmDataField gfield;
1918: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1919: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1920: }
1921: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1922: PetscFunctionReturn(PETSC_SUCCESS);
1923: }
1925: /*@
1926: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1928: Not Collective, No Fortran Support
1930: Input Parameters:
1931: + dm - a `DMSWARM`
1932: - fieldname - the textual name to identify this field
1934: Output Parameters:
1935: + blocksize - the number of each data type
1936: . type - the data type
1937: - data - pointer to raw array
1939: Level: beginner
1941: Notes:
1942: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1944: Fortran Note:
1945: Only works for `type` of `PETSC_SCALAR`
1947: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1948: @*/
1949: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1950: {
1951: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1952: DMSwarmDataField gfield;
1954: PetscFunctionBegin;
1956: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1957: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1958: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1959: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1960: if (blocksize) *blocksize = gfield->bs;
1961: if (type) *type = gfield->petsc_type;
1962: PetscFunctionReturn(PETSC_SUCCESS);
1963: }
1965: /*@
1966: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1968: Not Collective
1970: Input Parameters:
1971: + dm - a `DMSWARM`
1972: - fieldname - the textual name to identify this field
1974: Output Parameters:
1975: + blocksize - the number of each data type
1976: . type - the data type
1977: - data - pointer to raw array
1979: Level: beginner
1981: Notes:
1982: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1984: Fortran Note:
1985: Only works for `type` of `PETSC_SCALAR`
1987: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1988: @*/
1989: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1990: {
1991: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1992: DMSwarmDataField gfield;
1994: PetscFunctionBegin;
1996: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1997: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1998: if (data) *data = NULL;
1999: PetscFunctionReturn(PETSC_SUCCESS);
2000: }
2002: /*@
2003: DMSwarmGetFieldInfo - Return the block size and data type of a registered `DMSWARM` field without accessing its data.
2005: Not Collective
2007: Input Parameters:
2008: + dm - a `DMSWARM`
2009: - fieldname - the name of the registered field
2011: Output Parameters:
2012: + blocksize - the number of entries of `type` per particle, or `NULL`
2013: - type - the `PetscDataType` of a single entry, or `NULL`
2015: Level: intermediate
2017: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`, `DMSwarmRestoreField()`, `DMSwarmRegisterPetscDatatypeField()`
2018: @*/
2019: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
2020: {
2021: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2022: DMSwarmDataField gfield;
2024: PetscFunctionBegin;
2026: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
2027: if (blocksize) *blocksize = gfield->bs;
2028: if (type) *type = gfield->petsc_type;
2029: PetscFunctionReturn(PETSC_SUCCESS);
2030: }
2032: /*@
2033: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
2035: Not Collective
2037: Input Parameter:
2038: . dm - a `DMSWARM`
2040: Level: beginner
2042: Notes:
2043: The new point will have all fields initialized to zero.
2045: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
2046: @*/
2047: PetscErrorCode DMSwarmAddPoint(DM dm)
2048: {
2049: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2051: PetscFunctionBegin;
2052: if (!swarm->issetup) PetscCall(DMSetUp(dm));
2053: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
2054: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
2055: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
2056: PetscFunctionReturn(PETSC_SUCCESS);
2057: }
2059: /*@
2060: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
2062: Not Collective
2064: Input Parameters:
2065: + dm - a `DMSWARM`
2066: - npoints - the number of new points to add
2068: Level: beginner
2070: Notes:
2071: The new point will have all fields initialized to zero.
2073: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
2074: @*/
2075: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
2076: {
2077: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2078: PetscInt nlocal;
2080: PetscFunctionBegin;
2081: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
2082: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
2083: nlocal = PetscMax(nlocal, 0) + npoints;
2084: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
2085: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
2086: PetscFunctionReturn(PETSC_SUCCESS);
2087: }
2089: /*@
2090: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
2092: Not Collective
2094: Input Parameter:
2095: . dm - a `DMSWARM`
2097: Level: beginner
2099: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
2100: @*/
2101: PetscErrorCode DMSwarmRemovePoint(DM dm)
2102: {
2103: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2105: PetscFunctionBegin;
2106: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
2107: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
2108: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
2109: PetscFunctionReturn(PETSC_SUCCESS);
2110: }
2112: /*@
2113: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
2115: Not Collective
2117: Input Parameters:
2118: + dm - a `DMSWARM`
2119: - idx - index of point to remove
2121: Level: beginner
2123: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2124: @*/
2125: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
2126: {
2127: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2129: PetscFunctionBegin;
2130: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
2131: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
2132: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
2133: PetscFunctionReturn(PETSC_SUCCESS);
2134: }
2136: /*@
2137: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
2139: Not Collective
2141: Input Parameters:
2142: + dm - a `DMSWARM`
2143: . pi - the index of the point to copy
2144: - pj - the point index where the copy should be located
2146: Level: beginner
2148: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2149: @*/
2150: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
2151: {
2152: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2154: PetscFunctionBegin;
2155: if (!swarm->issetup) PetscCall(DMSetUp(dm));
2156: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
2157: PetscFunctionReturn(PETSC_SUCCESS);
2158: }
2160: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
2161: {
2162: PetscFunctionBegin;
2163: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
2164: PetscFunctionReturn(PETSC_SUCCESS);
2165: }
2167: /*@
2168: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
2170: Collective
2172: Input Parameters:
2173: + dm - the `DMSWARM`
2174: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
2176: Level: advanced
2178: Notes:
2179: The `DM` will be modified to accommodate received points.
2180: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
2181: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
2183: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
2184: @*/
2185: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
2186: {
2187: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2189: PetscFunctionBegin;
2190: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
2191: switch (swarm->migrate_type) {
2192: case DMSWARM_MIGRATE_BASIC:
2193: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
2194: break;
2195: case DMSWARM_MIGRATE_DMCELLNSCATTER:
2196: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
2197: break;
2198: case DMSWARM_MIGRATE_DMCELLEXACT:
2199: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
2200: case DMSWARM_MIGRATE_USER:
2201: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
2202: default:
2203: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
2204: }
2205: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
2206: PetscCall(DMClearGlobalVectors(dm));
2207: PetscFunctionReturn(PETSC_SUCCESS);
2208: }
2210: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
2212: /*
2213: DMSwarmCollectViewCreate
2215: * Applies a collection method and gathers point neighbour points into dm
2217: Notes:
2218: Users should call DMSwarmCollectViewDestroy() after
2219: they have finished computations associated with the collected points
2220: */
2222: /*@
2223: DMSwarmCollectViewCreate - Applies a collection method and gathers points
2224: in neighbour ranks into the `DMSWARM`
2226: Collective
2228: Input Parameter:
2229: . dm - the `DMSWARM`
2231: Level: advanced
2233: Notes:
2234: Users should call `DMSwarmCollectViewDestroy()` after
2235: they have finished computations associated with the collected points
2237: Different collect methods are supported. See `DMSwarmSetCollectType()`.
2239: Developer Note:
2240: Create and Destroy routines create new objects that can get destroyed, they do not change the state
2241: of the current object.
2243: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
2244: @*/
2245: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
2246: {
2247: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2248: PetscInt ng;
2250: PetscFunctionBegin;
2251: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
2252: PetscCall(DMSwarmGetLocalSize(dm, &ng));
2253: switch (swarm->collect_type) {
2254: case DMSWARM_COLLECT_BASIC:
2255: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
2256: break;
2257: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
2258: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
2259: case DMSWARM_COLLECT_GENERAL:
2260: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
2261: default:
2262: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
2263: }
2264: swarm->collect_view_active = PETSC_TRUE;
2265: swarm->collect_view_reset_nlocal = ng;
2266: PetscFunctionReturn(PETSC_SUCCESS);
2267: }
2269: /*@
2270: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
2272: Collective
2274: Input Parameters:
2275: . dm - the `DMSWARM`
2277: Notes:
2278: Users should call `DMSwarmCollectViewCreate()` before this function is called.
2280: Level: advanced
2282: Developer Note:
2283: Create and Destroy routines create new objects that can get destroyed, they do not change the state
2284: of the current object.
2286: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
2287: @*/
2288: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
2289: {
2290: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2292: PetscFunctionBegin;
2293: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
2294: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
2295: swarm->collect_view_active = PETSC_FALSE;
2296: PetscFunctionReturn(PETSC_SUCCESS);
2297: }
2299: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
2300: {
2301: PetscInt dim;
2303: PetscFunctionBegin;
2304: PetscCall(DMSwarmSetNumSpecies(dm, 1));
2305: PetscCall(DMGetDimension(dm, &dim));
2306: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2307: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2308: PetscFunctionReturn(PETSC_SUCCESS);
2309: }
2311: /*@
2312: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
2314: Collective
2316: Input Parameters:
2317: + dm - the `DMSWARM`
2318: - Npc - The number of particles per cell in the cell `DM`
2320: Level: intermediate
2322: Notes:
2323: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
2324: one particle is in each cell, it is placed at the centroid.
2326: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
2327: @*/
2328: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
2329: {
2330: DM cdm;
2331: DMSwarmCellDM celldm;
2332: PetscRandom rnd;
2333: DMPolytopeType ct;
2334: PetscBool simplex;
2335: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
2336: PetscInt dim, d, cStart, cEnd, c, p, Nfc;
2337: const char **coordFields;
2339: PetscFunctionBeginUser;
2340: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
2341: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
2342: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
2344: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2345: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2346: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2347: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2348: PetscCall(DMGetDimension(cdm, &dim));
2349: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
2350: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
2351: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
2353: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
2354: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
2355: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2356: for (c = cStart; c < cEnd; ++c) {
2357: if (Npc == 1) {
2358: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
2359: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
2360: } else {
2361: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
2362: for (p = 0; p < Npc; ++p) {
2363: const PetscInt n = c * Npc + p;
2364: PetscReal sum = 0.0, refcoords[3];
2366: for (d = 0; d < dim; ++d) {
2367: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
2368: sum += refcoords[d];
2369: }
2370: if (simplex && sum > 0.0)
2371: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
2372: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
2373: }
2374: }
2375: }
2376: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2377: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
2378: PetscCall(PetscRandomDestroy(&rnd));
2379: PetscFunctionReturn(PETSC_SUCCESS);
2380: }
2382: /*@
2383: DMSwarmGetType - Get particular flavor of `DMSWARM`
2385: Collective
2387: Input Parameter:
2388: . sw - the `DMSWARM`
2390: Output Parameter:
2391: . stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2393: Level: advanced
2395: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2396: @*/
2397: PetscErrorCode DMSwarmGetType(DM sw, DMSwarmType *stype)
2398: {
2399: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2401: PetscFunctionBegin;
2403: PetscAssertPointer(stype, 2);
2404: *stype = swarm->swarm_type;
2405: PetscFunctionReturn(PETSC_SUCCESS);
2406: }
2408: /*@
2409: DMSwarmSetType - Set particular flavor of `DMSWARM`
2411: Collective
2413: Input Parameters:
2414: + sw - the `DMSWARM`
2415: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2417: Level: advanced
2419: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2420: @*/
2421: PetscErrorCode DMSwarmSetType(DM sw, DMSwarmType stype)
2422: {
2423: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2425: PetscFunctionBegin;
2427: swarm->swarm_type = stype;
2428: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(sw));
2429: PetscFunctionReturn(PETSC_SUCCESS);
2430: }
2432: static PetscErrorCode DMSwarmCreateRemapDM_Private(DM sw, DM *rdm)
2433: {
2434: PetscFE fe;
2435: DMPolytopeType ct;
2436: PetscInt dim, cStart;
2437: const char *prefix = "remap_";
2439: PetscFunctionBegin;
2440: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), rdm));
2441: PetscCall(DMSetType(*rdm, DMPLEX));
2442: PetscCall(DMPlexSetOptionsPrefix(*rdm, prefix));
2443: PetscCall(DMSetFromOptions(*rdm));
2444: PetscCall(PetscObjectSetName((PetscObject)*rdm, "remap"));
2445: PetscCall(DMViewFromOptions(*rdm, NULL, "-dm_view"));
2447: PetscCall(DMGetDimension(*rdm, &dim));
2448: PetscCall(DMPlexGetHeightStratum(*rdm, 0, &cStart, NULL));
2449: PetscCall(DMPlexGetCellType(*rdm, cStart, &ct));
2450: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
2451: PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
2452: PetscCall(DMSetField(*rdm, 0, NULL, (PetscObject)fe));
2453: PetscCall(DMCreateDS(*rdm));
2454: PetscCall(PetscFEDestroy(&fe));
2455: PetscFunctionReturn(PETSC_SUCCESS);
2456: }
2458: static PetscErrorCode DMSetup_Swarm(DM sw)
2459: {
2460: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2462: PetscFunctionBegin;
2463: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
2464: swarm->issetup = PETSC_TRUE;
2466: if (swarm->remap_type != DMSWARM_REMAP_NONE) {
2467: DMSwarmCellDM celldm;
2469: PetscCall(DMSwarmGetCellDMByName(sw, "remap", &celldm));
2470: if (!celldm) {
2471: DM rdm;
2472: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
2473: const char *vfieldnames[1] = {"w_q"};
2475: PetscCall(DMSwarmCreateRemapDM_Private(sw, &rdm));
2476: PetscCall(DMSwarmCellDMCreate(rdm, 1, vfieldnames, 2, fieldnames, &celldm));
2477: PetscCall(DMSwarmAddCellDM(sw, celldm));
2478: PetscCall(DMSwarmCellDMDestroy(&celldm));
2479: PetscCall(DMDestroy(&rdm));
2480: }
2481: }
2483: if (swarm->swarm_type == DMSWARM_PIC) {
2484: DMSwarmCellDM celldm;
2486: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2487: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "No active cell DM. DMSWARM_PIC requires you call DMSwarmSetCellDM() or DMSwarmAddCellDM()");
2488: if (celldm->dm->ops->locatepointssubdomain) {
2489: /* check methods exists for exact ownership identificiation */
2490: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
2491: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
2492: } else {
2493: /* check methods exist for point location AND rank neighbor identification */
2494: PetscCheck(celldm->dm->ops->locatepoints, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
2495: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
2497: PetscCheck(celldm->dm->ops->getneighbors, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
2498: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
2500: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
2501: }
2502: }
2504: PetscCall(DMSwarmFinalizeFieldRegister(sw));
2506: /* check some fields were registered */
2507: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
2508: PetscFunctionReturn(PETSC_SUCCESS);
2509: }
2511: static PetscErrorCode DMDestroy_Swarm(DM dm)
2512: {
2513: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2515: PetscFunctionBegin;
2516: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2517: PetscCall(PetscObjectListDestroy(&swarm->cellDMs));
2518: PetscCall(PetscFree(swarm->activeCellDM));
2519: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
2520: PetscCall(PetscFree(swarm));
2521: PetscFunctionReturn(PETSC_SUCCESS);
2522: }
2524: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
2525: {
2526: DM cdm;
2527: DMSwarmCellDM celldm;
2528: PetscDraw draw;
2529: PetscReal *coords, oldPause, radius = 0.01;
2530: PetscInt Np, p, bs, Nfc;
2531: const char **coordFields;
2533: PetscFunctionBegin;
2534: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
2535: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2536: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2537: PetscCall(PetscDrawGetPause(draw, &oldPause));
2538: PetscCall(PetscDrawSetPause(draw, 0.0));
2539: PetscCall(DMView(cdm, viewer));
2540: PetscCall(PetscDrawSetPause(draw, oldPause));
2542: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2543: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2544: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2545: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2546: PetscCall(DMSwarmGetField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2547: for (p = 0; p < Np; ++p) {
2548: const PetscInt i = p * bs;
2550: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
2551: }
2552: PetscCall(DMSwarmRestoreField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2553: PetscCall(PetscDrawFlush(draw));
2554: PetscCall(PetscDrawPause(draw));
2555: PetscCall(PetscDrawSave(draw));
2556: PetscFunctionReturn(PETSC_SUCCESS);
2557: }
2559: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
2560: {
2561: PetscViewerFormat format;
2562: PetscDataType wtype = PETSC_DATATYPE_UNKNOWN;
2563: PetscInt *sizes;
2564: PetscInt dim, Np, wbs = 0, maxSize = 17;
2565: MPI_Comm comm;
2566: PetscMPIInt rank, size;
2567: const char *name, *cellid;
2569: PetscFunctionBegin;
2570: PetscCall(PetscViewerGetFormat(viewer, &format));
2571: PetscCall(DMGetDimension(dm, &dim));
2572: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2573: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2574: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2575: PetscCallMPI(MPI_Comm_size(comm, &size));
2576: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
2577: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
2578: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
2579: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
2580: else PetscCall(PetscCalloc1(3, &sizes));
2581: if (size < maxSize) {
2582: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
2583: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
2584: for (PetscInt p = 0; p < size; ++p) {
2585: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
2586: }
2587: } else {
2588: PetscInt locMinMax[2] = {Np, Np};
2590: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
2591: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
2592: }
2593: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2594: PetscCall(PetscFree(sizes));
2595: if (format == PETSC_VIEWER_ASCII_INFO) {
2596: DM_Swarm *sw = (DM_Swarm *)dm->data;
2597: DMSwarmCellDM celldm;
2598: PetscInt *cell;
2599: PetscBool hasWeight;
2600: const char *fname = "w_q";
2602: PetscCall(DMSwarmDataFieldStringInList(fname, sw->db->nfields, (const DMSwarmDataField *)sw->db->field, &hasWeight));
2603: if (hasWeight) {
2604: PetscCall(DMSwarmGetFieldInfo(dm, fname, &wbs, &wtype));
2605: PetscCheck(wtype == PETSC_REAL || wtype == PETSC_SCALAR, comm, PETSC_ERR_SUP, "Weight field %s must have type PETSC_REAL or PETSC_SCALAR, not %s", fname, PetscDataTypes[wtype]);
2606: PetscCheck(wbs == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must have block size 1, not %" PetscInt_FMT, fname, wbs);
2607: }
2608: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
2609: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2610: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2611: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
2612: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&cell));
2613: if (hasWeight) {
2614: PetscReal **coords;
2615: void *weight;
2616: PetscInt Ncf, *bsC;
2617: const char **coordNames;
2619: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncf, &coordNames));
2620: PetscCall(PetscMalloc2(Ncf, &coords, Ncf, &bsC));
2621: for (PetscInt n = 0; n < Ncf; ++n) PetscCall(DMSwarmGetField(dm, coordNames[n], &bsC[n], NULL, (void **)&coords[n]));
2622: PetscCall(DMSwarmGetField(dm, fname, NULL, NULL, &weight));
2623: for (PetscInt p = 0; p < Np; ++p) {
2624: PetscReal wp;
2626: if (wtype == PETSC_REAL) wp = ((PetscReal *)weight)[p];
2627: else wp = PetscRealPart(((PetscScalar *)weight)[p]);
2628: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%" PetscInt_FMT ": %" PetscInt_FMT " wt: %g x: (", p, cell[p], (double)wp));
2629: for (PetscInt n = 0; n < Ncf; ++n) {
2630: if (n > 0) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, ", "));
2631: for (PetscInt d = 0; d < bsC[n]; ++d) {
2632: if (d > 0) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, ", "));
2633: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)coords[n][p * bsC[n] + d]));
2634: }
2635: }
2636: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, ")\n"));
2637: }
2638: PetscCall(DMSwarmRestoreField(dm, fname, NULL, NULL, &weight));
2639: for (PetscInt n = 0; n < Ncf; ++n) PetscCall(DMSwarmRestoreField(dm, coordNames[n], &bsC[n], NULL, (void **)&coords[n]));
2640: PetscCall(PetscFree2(coords, bsC));
2641: } else {
2642: for (PetscInt p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%" PetscInt_FMT ": %" PetscInt_FMT "\n", p, cell[p]));
2643: }
2644: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&cell));
2645: PetscCall(PetscViewerFlush(viewer));
2646: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2647: }
2648: PetscFunctionReturn(PETSC_SUCCESS);
2649: }
2651: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
2652: {
2653: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2654: PetscBool isascii, ibinary, isvtk, isdraw, ispython;
2655: #if PetscDefined(HAVE_HDF5)
2656: PetscBool ishdf5;
2657: #endif
2659: PetscFunctionBegin;
2662: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2663: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
2664: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
2665: #if PetscDefined(HAVE_HDF5)
2666: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2667: #endif
2668: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2669: PetscCall(PetscObjectHasFunction((PetscObject)viewer, "PetscViewerPythonViewObject_C", &ispython));
2670: if (isascii) {
2671: PetscViewerFormat format;
2673: PetscCall(PetscViewerGetFormat(viewer, &format));
2674: switch (format) {
2675: case PETSC_VIEWER_ASCII_INFO_DETAIL:
2676: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
2677: break;
2678: default:
2679: PetscCall(DMView_Swarm_Ascii(dm, viewer));
2680: }
2681: } else {
2682: #if PetscDefined(HAVE_HDF5)
2683: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
2684: #endif
2685: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
2686: if (ispython) PetscCall(PetscViewerPythonViewObject(viewer, (PetscObject)dm));
2687: }
2688: PetscFunctionReturn(PETSC_SUCCESS);
2689: }
2691: /*@
2692: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
2693: 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.
2695: Noncollective
2697: Input Parameters:
2698: + sw - the `DMSWARM`
2699: . cellID - the integer id of the cell to be extracted and filtered
2700: - cellswarm - The `DMSWARM` to receive the cell
2702: Level: beginner
2704: Notes:
2705: This presently only supports `DMSWARM_PIC` type
2707: Should be restored with `DMSwarmRestoreCellSwarm()`
2709: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
2711: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
2712: @*/
2713: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2714: {
2715: DM_Swarm *original = (DM_Swarm *)sw->data;
2716: DMLabel label;
2717: DM dmc, subdmc;
2718: PetscInt *pids, particles, dim;
2719: const char *name;
2721: PetscFunctionBegin;
2722: /* Configure new swarm */
2723: PetscCall(DMSetType(cellswarm, DMSWARM));
2724: PetscCall(DMGetDimension(sw, &dim));
2725: PetscCall(DMSetDimension(cellswarm, dim));
2726: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
2727: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
2728: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
2729: PetscCall(DMSwarmSortGetAccess(sw));
2730: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
2731: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2732: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
2733: PetscCall(DMSwarmSortRestoreAccess(sw));
2734: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2735: PetscCall(DMSwarmGetCellDM(sw, &dmc));
2736: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
2737: PetscCall(DMAddLabel(dmc, label));
2738: PetscCall(DMLabelSetValue(label, cellID, 1));
2739: PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, PetscObjectComm((PetscObject)dmc), NULL, &subdmc));
2740: PetscCall(PetscObjectGetName((PetscObject)dmc, &name));
2741: PetscCall(PetscObjectSetName((PetscObject)subdmc, name));
2742: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
2743: PetscCall(DMLabelDestroy(&label));
2744: PetscFunctionReturn(PETSC_SUCCESS);
2745: }
2747: /*@
2748: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
2750: Noncollective
2752: Input Parameters:
2753: + sw - the parent `DMSWARM`
2754: . cellID - the integer id of the cell to be copied back into the parent swarm
2755: - cellswarm - the cell swarm object
2757: Level: beginner
2759: Note:
2760: This only supports `DMSWARM_PIC` types of `DMSWARM`s
2762: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
2763: @*/
2764: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2765: {
2766: DM dmc;
2767: PetscInt *pids, particles, p;
2769: PetscFunctionBegin;
2770: PetscCall(DMSwarmSortGetAccess(sw));
2771: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2772: PetscCall(DMSwarmSortRestoreAccess(sw));
2773: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
2774: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
2775: /* Free memory, destroy cell dm */
2776: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
2777: PetscCall(DMDestroy(&dmc));
2778: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2779: PetscFunctionReturn(PETSC_SUCCESS);
2780: }
2782: /*@
2783: DMSwarmComputeMoments - Compute the first three particle moments for a given field
2785: Noncollective
2787: Input Parameters:
2788: + sw - the `DMSWARM`
2789: . coordinate - the coordinate field name
2790: - weight - the weight field name
2792: Output Parameter:
2793: . moments - the field moments
2795: Level: intermediate
2797: Notes:
2798: The `moments` array should be of length bs + 2, where bs is the block size of the coordinate field.
2800: The weight field must have blocksize 1 and data type `PETSC_REAL` or `PETSC_SCALAR`. For complex scalars,
2801: only the real part of each weight is used.
2803: .seealso: `DM`, `DMSWARM`, `DMPlexComputeMoments()`
2804: @*/
2805: PetscErrorCode DMSwarmComputeMoments(DM sw, const char coordinate[], const char weight[], PetscReal moments[])
2806: {
2807: const PetscReal *coords;
2808: void *weights;
2809: PetscDataType dtc, dtw;
2810: PetscInt bsc, bsw, Np;
2811: MPI_Comm comm;
2813: PetscFunctionBegin;
2815: PetscAssertPointer(coordinate, 2);
2816: PetscAssertPointer(weight, 3);
2817: PetscAssertPointer(moments, 4);
2818: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
2819: PetscCall(DMSwarmGetFieldInfo(sw, coordinate, &bsc, &dtc));
2820: PetscCall(DMSwarmGetFieldInfo(sw, weight, &bsw, &dtw));
2821: PetscCheck(dtc == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s must be real, not %s", coordinate, PetscDataTypes[dtc]);
2822: PetscCheck(dtw == PETSC_REAL || dtw == PETSC_SCALAR, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must have type PETSC_REAL or PETSC_SCALAR, not %s", weight, PetscDataTypes[dtw]);
2823: PetscCheck(bsw == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be a scalar, not blocksize %" PetscInt_FMT, weight, bsw);
2824: PetscCall(DMSwarmGetField(sw, coordinate, NULL, NULL, (void **)&coords));
2825: PetscCall(DMSwarmGetField(sw, weight, NULL, NULL, &weights));
2826: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2827: PetscCall(PetscArrayzero(moments, bsc + 2));
2828: for (PetscInt p = 0; p < Np; ++p) {
2829: const PetscReal *c = &coords[p * bsc];
2830: PetscReal wp;
2832: if (dtw == PETSC_REAL) wp = ((PetscReal *)weights)[p];
2833: else wp = PetscRealPart(((PetscScalar *)weights)[p]);
2834: moments[0] += wp;
2835: for (PetscInt d = 0; d < bsc; ++d) {
2836: moments[d + 1] += wp * c[d];
2837: moments[bsc + 1] += wp * PetscSqr(c[d]);
2838: }
2839: }
2840: PetscCall(DMSwarmRestoreField(sw, coordinate, NULL, NULL, (void **)&coords));
2841: PetscCall(DMSwarmRestoreField(sw, weight, NULL, NULL, &weights));
2842: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, moments, bsc + 2, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
2843: PetscFunctionReturn(PETSC_SUCCESS);
2844: }
2846: static PetscErrorCode DMSetFromOptions_Swarm(DM dm, PetscOptionItems PetscOptionsObject)
2847: {
2848: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2850: PetscFunctionBegin;
2851: PetscOptionsHeadBegin(PetscOptionsObject, "DMSwarm Options");
2852: PetscCall(PetscOptionsEnum("-dm_swarm_remap_type", "Remap algorithm", "DMSwarmSetRemapType", DMSwarmRemapTypeNames, (PetscEnum)swarm->remap_type, (PetscEnum *)&swarm->remap_type, NULL));
2853: PetscOptionsHeadEnd();
2854: PetscFunctionReturn(PETSC_SUCCESS);
2855: }
2857: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
2859: static PetscErrorCode DMInitialize_Swarm(DM sw)
2860: {
2861: PetscFunctionBegin;
2862: sw->ops->view = DMView_Swarm;
2863: sw->ops->load = NULL;
2864: sw->ops->setfromoptions = DMSetFromOptions_Swarm;
2865: sw->ops->clone = DMClone_Swarm;
2866: sw->ops->setup = DMSetup_Swarm;
2867: sw->ops->createlocalsection = NULL;
2868: sw->ops->createsectionpermutation = NULL;
2869: sw->ops->createdefaultconstraints = NULL;
2870: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
2871: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
2872: sw->ops->getlocaltoglobalmapping = NULL;
2873: sw->ops->createfieldis = NULL;
2874: sw->ops->createcoordinatedm = NULL;
2875: sw->ops->createcellcoordinatedm = NULL;
2876: sw->ops->getcoloring = NULL;
2877: sw->ops->creatematrix = DMCreateMatrix_Swarm;
2878: sw->ops->createinterpolation = NULL;
2879: sw->ops->createinjection = NULL;
2880: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
2881: sw->ops->creategradientmatrix = DMCreateGradientMatrix_Swarm;
2882: sw->ops->refine = NULL;
2883: sw->ops->coarsen = NULL;
2884: sw->ops->refinehierarchy = NULL;
2885: sw->ops->coarsenhierarchy = NULL;
2886: sw->ops->globaltolocalbegin = DMGlobalToLocalBegin_Swarm;
2887: sw->ops->globaltolocalend = DMGlobalToLocalEnd_Swarm;
2888: sw->ops->localtoglobalbegin = DMLocalToGlobalBegin_Swarm;
2889: sw->ops->localtoglobalend = DMLocalToGlobalEnd_Swarm;
2890: sw->ops->destroy = DMDestroy_Swarm;
2891: sw->ops->createsubdm = NULL;
2892: sw->ops->getdimpoints = NULL;
2893: sw->ops->locatepoints = NULL;
2894: sw->ops->projectfieldlocal = DMProjectFieldLocal_Swarm;
2895: PetscFunctionReturn(PETSC_SUCCESS);
2896: }
2898: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
2899: {
2900: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2902: PetscFunctionBegin;
2903: swarm->refct++;
2904: (*newdm)->data = swarm;
2905: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
2906: PetscCall(DMInitialize_Swarm(*newdm));
2907: (*newdm)->dim = dm->dim;
2908: PetscFunctionReturn(PETSC_SUCCESS);
2909: }
2911: /*MC
2912: DMSWARM = "swarm" - A `DM` object for particle methods, such as particle-in-cell (PIC), in which the underlying
2913: data is both (i) dynamic in length, (ii) and of arbitrary data type.
2915: Level: intermediate
2917: Notes:
2918: User data can be represented by `DMSWARM` through a registering "fields" which are to be stored on particles.
2919: To register a field, the user must provide;
2920: (a) a unique name;
2921: (b) the data type (or size in bytes);
2922: (c) the block size of the data.
2924: For example, suppose the application requires a unique id, energy, momentum and density to be stored
2925: on a set of particles. Then the following code could be used
2926: .vb
2927: DMSwarmInitializeFieldRegister(dm)
2928: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
2929: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
2930: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
2931: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
2932: DMSwarmFinalizeFieldRegister(dm)
2933: .ve
2935: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
2936: The only restriction imposed by `DMSWARM` is that all fields contain the same number of particles.
2938: To support particle methods, "migration" techniques are provided. These methods migrate data
2939: between ranks.
2941: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
2942: As a `DMSWARM` may internally define and store values of different data types,
2943: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
2944: fields should be used to define a `Vec` object via `DMSwarmVectorDefineField()`
2945: The specified field can be changed at any time - thereby permitting vectors
2946: compatible with different fields to be created.
2948: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via `DMSwarmCreateGlobalVectorFromField()`
2949: Here the data defining the field in the `DMSWARM` is shared with a `Vec`.
2950: This is inherently unsafe if you alter the size of the field at any time between
2951: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
2952: If the local size of the `DMSWARM` does not match the local size of the global vector
2953: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
2955: Additional high-level support is provided for Particle-In-Cell methods. Refer to `DMSwarmSetType()`.
2957: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`, `DMSwarmSetType()`, `DMSwarmType`, `DMSwarmCreateGlobalVectorFromField()`,
2958: `DMCreateGlobalVector()`, `DMCreateLocalVector()`
2959: M*/
2961: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
2962: {
2963: DM_Swarm *swarm;
2965: PetscFunctionBegin;
2967: PetscCall(PetscNew(&swarm));
2968: dm->data = swarm;
2969: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
2970: PetscCall(DMSwarmInitializeFieldRegister(dm));
2971: dm->dim = 0;
2972: swarm->refct = 1;
2973: swarm->issetup = PETSC_FALSE;
2974: swarm->swarm_type = DMSWARM_BASIC;
2975: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
2976: swarm->collect_type = DMSWARM_COLLECT_BASIC;
2977: swarm->migrate_error_on_missing_point = PETSC_FALSE;
2978: swarm->collect_view_active = PETSC_FALSE;
2979: swarm->collect_view_reset_nlocal = -1;
2980: PetscCall(DMInitialize_Swarm(dm));
2981: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
2982: PetscFunctionReturn(PETSC_SUCCESS);
2983: }
2985: /* Replace dm with the contents of ndm, and then destroy ndm
2986: - Share the DM_Swarm structure
2987: */
2988: /*@
2989: DMSwarmReplace - Replace the internal state of a `DMSWARM` with that of another `DMSWARM`, sharing the underlying particle data and destroying the source `DM`.
2991: Collective
2993: Input Parameters:
2994: + dm - the destination `DMSWARM`, whose current contents are discarded
2995: - ndm - pointer to the source `DMSWARM`; destroyed and set to `NULL` on return
2997: Level: developer
2999: Note:
3000: The dimension, periodicity, name, and shared reference to the underlying particle bucket are transferred from the source to the destination.
3002: .seealso: `DM`, `DMSWARM`, `DMSwarmDuplicate()`, `DMDestroy()`
3003: @*/
3004: PetscErrorCode DMSwarmReplace(DM dm, DM *ndm)
3005: {
3006: DM dmNew = *ndm;
3007: const PetscReal *maxCell, *Lstart, *L;
3008: PetscInt dim;
3010: PetscFunctionBegin;
3011: if (dm == dmNew) {
3012: PetscCall(DMDestroy(ndm));
3013: PetscFunctionReturn(PETSC_SUCCESS);
3014: }
3015: dm->setupcalled = dmNew->setupcalled;
3016: if (!dm->hdr.name) {
3017: const char *name;
3019: PetscCall(PetscObjectGetName((PetscObject)*ndm, &name));
3020: PetscCall(PetscObjectSetName((PetscObject)dm, name));
3021: }
3022: PetscCall(DMGetDimension(dmNew, &dim));
3023: PetscCall(DMSetDimension(dm, dim));
3024: PetscCall(DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L));
3025: PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
3026: PetscCall(DMDestroy_Swarm(dm));
3027: PetscCall(DMInitialize_Swarm(dm));
3028: dm->data = dmNew->data;
3029: ((DM_Swarm *)dmNew->data)->refct++;
3030: PetscCall(DMDestroy(ndm));
3031: PetscFunctionReturn(PETSC_SUCCESS);
3032: }
3034: /*@
3035: DMSwarmDuplicate - Creates a new `DMSWARM` with the same fields and cell `DM`s but no particles
3037: Collective
3039: Input Parameter:
3040: . sw - the `DMSWARM`
3042: Output Parameter:
3043: . nsw - the new `DMSWARM`
3045: Level: beginner
3047: .seealso: `DM`, `DMSWARM`, `DMSwarmCreate()`, `DMClone()`
3048: @*/
3049: PetscErrorCode DMSwarmDuplicate(DM sw, DM *nsw)
3050: {
3051: DM_Swarm *swarm = (DM_Swarm *)sw->data;
3052: DMSwarmDataField *fields;
3053: DMSwarmCellDM celldm, ncelldm;
3054: DMSwarmType stype;
3055: const char *name, **celldmnames;
3056: void *ctx;
3057: PetscInt dim, Nf, Ndm;
3058: PetscBool flg;
3060: PetscFunctionBegin;
3061: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), nsw));
3062: PetscCall(DMSetType(*nsw, DMSWARM));
3063: PetscCall(PetscObjectGetName((PetscObject)sw, &name));
3064: PetscCall(PetscObjectSetName((PetscObject)*nsw, name));
3065: PetscCall(DMGetDimension(sw, &dim));
3066: PetscCall(DMSetDimension(*nsw, dim));
3067: PetscCall(DMSwarmGetType(sw, &stype));
3068: PetscCall(DMSwarmSetType(*nsw, stype));
3069: PetscCall(DMGetApplicationContext(sw, &ctx));
3070: PetscCall(DMSetApplicationContext(*nsw, ctx));
3072: PetscCall(DMSwarmDataBucketGetDMSwarmDataFields(swarm->db, &Nf, &fields));
3073: for (PetscInt f = 0; f < Nf; ++f) {
3074: PetscCall(DMSwarmDataFieldStringInList(fields[f]->name, ((DM_Swarm *)(*nsw)->data)->db->nfields, (const DMSwarmDataField *)((DM_Swarm *)(*nsw)->data)->db->field, &flg));
3075: if (!flg) PetscCall(DMSwarmRegisterPetscDatatypeField(*nsw, fields[f]->name, fields[f]->bs, fields[f]->petsc_type));
3076: }
3078: PetscCall(DMSwarmGetCellDMNames(sw, &Ndm, &celldmnames));
3079: for (PetscInt c = 0; c < Ndm; ++c) {
3080: DM dm;
3081: PetscInt Ncf;
3082: const char **coordfields, **fields;
3084: PetscCall(DMSwarmGetCellDMByName(sw, celldmnames[c], &celldm));
3085: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
3086: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncf, &coordfields));
3087: PetscCall(DMSwarmCellDMGetFields(celldm, &Nf, &fields));
3088: PetscCall(DMSwarmCellDMCreate(dm, Nf, fields, Ncf, coordfields, &ncelldm));
3089: PetscCall(DMSwarmAddCellDM(*nsw, ncelldm));
3090: PetscCall(DMSwarmCellDMDestroy(&ncelldm));
3091: }
3092: PetscCall(PetscFree(celldmnames));
3094: PetscCall(DMSetFromOptions(*nsw));
3095: PetscCall(DMSetUp(*nsw));
3096: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
3097: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
3098: PetscCall(DMSwarmSetCellDMActive(*nsw, name));
3099: PetscFunctionReturn(PETSC_SUCCESS);
3100: }
3102: PetscErrorCode DMLocalToGlobalBegin_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
3103: {
3104: PetscFunctionBegin;
3105: PetscFunctionReturn(PETSC_SUCCESS);
3106: }
3108: PetscErrorCode DMLocalToGlobalEnd_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
3109: {
3110: PetscFunctionBegin;
3111: switch (mode) {
3112: case INSERT_VALUES:
3113: PetscCall(VecCopy(l, g));
3114: break;
3115: case ADD_VALUES:
3116: PetscCall(VecAXPY(g, 1., l));
3117: break;
3118: default:
3119: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mode not supported: %d", mode);
3120: }
3121: PetscFunctionReturn(PETSC_SUCCESS);
3122: }
3124: PetscErrorCode DMGlobalToLocalBegin_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
3125: {
3126: PetscFunctionBegin;
3127: PetscFunctionReturn(PETSC_SUCCESS);
3128: }
3130: PetscErrorCode DMGlobalToLocalEnd_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
3131: {
3132: PetscFunctionBegin;
3133: PetscCall(DMLocalToGlobalEnd_Swarm(dm, g, mode, l));
3134: PetscFunctionReturn(PETSC_SUCCESS);
3135: }