Actual source code: swarm.c
1: #include "petscdmswarm.h"
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petsc/private/hashsetij.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petscviewer.h>
6: #include <petscdraw.h>
7: #include <petscdmplex.h>
8: #include <petscblaslapack.h>
9: #include "../src/dm/impls/swarm/data_bucket.h"
10: #include <petscdmlabel.h>
11: #include <petscsection.h>
13: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
14: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
15: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
17: const char *DMSwarmTypeNames[] = {"basic", "pic", NULL};
18: const char *DMSwarmMigrateTypeNames[] = {"basic", "dmcellnscatter", "dmcellexact", "user", NULL};
19: const char *DMSwarmCollectTypeNames[] = {"basic", "boundingbox", "general", "user", NULL};
20: const char *DMSwarmRemapTypeNames[] = {"none", "pfak", "colella", "DMSwarmRemapType", "DMSWARM_REMAP_", NULL};
21: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
23: const char DMSwarmField_pid[] = "DMSwarm_pid";
24: const char DMSwarmField_rank[] = "DMSwarm_rank";
25: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
27: PetscInt SwarmDataFieldId = -1;
29: #if PetscDefined(HAVE_HDF5)
30: #include <petscviewerhdf5.h>
32: static PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
33: {
34: DM dm;
35: PetscReal seqval;
36: PetscInt seqnum, bs;
37: PetscBool isseq, ists;
39: PetscFunctionBegin;
40: PetscCall(VecGetDM(v, &dm));
41: PetscCall(VecGetBlockSize(v, &bs));
42: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particle_fields"));
43: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
44: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
45: PetscCall(PetscViewerHDF5IsTimestepping(viewer, &ists));
46: if (ists) PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
47: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
48: PetscCall(VecViewNative(v, viewer));
49: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
50: PetscCall(PetscViewerHDF5PopGroup(viewer));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
55: {
56: DMSwarmCellDM celldm;
57: Vec coordinates;
58: PetscInt Np, Nfc;
59: PetscBool isseq;
60: const char **coordFields;
62: PetscFunctionBegin;
63: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
64: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
65: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
66: PetscCall(DMSwarmGetSize(dm, &Np));
67: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, coordFields[0], &coordinates));
68: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
69: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
70: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
71: PetscCall(VecViewNative(coordinates, viewer));
72: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
73: PetscCall(PetscViewerHDF5PopGroup(viewer));
74: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, coordFields[0], &coordinates));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
77: #endif
79: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
80: {
81: DM dm;
82: #if PetscDefined(HAVE_HDF5)
83: PetscBool ishdf5;
84: #endif
86: PetscFunctionBegin;
87: PetscCall(VecGetDM(v, &dm));
88: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
89: #if PetscDefined(HAVE_HDF5)
90: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
91: if (ishdf5) {
92: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
95: #endif
96: PetscCall(VecViewNative(v, viewer));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: /*@
101: DMSwarmVectorGetField - Gets the fields from which to define a `Vec` object
102: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
104: Not collective
106: Input Parameter:
107: . sw - a `DMSWARM`
109: Output Parameters:
110: + Nf - the number of fields
111: - fieldnames - the textual name given to each registered field, or NULL if it has not been set
113: Level: beginner
115: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
116: @*/
117: PetscErrorCode DMSwarmVectorGetField(DM sw, PetscInt *Nf, const char **fieldnames[])
118: {
119: DMSwarmCellDM celldm;
121: PetscFunctionBegin;
123: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
124: PetscCall(DMSwarmCellDMGetFields(celldm, Nf, fieldnames));
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*@
129: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
130: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
132: Collective
134: Input Parameters:
135: + dm - a `DMSWARM`
136: - fieldname - the textual name given to each registered field
138: Level: beginner
140: Notes:
141: The field with name `fieldname` must have data type `PETSC_REAL` or `PETSC_SCALAR`.
143: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
144: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
146: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineFields()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
147: @*/
148: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
149: {
150: PetscFunctionBegin;
151: PetscCall(DMSwarmVectorDefineFields(dm, 1, &fieldname));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: /*@
156: DMSwarmVectorDefineFields - Sets the fields from which to define a `Vec` object
157: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
159: Collective, No Fortran support
161: Input Parameters:
162: + sw - a `DMSWARM`
163: . Nf - the number of fields
164: - fieldnames - the textual name given to each registered field
166: Level: beginner
168: Notes:
169: Each field with name in `fieldnames` must have data type `PETSC_REAL` or `PETSC_SCALAR`.
171: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
172: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
174: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
175: @*/
176: PetscErrorCode DMSwarmVectorDefineFields(DM sw, PetscInt Nf, const char *fieldnames[])
177: {
178: DM_Swarm *swarm = (DM_Swarm *)sw->data;
179: DMSwarmCellDM celldm;
181: PetscFunctionBegin;
183: if (fieldnames) PetscAssertPointer(fieldnames, 3);
184: if (!swarm->issetup) PetscCall(DMSetUp(sw));
185: PetscCheck(Nf >= 0, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields must be non-negative, not %" PetscInt_FMT, Nf);
186: for (PetscInt f = 0; f < Nf; ++f) {
187: PetscDataType type;
189: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], NULL, &type));
190: 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]);
191: }
192: // Create a dummy cell DM if none has been specified (I think we should not support this mode)
193: if (!swarm->activeCellDM) {
194: DM dm;
195: DMSwarmCellDM celldm;
197: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), &dm));
198: PetscCall(DMSetType(dm, DMSHELL));
199: PetscCall(PetscObjectSetName((PetscObject)dm, "dummy"));
200: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 0, NULL, &celldm));
201: PetscCall(DMDestroy(&dm));
202: PetscCall(DMSwarmAddCellDM(sw, celldm));
203: PetscCall(DMSwarmCellDMDestroy(&celldm));
204: PetscCall(DMSwarmSetCellDMActive(sw, "dummy"));
205: }
206: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
207: for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscFree(celldm->dmFields[f]));
208: PetscCall(PetscFree(celldm->dmFields));
210: celldm->Nf = Nf;
211: PetscCall(PetscMalloc1(Nf, &celldm->dmFields));
212: for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscStrallocpy(fieldnames[f], (char **)&celldm->dmFields[f]));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /* requires DMSwarmDefineFieldVector has been called */
217: static PetscErrorCode DMCreateGlobalVector_Swarm(DM sw, Vec *vec)
218: {
219: DM_Swarm *swarm = (DM_Swarm *)sw->data;
220: DMSwarmCellDM celldm;
221: Vec x;
222: char name[PETSC_MAX_PATH_LEN];
223: PetscInt bs = 0, n;
225: PetscFunctionBegin;
226: if (!swarm->issetup) PetscCall(DMSetUp(sw));
227: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
228: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
229: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
231: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
232: for (PetscInt f = 0; f < celldm->Nf; ++f) {
233: PetscInt fbs;
234: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
235: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
236: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
237: bs += fbs;
238: }
239: PetscCall(VecCreate(PetscObjectComm((PetscObject)sw), &x));
240: PetscCall(PetscObjectSetName((PetscObject)x, name));
241: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
242: PetscCall(VecSetBlockSize(x, bs));
243: PetscCall(VecSetDM(x, sw));
244: PetscCall(VecSetFromOptions(x));
245: PetscCall(VecSetOperation(x, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
246: *vec = x;
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: /* requires DMSwarmDefineFieldVector has been called */
251: static PetscErrorCode DMCreateLocalVector_Swarm(DM sw, Vec *vec)
252: {
253: DM_Swarm *swarm = (DM_Swarm *)sw->data;
254: DMSwarmCellDM celldm;
255: Vec x;
256: char name[PETSC_MAX_PATH_LEN];
257: PetscInt bs = 0, n;
259: PetscFunctionBegin;
260: if (!swarm->issetup) PetscCall(DMSetUp(sw));
261: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
262: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
263: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
265: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
266: for (PetscInt f = 0; f < celldm->Nf; ++f) {
267: PetscInt fbs;
268: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
269: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
270: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
271: bs += fbs;
272: }
273: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
274: PetscCall(PetscObjectSetName((PetscObject)x, name));
275: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
276: PetscCall(VecSetBlockSize(x, bs));
277: PetscCall(VecSetDM(x, sw));
278: PetscCall(VecSetFromOptions(x));
279: *vec = x;
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
284: {
285: DM_Swarm *swarm = (DM_Swarm *)dm->data;
286: DMSwarmDataField gfield;
287: PetscInt bs, nlocal, fid = -1, cfid = -2;
288: PetscBool flg;
290: PetscFunctionBegin;
291: /* check vector is an inplace array */
292: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
293: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
294: (void)flg; /* avoid compiler warning */
295: 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);
296: PetscCall(VecGetLocalSize(*vec, &nlocal));
297: PetscCall(VecGetBlockSize(*vec, &bs));
298: 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");
299: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
300: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
301: PetscCall(VecResetArray(*vec));
302: PetscCall(VecDestroy(vec));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
307: {
308: DM_Swarm *swarm = (DM_Swarm *)dm->data;
309: PetscDataType type;
310: PetscScalar *array;
311: PetscInt bs, n, fid;
312: char name[PETSC_MAX_PATH_LEN];
313: PetscMPIInt size;
314: PetscBool iscuda, iskokkos, iship;
316: PetscFunctionBegin;
317: if (!swarm->issetup) PetscCall(DMSetUp(dm));
318: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
319: PetscCall(DMSwarmGetFieldInfo(dm, fieldname, &bs, &type));
320: PetscCheck(type == PETSC_SCALAR, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Field %s must have type PETSC_SCALAR, not %s", fieldname, PetscDataTypes[type]);
321: PetscCall(DMSwarmGetField(dm, fieldname, NULL, NULL, (void **)&array));
323: PetscCallMPI(MPI_Comm_size(comm, &size));
324: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
325: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
326: PetscCall(PetscStrcmp(dm->vectype, VECHIP, &iship));
327: PetscCall(VecCreate(comm, vec));
328: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
329: PetscCall(VecSetBlockSize(*vec, bs));
330: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
331: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
332: else if (iship) PetscCall(VecSetType(*vec, VECHIP));
333: else PetscCall(VecSetType(*vec, VECSTANDARD));
334: PetscCall(VecPlaceArray(*vec, array));
336: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
337: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
339: /* Set guard */
340: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
341: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
343: PetscCall(VecSetDM(*vec, dm));
344: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: static PetscErrorCode DMSwarmDestroyVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], Vec *vec)
349: {
350: DM_Swarm *swarm = (DM_Swarm *)sw->data;
351: const PetscScalar *array;
352: PetscInt bs, n, id = 0, cid = -2;
353: PetscBool flg;
355: PetscFunctionBegin;
356: for (PetscInt f = 0; f < Nf; ++f) {
357: PetscInt fid;
359: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
360: id += fid;
361: }
362: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cid, flg));
363: (void)flg; /* avoid compiler warning */
364: PetscCheck(cid == id, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)! %" PetscInt_FMT " != %" PetscInt_FMT, fieldnames[0], cid, id);
365: PetscCall(VecGetLocalSize(*vec, &n));
366: PetscCall(VecGetBlockSize(*vec, &bs));
367: n /= bs;
368: PetscCheck(n == swarm->db->L, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
369: PetscCall(VecGetArrayRead(*vec, &array));
370: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
371: void *farray;
372: PetscDataType ftype;
373: PetscInt fbs;
375: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
376: 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]);
377: PetscCall(DMSwarmGetField(sw, fieldnames[f], NULL, NULL, &farray));
378: PetscCheck(off + fbs <= bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid blocksize %" PetscInt_FMT " < %" PetscInt_FMT, bs, off + fbs);
379: if (ftype == PETSC_REAL) {
380: PetscReal *rarray = (PetscReal *)farray;
382: for (PetscInt i = 0; i < n; ++i) {
383: for (PetscInt b = 0; b < fbs; ++b) rarray[i * fbs + b] = PetscRealPart(array[i * bs + off + b]);
384: }
385: } else {
386: PetscScalar *sarray = (PetscScalar *)farray;
388: for (PetscInt i = 0; i < n; ++i) {
389: for (PetscInt b = 0; b < fbs; ++b) sarray[i * fbs + b] = array[i * bs + off + b];
390: }
391: }
392: off += fbs;
393: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], NULL, NULL, &farray));
394: }
395: PetscCall(VecRestoreArrayRead(*vec, &array));
396: PetscCall(VecDestroy(vec));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: static PetscErrorCode DMSwarmCreateVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], MPI_Comm comm, Vec *vec)
401: {
402: DM_Swarm *swarm = (DM_Swarm *)sw->data;
403: PetscScalar *array;
404: PetscInt n, bs = 0, id = 0;
405: char name[PETSC_MAX_PATH_LEN];
407: PetscFunctionBegin;
408: if (!swarm->issetup) PetscCall(DMSetUp(sw));
409: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
410: for (PetscInt f = 0; f < Nf; ++f) {
411: PetscDataType ftype;
412: PetscInt fbs;
414: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
415: 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]);
416: bs += fbs;
417: }
419: PetscCall(VecCreate(comm, vec));
420: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
421: PetscCall(VecSetBlockSize(*vec, bs));
422: PetscCall(VecSetType(*vec, sw->vectype));
424: PetscCall(VecGetArrayWrite(*vec, &array));
425: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
426: void *farray;
427: PetscDataType ftype;
428: PetscInt fbs;
430: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, &farray));
431: if (ftype == PETSC_REAL) {
432: const PetscReal *rarray = (const PetscReal *)farray;
434: for (PetscInt i = 0; i < n; ++i) {
435: for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = rarray[i * fbs + b];
436: }
437: } else {
438: const PetscScalar *sarray = (const PetscScalar *)farray;
440: for (PetscInt i = 0; i < n; ++i) {
441: for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = sarray[i * fbs + b];
442: }
443: }
444: off += fbs;
445: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, &farray));
446: }
447: PetscCall(VecRestoreArrayWrite(*vec, &array));
449: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
450: for (PetscInt f = 0; f < Nf; ++f) {
451: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
452: PetscCall(PetscStrlcat(name, fieldnames[f], PETSC_MAX_PATH_LEN));
453: }
454: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
456: for (PetscInt f = 0; f < Nf; ++f) {
457: PetscInt fid;
459: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
460: id += fid;
461: }
462: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, id));
464: PetscCall(VecSetDM(*vec, sw));
465: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@
470: DMSwarmPreallocateMassMatrix - Preallocate particle mass matrix between a `DMSWARM` and a `DMPLEX`
472: Collective
474: Input Parameters:
475: + dmc - The `DMSWARM` object
476: . dmf - The `DMPLEX` object
477: . mass - The mass matrix to preallocate
478: . rStart - The starting row index for this process
479: . maxC - The maximum number of columns per row
480: - ctx - The user context
482: Level: developer
484: .seealso: [](ch_unstructured), `DM`, `DMSWARM`, `DMPLEX`, `DMCreateMassMatrix()`, `DMSwarmFillMassMatrix()`
485: @*/
486: PetscErrorCode DMSwarmPreallocateMassMatrix(DM dmc, DM dmf, Mat mass, PetscInt *rStart, PetscInt *maxC, PetscCtx ctx)
487: {
488: MPI_Comm comm;
489: PetscHSetIJ ht;
490: PetscDS ds;
491: PetscSection fsection, globalFSection;
492: PetscLayout rLayout, colLayout;
493: PetscInt *dnz, *onz;
494: PetscInt cStart, cEnd, locRows, locCols, colStart, colEnd, Nf, totNc = 0;
496: PetscFunctionBegin;
500: PetscAssertPointer(rStart, 4);
501: PetscAssertPointer(maxC, 5);
502: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
503: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
504: PetscCall(DMGetLocalSection(dmf, &fsection));
505: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
506: PetscCall(DMGetDS(dmf, &ds));
507: PetscCall(PetscDSGetNumFields(ds, &Nf));
508: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
509: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
510: PetscCall(PetscHSetIJCreate(&ht));
512: PetscCall(PetscLayoutCreate(comm, &colLayout));
513: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
514: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
515: PetscCall(PetscLayoutSetUp(colLayout));
516: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
517: PetscCall(PetscLayoutDestroy(&colLayout));
519: PetscCall(PetscLayoutCreate(comm, &rLayout));
520: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
521: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
522: PetscCall(PetscLayoutSetUp(rLayout));
523: PetscCall(PetscLayoutGetRange(rLayout, rStart, NULL));
524: PetscCall(PetscLayoutDestroy(&rLayout));
526: for (PetscInt field = 0; field < Nf; ++field) {
527: PetscObject obj;
528: PetscClassId id;
529: PetscInt Nc;
531: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
532: PetscCall(PetscObjectGetClassId(obj, &id));
533: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
534: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
535: totNc += Nc;
536: }
537: PetscCall(DMSwarmSortGetAccess(dmc));
538: for (PetscInt field = 0; field < Nf; ++field) {
539: PetscObject obj;
540: PetscClassId id;
541: PetscInt Nc;
543: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
544: PetscCall(PetscObjectGetClassId(obj, &id));
545: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
546: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
548: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
549: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
550: PetscInt numFIndices, numCIndices;
552: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
553: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
554: *maxC = PetscMax(*maxC, numCIndices);
555: {
556: PetscHashIJKey key;
557: PetscBool missing;
558: for (PetscInt i = 0; i < numFIndices; ++i) {
559: key.j = findices[i]; /* global column (from Plex) */
560: if (key.j >= 0) {
561: /* Get indices for coarse elements */
562: for (PetscInt j = 0; j < numCIndices; ++j) {
563: for (PetscInt c = 0; c < Nc; ++c) {
564: // TODO Need field offset on particle here
565: key.i = cindices[j] * totNc + c + *rStart; /* global cols (from Swarm) */
566: if (key.i < 0) continue;
567: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
568: PetscCheck(missing, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
569: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - *rStart];
570: else ++onz[key.i - *rStart];
571: }
572: }
573: }
574: }
575: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
576: }
577: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
578: }
579: }
580: PetscCall(DMSwarmSortRestoreAccess(dmc));
581: PetscCall(PetscHSetIJDestroy(&ht));
582: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
583: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
584: PetscCall(PetscFree2(dnz, onz));
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
588: /*@
589: DMSwarmFillMassMatrix - Assemble the particle mass matrix between a `DMSWARM` and a `DMPLEX`
591: Collective
593: Input Parameters:
594: + dmc - The `DMSWARM` object
595: . dmf - The `DMPLEX` object
596: . mass - The mass matrix to fill up
597: . rStart - The starting row index for this process
598: . maxC - The maximum number of columns per row
599: . useDeltaFunction - Flag to use a delta function for the particle shape function
600: . Nfc - The number of swarm coordinate fields
601: . bs - The block size for each swarm coordinate field
602: . coordVals - The values for each particle for each swarm coordinate field
603: - ctx - The user context
605: Level: developer
607: .seealso: [](ch_unstructured), `DM`, `DMSWARM`, `DMPLEX`, `DMCreateMassMatrix()`, `DMSwarmPreallocateMassMatrix()`
608: @*/
609: PetscErrorCode DMSwarmFillMassMatrix(DM dmc, DM dmf, Mat mass, PetscInt rStart, PetscInt maxC, PetscBool useDeltaFunction, PetscInt Nfc, const PetscInt bs[], PetscReal *coordVals[], PetscCtx ctx)
610: {
611: const char *name = "Mass Matrix";
612: PetscDS ds;
613: PetscSection fsection, globalFSection;
614: PetscInt dim, cStart, cEnd, Nf, totDim, totNc = 0, *rowIDXs;
615: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
616: PetscScalar *elemMat;
618: PetscFunctionBegin;
622: PetscCall(DMGetCoordinateDim(dmf, &dim));
623: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
624: PetscCall(DMGetLocalSection(dmf, &fsection));
625: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
626: PetscCall(DMGetDS(dmf, &ds));
627: PetscCall(PetscDSGetNumFields(ds, &Nf));
628: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
629: for (PetscInt field = 0; field < Nf; ++field) {
630: PetscObject obj;
631: PetscClassId id;
632: PetscInt Nc;
634: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
635: PetscCall(PetscObjectGetClassId(obj, &id));
636: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
637: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
638: totNc += Nc;
639: }
641: PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
642: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
643: PetscCall(DMSwarmSortGetAccess(dmc));
644: for (PetscInt field = 0; field < Nf; ++field) {
645: PetscTabulation Tcoarse;
646: PetscObject obj;
647: PetscClassId id;
648: PetscInt Nc;
650: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
651: PetscCall(PetscObjectGetClassId(obj, &id));
652: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
653: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
655: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
656: PetscInt *findices, *cindices;
657: PetscInt numFIndices, numCIndices;
659: /* TODO: Use DMField instead of assuming affine */
660: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
661: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
662: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
663: for (PetscInt j = 0; j < numCIndices; ++j) {
664: PetscReal xr[8];
665: PetscInt off = 0;
667: for (PetscInt i = 0; i < Nfc; ++i) {
668: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[j] * bs[i] + b];
669: }
670: 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);
671: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, xr, &xi[j * dim]);
672: }
673: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
674: else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
675: /* Get elemMat entries by multiplying by weight */
676: PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
677: for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
678: for (PetscInt j = 0; j < numCIndices; ++j) {
679: for (PetscInt c = 0; c < Nc; ++c) {
680: // TODO Need field offset on particle and field here
681: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
682: elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
683: }
684: }
685: }
686: for (PetscInt j = 0; j < numCIndices; ++j)
687: // TODO Need field offset on particle here
688: for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
689: if (0) PetscCall(DMPrintCellMatrix(cell, name, numCIndices * Nc, numFIndices, elemMat));
690: PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
691: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
692: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
693: PetscCall(PetscTabulationDestroy(&Tcoarse));
694: }
695: }
696: PetscCall(DMSwarmSortRestoreAccess(dmc));
697: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
698: PetscCall(PetscFree3(v0, J, invJ));
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
704: \hat f = \sum_i f_i \phi_i
706: and a particle function is given by
708: f = \sum_i w_i \delta(x - x_i)
710: then we want to require that
712: M \hat f = M_p f
714: where the particle mass matrix is given by
716: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
718: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
719: his integral. We allow this with the boolean flag.
720: */
721: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, PetscCtx ctx)
722: {
723: DMSwarmCellDM celldm;
724: PetscInt rStart, maxC = 0;
725: PetscInt Nfc;
726: const char **coordFields;
727: PetscReal **coordVals;
728: PetscInt *bs;
730: PetscFunctionBegin;
731: PetscCall(DMSwarmPreallocateMassMatrix(dmc, dmf, mass, &rStart, &maxC, ctx));
733: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
734: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
735: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
736: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
737: PetscCall(DMSwarmFillMassMatrix(dmc, dmf, mass, rStart, maxC, useDeltaFunction, Nfc, bs, coordVals, ctx));
738: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
739: PetscCall(PetscFree2(coordVals, bs));
741: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
742: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: /* Returns empty matrix for use with SNES FD */
747: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
748: {
749: Vec field;
750: PetscInt size;
752: PetscFunctionBegin;
753: PetscCall(DMGetGlobalVector(sw, &field));
754: PetscCall(VecGetLocalSize(field, &size));
755: PetscCall(DMRestoreGlobalVector(sw, &field));
756: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
757: PetscCall(MatSetFromOptions(*m));
758: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
759: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
760: PetscCall(MatZeroEntries(*m));
761: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
762: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
763: PetscCall(MatShift(*m, 1.0));
764: PetscCall(MatSetDM(*m, sw));
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: /* FEM cols, Particle rows */
769: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
770: {
771: DMSwarmCellDM celldm;
772: PetscSection gsf;
773: PetscInt m, n, Np, bs;
774: void *ctx;
776: PetscFunctionBegin;
777: PetscCall(DMSwarmGetCellDMActive(dmCoarse, &celldm));
778: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)dmCoarse), PETSC_ERR_USER, "Active cell DM does not define any fields");
779: PetscCall(DMGetGlobalSection(dmFine, &gsf));
780: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
781: PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
782: PetscCall(DMSwarmCellDMGetBlockSize(celldm, dmCoarse, &bs));
783: n = Np * bs;
784: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
785: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
786: PetscCall(MatSetType(*mass, dmCoarse->mattype));
787: PetscCall(DMGetApplicationContext(dmFine, &ctx));
789: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
790: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, PetscCtx ctx)
795: {
796: const char *name = "Mass Matrix Square";
797: MPI_Comm comm;
798: DMSwarmCellDM celldm;
799: PetscDS prob;
800: PetscSection fsection, globalFSection;
801: PetscHSetIJ ht;
802: PetscLayout rLayout, colLayout;
803: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
804: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
805: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
806: PetscScalar *elemMat, *elemMatSq;
807: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0;
808: const char **coordFields;
809: PetscReal **coordVals;
810: PetscInt *bs;
812: PetscFunctionBegin;
813: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
814: PetscCall(DMGetCoordinateDim(dmf, &cdim));
815: PetscCall(DMGetDS(dmf, &prob));
816: PetscCall(PetscDSGetNumFields(prob, &Nf));
817: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
818: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
819: PetscCall(DMGetLocalSection(dmf, &fsection));
820: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
821: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
822: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
824: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
825: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
826: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
828: PetscCall(PetscLayoutCreate(comm, &colLayout));
829: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
830: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
831: PetscCall(PetscLayoutSetUp(colLayout));
832: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
833: PetscCall(PetscLayoutDestroy(&colLayout));
835: PetscCall(PetscLayoutCreate(comm, &rLayout));
836: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
837: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
838: PetscCall(PetscLayoutSetUp(rLayout));
839: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
840: PetscCall(PetscLayoutDestroy(&rLayout));
842: PetscCall(DMPlexGetDepth(dmf, &depth));
843: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
844: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
845: PetscCall(PetscMalloc1(maxAdjSize, &adj));
847: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
848: PetscCall(PetscHSetIJCreate(&ht));
849: /* Count nonzeros
850: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
851: */
852: PetscCall(DMSwarmSortGetAccess(dmc));
853: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
854: PetscInt *cindices;
855: PetscInt numCIndices;
856: #if 0
857: PetscInt adjSize = maxAdjSize, a, j;
858: #endif
860: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
861: maxC = PetscMax(maxC, numCIndices);
862: /* Diagonal block */
863: for (PetscInt i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
864: #if 0
865: /* Off-diagonal blocks */
866: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
867: for (a = 0; a < adjSize; ++a) {
868: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
869: const PetscInt ncell = adj[a];
870: PetscInt *ncindices;
871: PetscInt numNCIndices;
873: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
874: {
875: PetscHashIJKey key;
876: PetscBool missing;
878: for (i = 0; i < numCIndices; ++i) {
879: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
880: if (key.i < 0) continue;
881: for (j = 0; j < numNCIndices; ++j) {
882: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
883: if (key.j < 0) continue;
884: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
885: if (missing) {
886: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
887: else ++onz[key.i - rStart];
888: }
889: }
890: }
891: }
892: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
893: }
894: }
895: #endif
896: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
897: }
898: PetscCall(PetscHSetIJDestroy(&ht));
899: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
900: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
901: PetscCall(PetscFree2(dnz, onz));
902: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
903: /* Fill in values
904: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
905: Start just by producing block diagonal
906: Could loop over adjacent cells
907: Produce neighboring element matrix
908: TODO Determine which columns and rows correspond to shared dual vector
909: Do MatMatMult with rectangular matrices
910: Insert block
911: */
912: for (PetscInt field = 0; field < Nf; ++field) {
913: PetscTabulation Tcoarse;
914: PetscObject obj;
915: PetscInt Nc;
917: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
918: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
919: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
920: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
921: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
922: PetscInt *findices, *cindices;
923: PetscInt numFIndices, numCIndices;
925: /* TODO: Use DMField instead of assuming affine */
926: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
927: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
928: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
929: for (PetscInt p = 0; p < numCIndices; ++p) {
930: PetscReal xr[8];
931: PetscInt off = 0;
933: for (PetscInt i = 0; i < Nfc; ++i) {
934: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[p] * bs[i] + b];
935: }
936: 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);
937: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[p * cdim]);
938: }
939: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
940: /* Get elemMat entries by multiplying by weight */
941: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
942: for (PetscInt i = 0; i < numFIndices; ++i) {
943: for (PetscInt p = 0; p < numCIndices; ++p) {
944: for (PetscInt c = 0; c < Nc; ++c) {
945: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
946: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
947: }
948: }
949: }
950: PetscCall(PetscTabulationDestroy(&Tcoarse));
951: for (PetscInt p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
952: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
953: /* Block diagonal */
954: if (numCIndices) {
955: PetscBLASInt blasn, blask;
956: PetscScalar one = 1.0, zero = 0.0;
958: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
959: PetscCall(PetscBLASIntCast(numFIndices, &blask));
960: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
961: }
962: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
963: /* TODO off-diagonal */
964: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
965: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
966: }
967: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
968: }
969: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
970: PetscCall(PetscFree(adj));
971: PetscCall(DMSwarmSortRestoreAccess(dmc));
972: PetscCall(PetscFree3(v0, J, invJ));
973: PetscCall(PetscFree2(coordVals, bs));
974: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
975: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: /*@
980: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
982: Collective
984: Input Parameters:
985: + dmCoarse - a `DMSWARM`
986: - dmFine - a `DMPLEX`
988: Output Parameter:
989: . mass - the square of the particle mass matrix
991: Level: advanced
993: Note:
994: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
995: future to compute the full normal equations.
997: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
998: @*/
999: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
1000: {
1001: PetscInt n;
1002: void *ctx;
1004: PetscFunctionBegin;
1005: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
1006: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
1007: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
1008: PetscCall(MatSetType(*mass, dmCoarse->mattype));
1009: PetscCall(DMGetApplicationContext(dmFine, &ctx));
1011: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
1012: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: /* 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
1018: \int_X \psi_i \nabla \cdot \hat f = \int_X \psi_i \nabla \cdot f
1020: and then integrate by parts
1022: \int_X \nabla \psi_i \cdot \hat f = \int_X \nabla \psi_i \cdot f
1024: where \psi is from a scalar FE space. If a finite element interpolant is given by
1026: \hat f^c = \sum_i f_i \phi^c_i
1028: and a particle function is given by
1030: f^c = \sum_p f^c_p \delta(x - x_p)
1032: then we want to require that
1034: D_f \hat f = D_p f
1036: where the gradient matrices are given by
1038: (D_f)_{i(jc)} = \int \partial_c \psi_i \phi_j
1039: (D_p)_{i(jc)} = \int \partial_c \psi_i \delta(x - x_j)
1041: Thus we need two finite element spaces, a scalar and a vector. The vector space holds the representer for the
1042: vector particle field. The scalar space holds the output of D_p or D_f, which is the weak divergence of the field.
1044: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
1045: his integral. We allow this with the boolean flag.
1046: */
1047: static PetscErrorCode DMSwarmComputeGradientMatrix_Private(DM sw, DM dm, Mat derv, PetscBool useDeltaFunction, PetscCtx ctx)
1048: {
1049: const char *name = "Derivative Matrix";
1050: MPI_Comm comm;
1051: DMSwarmCellDM celldm;
1052: PetscDS ds;
1053: PetscSection fsection, globalFSection;
1054: PetscLayout rLayout;
1055: PetscInt locRows, rStart, *rowIDXs;
1056: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
1057: PetscScalar *elemMat;
1058: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxNpc = 0, totNc = 0;
1059: const char **coordFields;
1060: PetscReal **coordVals;
1061: PetscInt *bs;
1063: PetscFunctionBegin;
1064: PetscCall(PetscObjectGetComm((PetscObject)derv, &comm));
1065: PetscCall(DMGetCoordinateDim(dm, &cdim));
1066: PetscCall(DMGetDS(dm, &ds));
1067: PetscCall(PetscDSGetNumFields(ds, &Nf));
1068: PetscCheck(Nf == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
1069: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
1070: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
1071: PetscCall(DMGetLocalSection(dm, &fsection));
1072: PetscCall(DMGetGlobalSection(dm, &globalFSection));
1073: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1074: PetscCall(MatGetLocalSize(derv, &locRows, NULL));
1076: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1077: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1078: PetscCheck(Nfc == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
1079: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
1081: PetscCall(PetscLayoutCreate(comm, &rLayout));
1082: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
1083: PetscCall(PetscLayoutSetBlockSize(rLayout, cdim));
1084: PetscCall(PetscLayoutSetUp(rLayout));
1085: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
1086: PetscCall(PetscLayoutDestroy(&rLayout));
1088: for (PetscInt field = 0; field < Nf; ++field) {
1089: PetscObject obj;
1090: PetscInt Nc;
1092: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
1093: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
1094: totNc += Nc;
1095: }
1096: PetscCheck(totNc == 1, comm, PETSC_ERR_ARG_WRONG, "The number of field components %" PetscInt_FMT " != 1", totNc);
1097: /* count non-zeros */
1098: PetscCall(DMSwarmSortGetAccess(sw));
1099: for (PetscInt field = 0; field < Nf; ++field) {
1100: PetscObject obj;
1101: PetscInt Nc;
1103: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
1104: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
1105: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1106: PetscInt *pind;
1107: PetscInt Npc;
1109: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1110: maxNpc = PetscMax(maxNpc, Npc);
1111: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1112: }
1113: }
1114: PetscCall(PetscMalloc3(maxNpc * cdim * totDim, &elemMat, maxNpc * cdim, &rowIDXs, maxNpc * cdim, &xi));
1115: for (PetscInt field = 0; field < Nf; ++field) {
1116: PetscTabulation Tcoarse;
1117: PetscFE fe;
1119: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1120: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1121: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1122: PetscInt *findices, *pind;
1123: PetscInt numFIndices, Npc;
1125: /* TODO: Use DMField instead of assuming affine */
1126: PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
1127: PetscCall(DMPlexGetClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1128: PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1129: for (PetscInt j = 0; j < Npc; ++j) {
1130: PetscReal xr[8];
1131: PetscInt off = 0;
1133: for (PetscInt i = 0; i < Nfc; ++i) {
1134: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][pind[j] * bs[i] + b];
1135: }
1136: 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);
1137: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[j * cdim]);
1138: }
1139: PetscCall(PetscFECreateTabulation(fe, 1, Npc, xi, 1, &Tcoarse));
1140: /* Get elemMat entries by multiplying by weight */
1141: PetscCall(PetscArrayzero(elemMat, Npc * cdim * totDim));
1142: for (PetscInt i = 0; i < numFIndices; ++i) {
1143: for (PetscInt j = 0; j < Npc; ++j) {
1144: /* D[((p*pdim + i)*Nc + c)*cdim + d] is the value at point p for basis function i, component c, derivative d */
1145: for (PetscInt d = 0; d < cdim; ++d) {
1146: xi[d] = 0.;
1147: for (PetscInt e = 0; e < cdim; ++e) xi[d] += invJ[e * cdim + d] * Tcoarse->T[1][(j * numFIndices + i) * cdim + e];
1148: elemMat[(j * cdim + d) * numFIndices + i] += xi[d] * (useDeltaFunction ? 1.0 : detJ);
1149: }
1150: }
1151: }
1152: for (PetscInt j = 0; j < Npc; ++j)
1153: for (PetscInt d = 0; d < cdim; ++d) rowIDXs[j * cdim + d] = pind[j] * cdim + d + rStart;
1154: if (0) PetscCall(DMPrintCellMatrix(cell, name, Npc * cdim, numFIndices, elemMat));
1155: PetscCall(MatSetValues(derv, Npc * cdim, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
1156: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1157: PetscCall(DMPlexRestoreClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1158: PetscCall(PetscTabulationDestroy(&Tcoarse));
1159: }
1160: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1161: }
1162: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
1163: PetscCall(DMSwarmSortRestoreAccess(sw));
1164: PetscCall(PetscFree3(v0, J, invJ));
1165: PetscCall(PetscFree2(coordVals, bs));
1166: PetscCall(MatAssemblyBegin(derv, MAT_FINAL_ASSEMBLY));
1167: PetscCall(MatAssemblyEnd(derv, MAT_FINAL_ASSEMBLY));
1168: PetscFunctionReturn(PETSC_SUCCESS);
1169: }
1171: /* FEM cols: this is a scalar space
1172: Particle rows: this is a vector space that contracts with the derivative
1173: */
1174: static PetscErrorCode DMCreateGradientMatrix_Swarm(DM sw, DM dm, Mat *derv)
1175: {
1176: DMSwarmCellDM celldm;
1177: PetscSection gs;
1178: PetscInt cdim, m, n, Np, bs;
1179: void *ctx;
1180: MPI_Comm comm;
1182: PetscFunctionBegin;
1183: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1184: PetscCall(DMGetCoordinateDim(dm, &cdim));
1185: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1186: PetscCheck(celldm->Nf, comm, PETSC_ERR_USER, "Active cell DM does not define any fields");
1187: PetscCall(DMGetGlobalSection(dm, &gs));
1188: PetscCall(PetscSectionGetConstrainedStorageSize(gs, &n));
1189: PetscCall(DMSwarmGetLocalSize(sw, &Np));
1190: PetscCall(DMSwarmCellDMGetBlockSize(celldm, sw, &bs));
1191: PetscCheck(cdim == bs, comm, PETSC_ERR_ARG_WRONG, "Coordinate dimension %" PetscInt_FMT " != %" PetscInt_FMT " swarm field block size", cdim, bs);
1192: m = Np * bs;
1193: PetscCall(MatCreate(PetscObjectComm((PetscObject)sw), derv));
1194: PetscCall(PetscObjectSetName((PetscObject)*derv, "Swarm Derivative Matrix"));
1195: PetscCall(MatSetSizes(*derv, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
1196: PetscCall(MatSetType(*derv, sw->mattype));
1197: PetscCall(DMGetApplicationContext(dm, &ctx));
1199: PetscCall(DMSwarmComputeGradientMatrix_Private(sw, dm, *derv, PETSC_TRUE, ctx));
1200: PetscCall(MatViewFromOptions(*derv, NULL, "-gradient_mat_view"));
1201: PetscFunctionReturn(PETSC_SUCCESS);
1202: }
1204: /*@
1205: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
1207: Collective
1209: Input Parameters:
1210: + dm - a `DMSWARM`
1211: - fieldname - the textual name given to a registered field
1213: Output Parameter:
1214: . vec - the vector
1216: Level: beginner
1218: Note:
1219: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
1221: The field must have data type `PETSC_SCALAR`.
1223: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
1224: @*/
1225: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1226: {
1227: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1229: PetscFunctionBegin;
1231: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
1232: PetscFunctionReturn(PETSC_SUCCESS);
1233: }
1235: /*@
1236: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
1238: Collective
1240: Input Parameters:
1241: + dm - a `DMSWARM`
1242: - fieldname - the textual name given to a registered field
1244: Output Parameter:
1245: . vec - the vector
1247: Level: beginner
1249: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1250: @*/
1251: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1252: {
1253: PetscFunctionBegin;
1255: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1256: PetscFunctionReturn(PETSC_SUCCESS);
1257: }
1259: /*@
1260: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
1262: Collective
1264: Input Parameters:
1265: + dm - a `DMSWARM`
1266: - fieldname - the textual name given to a registered field
1268: Output Parameter:
1269: . vec - the vector
1271: Level: beginner
1273: Note:
1274: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1276: The field must have data type `PETSC_SCALAR`.
1278: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1279: @*/
1280: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1281: {
1282: PetscFunctionBegin;
1283: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, PETSC_COMM_SELF, vec));
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: /*@
1288: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
1290: Collective
1292: Input Parameters:
1293: + dm - a `DMSWARM`
1294: - fieldname - the textual name given to a registered field
1296: Output Parameter:
1297: . vec - the vector
1299: Level: beginner
1301: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
1302: @*/
1303: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1304: {
1305: PetscFunctionBegin;
1307: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1308: PetscFunctionReturn(PETSC_SUCCESS);
1309: }
1311: /*@
1312: DMSwarmCreateGlobalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1314: Collective
1316: Input Parameters:
1317: + dm - a `DMSWARM`
1318: . Nf - the number of fields
1319: - fieldnames - the textual names given to the registered fields
1321: Output Parameter:
1322: . vec - the vector
1324: Level: beginner
1326: Notes:
1327: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromFields()`.
1329: All fields must have data type `PETSC_REAL` or `PETSC_SCALAR`. When PETSc is configured with complex scalars,
1330: only the real part of vector entries is copied back to `PETSC_REAL` fields.
1332: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateGlobalVectorFromField()`
1334: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromFields()`
1335: @*/
1336: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1337: {
1338: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1340: PetscFunctionBegin;
1342: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1343: PetscFunctionReturn(PETSC_SUCCESS);
1344: }
1346: /*@
1347: DMSwarmDestroyGlobalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1349: Collective
1351: Input Parameters:
1352: + dm - a `DMSWARM`
1353: . Nf - the number of fields
1354: - fieldnames - the textual names given to the registered fields
1356: Output Parameter:
1357: . vec - the vector
1359: Level: beginner
1361: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1362: @*/
1363: PetscErrorCode DMSwarmDestroyGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1364: {
1365: PetscFunctionBegin;
1367: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1368: PetscFunctionReturn(PETSC_SUCCESS);
1369: }
1371: /*@
1372: DMSwarmCreateLocalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1374: Collective
1376: Input Parameters:
1377: + dm - a `DMSWARM`
1378: . Nf - the number of fields
1379: - fieldnames - the textual names given to the registered fields
1381: Output Parameter:
1382: . vec - the vector
1384: Level: beginner
1386: Notes:
1387: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1389: All fields must have data type `PETSC_REAL` or `PETSC_SCALAR`. When PETSc is configured with complex scalars,
1390: only the real part of vector entries is copied back to `PETSC_REAL` fields.
1392: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateLocalVectorFromField()`
1394: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1395: @*/
1396: PetscErrorCode DMSwarmCreateLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1397: {
1398: PetscFunctionBegin;
1399: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, PETSC_COMM_SELF, vec));
1400: PetscFunctionReturn(PETSC_SUCCESS);
1401: }
1403: /*@
1404: DMSwarmDestroyLocalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1406: Collective
1408: Input Parameters:
1409: + dm - a `DMSWARM`
1410: . Nf - the number of fields
1411: - fieldnames - the textual names given to the registered fields
1413: Output Parameter:
1414: . vec - the vector
1416: Level: beginner
1418: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromFields()`
1419: @*/
1420: PetscErrorCode DMSwarmDestroyLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1421: {
1422: PetscFunctionBegin;
1424: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1425: PetscFunctionReturn(PETSC_SUCCESS);
1426: }
1428: /*@
1429: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
1431: Collective
1433: Input Parameter:
1434: . dm - a `DMSWARM`
1436: Level: beginner
1438: Note:
1439: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
1441: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1442: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1443: @*/
1444: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
1445: {
1446: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1448: PetscFunctionBegin;
1449: if (!swarm->field_registration_initialized) {
1450: swarm->field_registration_initialized = PETSC_TRUE;
1451: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
1452: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
1453: }
1454: PetscFunctionReturn(PETSC_SUCCESS);
1455: }
1457: /*@
1458: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
1460: Collective
1462: Input Parameter:
1463: . dm - a `DMSWARM`
1465: Level: beginner
1467: Note:
1468: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
1470: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1471: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1472: @*/
1473: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
1474: {
1475: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1477: PetscFunctionBegin;
1478: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
1479: swarm->field_registration_finalized = PETSC_TRUE;
1480: PetscFunctionReturn(PETSC_SUCCESS);
1481: }
1483: /*@
1484: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
1486: Not Collective
1488: Input Parameters:
1489: + sw - a `DMSWARM`
1490: . nlocal - the length of each registered field
1491: - buffer - the length of the buffer used to efficient dynamic re-sizing
1493: Level: beginner
1495: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
1496: @*/
1497: PetscErrorCode DMSwarmSetLocalSizes(DM sw, PetscInt nlocal, PetscInt buffer)
1498: {
1499: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1500: PetscMPIInt rank;
1501: PetscInt *rankval;
1503: PetscFunctionBegin;
1504: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
1505: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
1506: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
1508: // Initialize values in pid and rank placeholders
1509: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1510: PetscCall(DMSwarmGetField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1511: for (PetscInt p = 0; p < nlocal; p++) rankval[p] = rank;
1512: PetscCall(DMSwarmRestoreField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1513: /* TODO: [pid - use MPI_Scan] */
1514: PetscFunctionReturn(PETSC_SUCCESS);
1515: }
1517: /*@
1518: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
1520: Collective
1522: Input Parameters:
1523: + sw - a `DMSWARM`
1524: - dm - the `DM` to attach to the `DMSWARM`
1526: Level: beginner
1528: Note:
1529: The attached `DM` (dm) will be queried for point location and
1530: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1532: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
1533: @*/
1534: PetscErrorCode DMSwarmSetCellDM(DM sw, DM dm)
1535: {
1536: DMSwarmCellDM celldm;
1537: const char *name;
1538: char *coordName;
1540: PetscFunctionBegin;
1543: PetscCall(PetscStrallocpy(DMSwarmPICField_coor, &coordName));
1544: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 1, (const char **)&coordName, &celldm));
1545: PetscCall(PetscFree(coordName));
1546: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1547: PetscCall(DMSwarmAddCellDM(sw, celldm));
1548: PetscCall(DMSwarmCellDMDestroy(&celldm));
1549: PetscCall(DMSwarmSetCellDMActive(sw, name));
1550: PetscFunctionReturn(PETSC_SUCCESS);
1551: }
1553: /*@
1554: DMSwarmGetCellDM - Fetches the active cell `DM`
1556: Collective
1558: Input Parameter:
1559: . sw - a `DMSWARM`
1561: Output Parameter:
1562: . dm - the active `DM` for the `DMSWARM`
1564: Level: beginner
1566: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1567: @*/
1568: PetscErrorCode DMSwarmGetCellDM(DM sw, DM *dm)
1569: {
1570: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1571: DMSwarmCellDM celldm;
1573: PetscFunctionBegin;
1575: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)&celldm));
1576: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "There is no cell DM named %s in this Swarm", swarm->activeCellDM);
1577: PetscCall(DMSwarmCellDMGetDM(celldm, dm));
1578: PetscFunctionReturn(PETSC_SUCCESS);
1579: }
1581: /*@
1582: DMSwarmGetCellDMNames - Get the list of cell `DM` names
1584: Not collective
1586: Input Parameter:
1587: . sw - a `DMSWARM`
1589: Output Parameters:
1590: + Ndm - the number of `DMSwarmCellDM` in the `DMSWARM`
1591: - celldms - the name of each `DMSwarmCellDM`
1593: Level: beginner
1595: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMByName()`
1596: @*/
1597: PetscErrorCode DMSwarmGetCellDMNames(DM sw, PetscInt *Ndm, const char **celldms[])
1598: {
1599: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1600: PetscObjectList next = swarm->cellDMs;
1601: PetscInt n = 0;
1603: PetscFunctionBegin;
1605: PetscAssertPointer(Ndm, 2);
1606: PetscAssertPointer(celldms, 3);
1607: while (next) {
1608: next = next->next;
1609: ++n;
1610: }
1611: PetscCall(PetscMalloc1(n, celldms));
1612: next = swarm->cellDMs;
1613: n = 0;
1614: while (next) {
1615: (*celldms)[n] = (const char *)next->obj->name;
1616: next = next->next;
1617: ++n;
1618: }
1619: *Ndm = n;
1620: PetscFunctionReturn(PETSC_SUCCESS);
1621: }
1623: /*@
1624: DMSwarmSetCellDMActive - Activates a cell `DM` for a `DMSWARM`
1626: Collective
1628: Input Parameters:
1629: + sw - a `DMSWARM`
1630: - name - name of the cell `DM` to active for the `DMSWARM`
1632: Level: beginner
1634: Note:
1635: The attached `DM` (dmcell) will be queried for point location and
1636: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1638: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1639: @*/
1640: PetscErrorCode DMSwarmSetCellDMActive(DM sw, const char name[])
1641: {
1642: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1643: DMSwarmCellDM celldm;
1645: PetscFunctionBegin;
1647: PetscCall(PetscInfo(sw, "Setting cell DM to %s\n", name));
1648: PetscCall(PetscFree(swarm->activeCellDM));
1649: PetscCall(PetscStrallocpy(name, (char **)&swarm->activeCellDM));
1650: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1651: PetscFunctionReturn(PETSC_SUCCESS);
1652: }
1654: /*@
1655: DMSwarmGetCellDMActive - Returns the active cell `DM` for a `DMSWARM`
1657: Collective
1659: Input Parameter:
1660: . sw - a `DMSWARM`
1662: Output Parameter:
1663: . celldm - the active `DMSwarmCellDM`
1665: Level: beginner
1667: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1668: @*/
1669: PetscErrorCode DMSwarmGetCellDMActive(DM sw, DMSwarmCellDM *celldm)
1670: {
1671: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1673: PetscFunctionBegin;
1675: PetscAssertPointer(celldm, 2);
1676: PetscCheck(swarm->activeCellDM, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no active cell DM");
1677: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)celldm));
1678: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", swarm->activeCellDM);
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: /*@
1683: DMSwarmGetCellDMByName - Get a `DMSwarmCellDM` from its name
1685: Not collective
1687: Input Parameters:
1688: + sw - a `DMSWARM`
1689: - name - the name
1691: Output Parameter:
1692: . celldm - the `DMSwarmCellDM`, or `NULL` if the name is unknown
1694: Level: beginner
1696: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMNames()`
1697: @*/
1698: PetscErrorCode DMSwarmGetCellDMByName(DM sw, const char name[], DMSwarmCellDM *celldm)
1699: {
1700: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1702: PetscFunctionBegin;
1704: PetscAssertPointer(name, 2);
1705: PetscAssertPointer(celldm, 3);
1706: PetscCall(PetscObjectListFind(swarm->cellDMs, name, (PetscObject *)celldm));
1707: PetscFunctionReturn(PETSC_SUCCESS);
1708: }
1710: /*@
1711: DMSwarmAddCellDM - Adds a cell `DM` to the `DMSWARM`
1713: Collective
1715: Input Parameters:
1716: + sw - a `DMSWARM`
1717: - celldm - the `DMSwarmCellDM`
1719: Level: beginner
1721: Note:
1722: Cell DMs with the same name will share the cellid field
1724: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPushCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1725: @*/
1726: PetscErrorCode DMSwarmAddCellDM(DM sw, DMSwarmCellDM celldm)
1727: {
1728: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1729: const char *name;
1730: PetscInt dim;
1731: PetscBool flg;
1732: MPI_Comm comm;
1734: PetscFunctionBegin;
1736: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1738: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1739: PetscCall(PetscObjectListAdd(&swarm->cellDMs, name, (PetscObject)celldm));
1740: PetscCall(DMGetDimension(sw, &dim));
1741: for (PetscInt f = 0; f < celldm->Nfc; ++f) {
1742: PetscCall(DMSwarmDataFieldStringInList(celldm->coordFields[f], swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1743: if (!flg) {
1744: PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->coordFields[f], dim, PETSC_REAL));
1745: } else {
1746: PetscDataType dt;
1747: PetscInt bs;
1749: PetscCall(DMSwarmGetFieldInfo(sw, celldm->coordFields[f], &bs, &dt));
1750: PetscCheck(bs == dim, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has blocksize %" PetscInt_FMT " != %" PetscInt_FMT " spatial dimension", celldm->coordFields[f], bs, dim);
1751: PetscCheck(dt == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has datatype %s != PETSC_REAL", celldm->coordFields[f], PetscDataTypes[dt]);
1752: }
1753: }
1754: // Assume that DMs with the same name share the cellid field
1755: PetscCall(DMSwarmDataFieldStringInList(celldm->cellid, swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1756: if (!flg) {
1757: PetscBool isShell, isDummy;
1758: const char *name;
1760: // Allow dummy DMSHELL (I don't think we should support this mode)
1761: PetscCall(PetscObjectTypeCompare((PetscObject)celldm->dm, DMSHELL, &isShell));
1762: PetscCall(PetscObjectGetName((PetscObject)celldm->dm, &name));
1763: PetscCall(PetscStrcmp(name, "dummy", &isDummy));
1764: if (!isShell || !isDummy) PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->cellid, 1, PETSC_INT));
1765: }
1766: PetscCall(DMSwarmSetCellDMActive(sw, name));
1767: PetscFunctionReturn(PETSC_SUCCESS);
1768: }
1770: /*@
1771: DMSwarmGetLocalSize - Retrieves the local length of fields registered
1773: Not Collective
1775: Input Parameter:
1776: . dm - a `DMSWARM`
1778: Output Parameter:
1779: . nlocal - the length of each registered field
1781: Level: beginner
1783: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
1784: @*/
1785: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
1786: {
1787: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1789: PetscFunctionBegin;
1790: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
1791: PetscFunctionReturn(PETSC_SUCCESS);
1792: }
1794: /*@
1795: DMSwarmGetSize - Retrieves the total length of fields registered
1797: Collective
1799: Input Parameter:
1800: . dm - a `DMSWARM`
1802: Output Parameter:
1803: . n - the total length of each registered field
1805: Level: beginner
1807: Note:
1808: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
1810: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
1811: @*/
1812: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
1813: {
1814: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1816: PetscFunctionBegin;
1817: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, n, NULL, NULL));
1818: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1819: PetscFunctionReturn(PETSC_SUCCESS);
1820: }
1822: /*@
1823: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
1825: Collective
1827: Input Parameters:
1828: + dm - a `DMSWARM`
1829: . fieldname - the textual name to identify this field
1830: . blocksize - the number of each data type
1831: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
1833: Level: beginner
1835: Notes:
1836: The textual name for each registered field must be unique.
1838: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1839: @*/
1840: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1841: {
1842: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1843: size_t size;
1845: PetscFunctionBegin;
1846: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1847: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1849: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1850: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1851: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1852: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1853: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1855: PetscCall(PetscDataTypeGetSize(type, &size));
1856: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1857: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1858: {
1859: DMSwarmDataField gfield;
1861: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1862: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1863: }
1864: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1865: PetscFunctionReturn(PETSC_SUCCESS);
1866: }
1868: /*@
1869: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
1871: Collective
1873: Input Parameters:
1874: + dm - a `DMSWARM`
1875: . fieldname - the textual name to identify this field
1876: - size - the size in bytes of the user struct of each data type
1878: Level: beginner
1880: Note:
1881: The textual name for each registered field must be unique.
1883: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1884: @*/
1885: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1886: {
1887: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1889: PetscFunctionBegin;
1890: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1891: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1892: PetscFunctionReturn(PETSC_SUCCESS);
1893: }
1895: /*@
1896: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1898: Collective
1900: Input Parameters:
1901: + dm - a `DMSWARM`
1902: . fieldname - the textual name to identify this field
1903: . size - the size in bytes of the user data type
1904: - blocksize - the number of each data type
1906: Level: beginner
1908: Note:
1909: The textual name for each registered field must be unique.
1911: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1912: @*/
1913: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1914: {
1915: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1917: PetscFunctionBegin;
1918: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1919: {
1920: DMSwarmDataField gfield;
1922: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1923: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1924: }
1925: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1926: PetscFunctionReturn(PETSC_SUCCESS);
1927: }
1929: /*@
1930: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1932: Not Collective, No Fortran Support
1934: Input Parameters:
1935: + dm - a `DMSWARM`
1936: - fieldname - the textual name to identify this field
1938: Output Parameters:
1939: + blocksize - the number of each data type
1940: . type - the data type
1941: - data - pointer to raw array
1943: Level: beginner
1945: Notes:
1946: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1948: Fortran Note:
1949: Only works for `type` of `PETSC_SCALAR`
1951: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1952: @*/
1953: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1954: {
1955: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1956: DMSwarmDataField gfield;
1958: PetscFunctionBegin;
1960: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1961: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1962: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1963: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1964: if (blocksize) *blocksize = gfield->bs;
1965: if (type) *type = gfield->petsc_type;
1966: PetscFunctionReturn(PETSC_SUCCESS);
1967: }
1969: /*@
1970: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1972: Not Collective
1974: Input Parameters:
1975: + dm - a `DMSWARM`
1976: - fieldname - the textual name to identify this field
1978: Output Parameters:
1979: + blocksize - the number of each data type
1980: . type - the data type
1981: - data - pointer to raw array
1983: Level: beginner
1985: Notes:
1986: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1988: Fortran Note:
1989: Only works for `type` of `PETSC_SCALAR`
1991: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1992: @*/
1993: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1994: {
1995: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1996: DMSwarmDataField gfield;
1998: PetscFunctionBegin;
2000: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
2001: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
2002: if (data) *data = NULL;
2003: PetscFunctionReturn(PETSC_SUCCESS);
2004: }
2006: /*@
2007: DMSwarmGetFieldInfo - Return the block size and data type of a registered `DMSWARM` field without accessing its data.
2009: Not Collective
2011: Input Parameters:
2012: + dm - a `DMSWARM`
2013: - fieldname - the name of the registered field
2015: Output Parameters:
2016: + blocksize - the number of entries of `type` per particle, or `NULL`
2017: - type - the `PetscDataType` of a single entry, or `NULL`
2019: Level: intermediate
2021: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`, `DMSwarmRestoreField()`, `DMSwarmRegisterPetscDatatypeField()`
2022: @*/
2023: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
2024: {
2025: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2026: DMSwarmDataField gfield;
2028: PetscFunctionBegin;
2030: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
2031: if (blocksize) *blocksize = gfield->bs;
2032: if (type) *type = gfield->petsc_type;
2033: PetscFunctionReturn(PETSC_SUCCESS);
2034: }
2036: /*@
2037: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
2039: Not Collective
2041: Input Parameter:
2042: . dm - a `DMSWARM`
2044: Level: beginner
2046: Notes:
2047: The new point will have all fields initialized to zero.
2049: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
2050: @*/
2051: PetscErrorCode DMSwarmAddPoint(DM dm)
2052: {
2053: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2055: PetscFunctionBegin;
2056: if (!swarm->issetup) PetscCall(DMSetUp(dm));
2057: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
2058: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
2059: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
2060: PetscFunctionReturn(PETSC_SUCCESS);
2061: }
2063: /*@
2064: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
2066: Not Collective
2068: Input Parameters:
2069: + dm - a `DMSWARM`
2070: - npoints - the number of new points to add
2072: Level: beginner
2074: Notes:
2075: The new point will have all fields initialized to zero.
2077: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
2078: @*/
2079: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
2080: {
2081: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2082: PetscInt nlocal;
2084: PetscFunctionBegin;
2085: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
2086: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
2087: nlocal = PetscMax(nlocal, 0) + npoints;
2088: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
2089: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
2090: PetscFunctionReturn(PETSC_SUCCESS);
2091: }
2093: /*@
2094: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
2096: Not Collective
2098: Input Parameter:
2099: . dm - a `DMSWARM`
2101: Level: beginner
2103: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
2104: @*/
2105: PetscErrorCode DMSwarmRemovePoint(DM dm)
2106: {
2107: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2109: PetscFunctionBegin;
2110: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
2111: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
2112: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
2113: PetscFunctionReturn(PETSC_SUCCESS);
2114: }
2116: /*@
2117: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
2119: Not Collective
2121: Input Parameters:
2122: + dm - a `DMSWARM`
2123: - idx - index of point to remove
2125: Level: beginner
2127: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2128: @*/
2129: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
2130: {
2131: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2133: PetscFunctionBegin;
2134: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
2135: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
2136: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
2137: PetscFunctionReturn(PETSC_SUCCESS);
2138: }
2140: /*@
2141: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
2143: Not Collective
2145: Input Parameters:
2146: + dm - a `DMSWARM`
2147: . pi - the index of the point to copy
2148: - pj - the point index where the copy should be located
2150: Level: beginner
2152: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2153: @*/
2154: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
2155: {
2156: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2158: PetscFunctionBegin;
2159: if (!swarm->issetup) PetscCall(DMSetUp(dm));
2160: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
2161: PetscFunctionReturn(PETSC_SUCCESS);
2162: }
2164: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
2165: {
2166: PetscFunctionBegin;
2167: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
2168: PetscFunctionReturn(PETSC_SUCCESS);
2169: }
2171: /*@
2172: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
2174: Collective
2176: Input Parameters:
2177: + dm - the `DMSWARM`
2178: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
2180: Level: advanced
2182: Notes:
2183: The `DM` will be modified to accommodate received points.
2184: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
2185: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
2187: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
2188: @*/
2189: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
2190: {
2191: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2193: PetscFunctionBegin;
2194: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
2195: switch (swarm->migrate_type) {
2196: case DMSWARM_MIGRATE_BASIC:
2197: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
2198: break;
2199: case DMSWARM_MIGRATE_DMCELLNSCATTER:
2200: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
2201: break;
2202: case DMSWARM_MIGRATE_DMCELLEXACT:
2203: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
2204: case DMSWARM_MIGRATE_USER:
2205: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
2206: default:
2207: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
2208: }
2209: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
2210: PetscCall(DMClearGlobalVectors(dm));
2211: PetscFunctionReturn(PETSC_SUCCESS);
2212: }
2214: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
2216: /*
2217: DMSwarmCollectViewCreate
2219: * Applies a collection method and gathers point neighbour points into dm
2221: Notes:
2222: Users should call DMSwarmCollectViewDestroy() after
2223: they have finished computations associated with the collected points
2224: */
2226: /*@
2227: DMSwarmCollectViewCreate - Applies a collection method and gathers points
2228: in neighbour ranks into the `DMSWARM`
2230: Collective
2232: Input Parameter:
2233: . dm - the `DMSWARM`
2235: Level: advanced
2237: Notes:
2238: Users should call `DMSwarmCollectViewDestroy()` after
2239: they have finished computations associated with the collected points
2241: Different collect methods are supported. See `DMSwarmSetCollectType()`.
2243: Developer Note:
2244: Create and Destroy routines create new objects that can get destroyed, they do not change the state
2245: of the current object.
2247: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
2248: @*/
2249: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
2250: {
2251: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2252: PetscInt ng;
2254: PetscFunctionBegin;
2255: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
2256: PetscCall(DMSwarmGetLocalSize(dm, &ng));
2257: switch (swarm->collect_type) {
2258: case DMSWARM_COLLECT_BASIC:
2259: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
2260: break;
2261: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
2262: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
2263: case DMSWARM_COLLECT_GENERAL:
2264: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
2265: default:
2266: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
2267: }
2268: swarm->collect_view_active = PETSC_TRUE;
2269: swarm->collect_view_reset_nlocal = ng;
2270: PetscFunctionReturn(PETSC_SUCCESS);
2271: }
2273: /*@
2274: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
2276: Collective
2278: Input Parameters:
2279: . dm - the `DMSWARM`
2281: Notes:
2282: Users should call `DMSwarmCollectViewCreate()` before this function is called.
2284: Level: advanced
2286: Developer Note:
2287: Create and Destroy routines create new objects that can get destroyed, they do not change the state
2288: of the current object.
2290: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
2291: @*/
2292: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
2293: {
2294: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2296: PetscFunctionBegin;
2297: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
2298: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
2299: swarm->collect_view_active = PETSC_FALSE;
2300: PetscFunctionReturn(PETSC_SUCCESS);
2301: }
2303: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
2304: {
2305: PetscInt dim;
2307: PetscFunctionBegin;
2308: PetscCall(DMSwarmSetNumSpecies(dm, 1));
2309: PetscCall(DMGetDimension(dm, &dim));
2310: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2311: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2312: PetscFunctionReturn(PETSC_SUCCESS);
2313: }
2315: /*@
2316: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
2318: Collective
2320: Input Parameters:
2321: + dm - the `DMSWARM`
2322: - Npc - The number of particles per cell in the cell `DM`
2324: Level: intermediate
2326: Notes:
2327: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
2328: one particle is in each cell, it is placed at the centroid.
2330: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
2331: @*/
2332: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
2333: {
2334: DM cdm;
2335: DMSwarmCellDM celldm;
2336: PetscRandom rnd;
2337: DMPolytopeType ct;
2338: PetscBool simplex;
2339: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
2340: PetscInt dim, d, cStart, cEnd, c, p, Nfc;
2341: const char **coordFields;
2343: PetscFunctionBeginUser;
2344: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
2345: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
2346: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
2348: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2349: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2350: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2351: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2352: PetscCall(DMGetDimension(cdm, &dim));
2353: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
2354: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
2355: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
2357: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
2358: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
2359: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2360: for (c = cStart; c < cEnd; ++c) {
2361: if (Npc == 1) {
2362: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
2363: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
2364: } else {
2365: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
2366: for (p = 0; p < Npc; ++p) {
2367: const PetscInt n = c * Npc + p;
2368: PetscReal sum = 0.0, refcoords[3];
2370: for (d = 0; d < dim; ++d) {
2371: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
2372: sum += refcoords[d];
2373: }
2374: if (simplex && sum > 0.0)
2375: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
2376: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
2377: }
2378: }
2379: }
2380: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2381: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
2382: PetscCall(PetscRandomDestroy(&rnd));
2383: PetscFunctionReturn(PETSC_SUCCESS);
2384: }
2386: /*@
2387: DMSwarmGetType - Get particular flavor of `DMSWARM`
2389: Collective
2391: Input Parameter:
2392: . sw - the `DMSWARM`
2394: Output Parameter:
2395: . stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2397: Level: advanced
2399: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2400: @*/
2401: PetscErrorCode DMSwarmGetType(DM sw, DMSwarmType *stype)
2402: {
2403: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2405: PetscFunctionBegin;
2407: PetscAssertPointer(stype, 2);
2408: *stype = swarm->swarm_type;
2409: PetscFunctionReturn(PETSC_SUCCESS);
2410: }
2412: /*@
2413: DMSwarmSetType - Set particular flavor of `DMSWARM`
2415: Collective
2417: Input Parameters:
2418: + sw - the `DMSWARM`
2419: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
2421: Level: advanced
2423: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2424: @*/
2425: PetscErrorCode DMSwarmSetType(DM sw, DMSwarmType stype)
2426: {
2427: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2429: PetscFunctionBegin;
2431: swarm->swarm_type = stype;
2432: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(sw));
2433: PetscFunctionReturn(PETSC_SUCCESS);
2434: }
2436: static PetscErrorCode DMSwarmCreateRemapDM_Private(DM sw, DM *rdm)
2437: {
2438: PetscFE fe;
2439: DMPolytopeType ct;
2440: PetscInt dim, cStart;
2441: const char *prefix = "remap_";
2443: PetscFunctionBegin;
2444: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), rdm));
2445: PetscCall(DMSetType(*rdm, DMPLEX));
2446: PetscCall(DMPlexSetOptionsPrefix(*rdm, prefix));
2447: PetscCall(DMSetFromOptions(*rdm));
2448: PetscCall(PetscObjectSetName((PetscObject)*rdm, "remap"));
2449: PetscCall(DMViewFromOptions(*rdm, NULL, "-dm_view"));
2451: PetscCall(DMGetDimension(*rdm, &dim));
2452: PetscCall(DMPlexGetHeightStratum(*rdm, 0, &cStart, NULL));
2453: PetscCall(DMPlexGetCellType(*rdm, cStart, &ct));
2454: PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
2455: PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
2456: PetscCall(DMSetField(*rdm, 0, NULL, (PetscObject)fe));
2457: PetscCall(DMCreateDS(*rdm));
2458: PetscCall(PetscFEDestroy(&fe));
2459: PetscFunctionReturn(PETSC_SUCCESS);
2460: }
2462: static PetscErrorCode DMSetup_Swarm(DM sw)
2463: {
2464: DM_Swarm *swarm = (DM_Swarm *)sw->data;
2466: PetscFunctionBegin;
2467: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
2468: swarm->issetup = PETSC_TRUE;
2470: if (swarm->remap_type != DMSWARM_REMAP_NONE) {
2471: DMSwarmCellDM celldm;
2473: PetscCall(DMSwarmGetCellDMByName(sw, "remap", &celldm));
2474: if (!celldm) {
2475: DM rdm;
2476: const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
2477: const char *vfieldnames[1] = {"w_q"};
2479: PetscCall(DMSwarmCreateRemapDM_Private(sw, &rdm));
2480: PetscCall(DMSwarmCellDMCreate(rdm, 1, vfieldnames, 2, fieldnames, &celldm));
2481: PetscCall(DMSwarmAddCellDM(sw, celldm));
2482: PetscCall(DMSwarmCellDMDestroy(&celldm));
2483: PetscCall(DMDestroy(&rdm));
2484: }
2485: }
2487: if (swarm->swarm_type == DMSWARM_PIC) {
2488: DMSwarmCellDM celldm;
2490: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2491: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "No active cell DM. DMSWARM_PIC requires you call DMSwarmSetCellDM() or DMSwarmAddCellDM()");
2492: if (celldm->dm->ops->locatepointssubdomain) {
2493: /* check methods exists for exact ownership identificiation */
2494: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
2495: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
2496: } else {
2497: /* check methods exist for point location AND rank neighbor identification */
2498: PetscCheck(celldm->dm->ops->locatepoints, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
2499: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
2501: PetscCheck(celldm->dm->ops->getneighbors, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
2502: PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
2504: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
2505: }
2506: }
2508: PetscCall(DMSwarmFinalizeFieldRegister(sw));
2510: /* check some fields were registered */
2511: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
2512: PetscFunctionReturn(PETSC_SUCCESS);
2513: }
2515: static PetscErrorCode DMDestroy_Swarm(DM dm)
2516: {
2517: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2519: PetscFunctionBegin;
2520: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2521: PetscCall(PetscObjectListDestroy(&swarm->cellDMs));
2522: PetscCall(PetscFree(swarm->activeCellDM));
2523: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
2524: PetscCall(PetscFree(swarm));
2525: PetscFunctionReturn(PETSC_SUCCESS);
2526: }
2528: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
2529: {
2530: DM cdm;
2531: DMSwarmCellDM celldm;
2532: PetscDraw draw;
2533: PetscReal *coords, oldPause, radius = 0.01;
2534: PetscInt Np, p, bs, Nfc;
2535: const char **coordFields;
2537: PetscFunctionBegin;
2538: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
2539: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2540: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2541: PetscCall(PetscDrawGetPause(draw, &oldPause));
2542: PetscCall(PetscDrawSetPause(draw, 0.0));
2543: PetscCall(DMView(cdm, viewer));
2544: PetscCall(PetscDrawSetPause(draw, oldPause));
2546: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2547: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2548: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2549: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2550: PetscCall(DMSwarmGetField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2551: for (p = 0; p < Np; ++p) {
2552: const PetscInt i = p * bs;
2554: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
2555: }
2556: PetscCall(DMSwarmRestoreField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2557: PetscCall(PetscDrawFlush(draw));
2558: PetscCall(PetscDrawPause(draw));
2559: PetscCall(PetscDrawSave(draw));
2560: PetscFunctionReturn(PETSC_SUCCESS);
2561: }
2563: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
2564: {
2565: PetscViewerFormat format;
2566: PetscDataType wtype = PETSC_DATATYPE_UNKNOWN;
2567: PetscInt *sizes;
2568: PetscInt dim, Np, wbs = 0, maxSize = 17;
2569: MPI_Comm comm;
2570: PetscMPIInt rank, size;
2571: const char *name, *cellid;
2573: PetscFunctionBegin;
2574: PetscCall(PetscViewerGetFormat(viewer, &format));
2575: PetscCall(DMGetDimension(dm, &dim));
2576: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2577: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2578: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2579: PetscCallMPI(MPI_Comm_size(comm, &size));
2580: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
2581: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
2582: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
2583: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
2584: else PetscCall(PetscCalloc1(3, &sizes));
2585: if (size < maxSize) {
2586: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
2587: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
2588: for (PetscInt p = 0; p < size; ++p) {
2589: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
2590: }
2591: } else {
2592: PetscInt locMinMax[2] = {Np, Np};
2594: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
2595: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
2596: }
2597: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2598: PetscCall(PetscFree(sizes));
2599: if (format == PETSC_VIEWER_ASCII_INFO) {
2600: DM_Swarm *sw = (DM_Swarm *)dm->data;
2601: DMSwarmCellDM celldm;
2602: PetscInt *cell;
2603: PetscBool hasWeight;
2604: const char *fname = "w_q";
2606: PetscCall(DMSwarmDataFieldStringInList(fname, sw->db->nfields, (const DMSwarmDataField *)sw->db->field, &hasWeight));
2607: if (hasWeight) {
2608: PetscCall(DMSwarmGetFieldInfo(dm, fname, &wbs, &wtype));
2609: 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]);
2610: PetscCheck(wbs == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must have block size 1, not %" PetscInt_FMT, fname, wbs);
2611: }
2612: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
2613: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2614: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2615: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
2616: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&cell));
2617: if (hasWeight) {
2618: PetscReal **coords;
2619: void *weight;
2620: PetscInt Ncf, *bsC;
2621: const char **coordNames;
2623: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncf, &coordNames));
2624: PetscCall(PetscMalloc2(Ncf, &coords, Ncf, &bsC));
2625: for (PetscInt n = 0; n < Ncf; ++n) PetscCall(DMSwarmGetField(dm, coordNames[n], &bsC[n], NULL, (void **)&coords[n]));
2626: PetscCall(DMSwarmGetField(dm, fname, NULL, NULL, &weight));
2627: for (PetscInt p = 0; p < Np; ++p) {
2628: PetscReal wp;
2630: if (wtype == PETSC_REAL) wp = ((PetscReal *)weight)[p];
2631: else wp = PetscRealPart(((PetscScalar *)weight)[p]);
2632: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%" PetscInt_FMT ": %" PetscInt_FMT " wt: %g x: (", p, cell[p], (double)wp));
2633: for (PetscInt n = 0; n < Ncf; ++n) {
2634: if (n > 0) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, ", "));
2635: for (PetscInt d = 0; d < bsC[n]; ++d) {
2636: if (d > 0) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, ", "));
2637: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%g", (double)coords[n][p * bsC[n] + d]));
2638: }
2639: }
2640: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, ")\n"));
2641: }
2642: PetscCall(DMSwarmRestoreField(dm, fname, NULL, NULL, &weight));
2643: for (PetscInt n = 0; n < Ncf; ++n) PetscCall(DMSwarmRestoreField(dm, coordNames[n], &bsC[n], NULL, (void **)&coords[n]));
2644: PetscCall(PetscFree2(coords, bsC));
2645: } else {
2646: for (PetscInt p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%" PetscInt_FMT ": %" PetscInt_FMT "\n", p, cell[p]));
2647: }
2648: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&cell));
2649: PetscCall(PetscViewerFlush(viewer));
2650: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2651: }
2652: PetscFunctionReturn(PETSC_SUCCESS);
2653: }
2655: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
2656: {
2657: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2658: PetscBool isascii, ibinary, isvtk, isdraw, ispython;
2659: #if PetscDefined(HAVE_HDF5)
2660: PetscBool ishdf5;
2661: #endif
2663: PetscFunctionBegin;
2666: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2667: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
2668: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
2669: #if PetscDefined(HAVE_HDF5)
2670: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2671: #endif
2672: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2673: PetscCall(PetscObjectHasFunction((PetscObject)viewer, "PetscViewerPythonViewObject_C", &ispython));
2674: if (isascii) {
2675: PetscViewerFormat format;
2677: PetscCall(PetscViewerGetFormat(viewer, &format));
2678: switch (format) {
2679: case PETSC_VIEWER_ASCII_INFO_DETAIL:
2680: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
2681: break;
2682: default:
2683: PetscCall(DMView_Swarm_Ascii(dm, viewer));
2684: }
2685: } else {
2686: #if PetscDefined(HAVE_HDF5)
2687: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
2688: #endif
2689: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
2690: if (ispython) PetscCall(PetscViewerPythonViewObject(viewer, (PetscObject)dm));
2691: }
2692: PetscFunctionReturn(PETSC_SUCCESS);
2693: }
2695: /*@
2696: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
2697: 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.
2699: Noncollective
2701: Input Parameters:
2702: + sw - the `DMSWARM`
2703: . cellID - the integer id of the cell to be extracted and filtered
2704: - cellswarm - The `DMSWARM` to receive the cell
2706: Level: beginner
2708: Notes:
2709: This presently only supports `DMSWARM_PIC` type
2711: Should be restored with `DMSwarmRestoreCellSwarm()`
2713: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
2715: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
2716: @*/
2717: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2718: {
2719: DM_Swarm *original = (DM_Swarm *)sw->data;
2720: DMLabel label;
2721: DM dmc, subdmc;
2722: PetscInt *pids, particles, dim;
2723: const char *name;
2725: PetscFunctionBegin;
2726: /* Configure new swarm */
2727: PetscCall(DMSetType(cellswarm, DMSWARM));
2728: PetscCall(DMGetDimension(sw, &dim));
2729: PetscCall(DMSetDimension(cellswarm, dim));
2730: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
2731: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
2732: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
2733: PetscCall(DMSwarmSortGetAccess(sw));
2734: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
2735: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2736: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
2737: PetscCall(DMSwarmSortRestoreAccess(sw));
2738: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2739: PetscCall(DMSwarmGetCellDM(sw, &dmc));
2740: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
2741: PetscCall(DMAddLabel(dmc, label));
2742: PetscCall(DMLabelSetValue(label, cellID, 1));
2743: PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, PetscObjectComm((PetscObject)dmc), NULL, &subdmc));
2744: PetscCall(PetscObjectGetName((PetscObject)dmc, &name));
2745: PetscCall(PetscObjectSetName((PetscObject)subdmc, name));
2746: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
2747: PetscCall(DMLabelDestroy(&label));
2748: PetscFunctionReturn(PETSC_SUCCESS);
2749: }
2751: /*@
2752: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
2754: Noncollective
2756: Input Parameters:
2757: + sw - the parent `DMSWARM`
2758: . cellID - the integer id of the cell to be copied back into the parent swarm
2759: - cellswarm - the cell swarm object
2761: Level: beginner
2763: Note:
2764: This only supports `DMSWARM_PIC` types of `DMSWARM`s
2766: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
2767: @*/
2768: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2769: {
2770: DM dmc;
2771: PetscInt *pids, particles, p;
2773: PetscFunctionBegin;
2774: PetscCall(DMSwarmSortGetAccess(sw));
2775: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2776: PetscCall(DMSwarmSortRestoreAccess(sw));
2777: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
2778: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
2779: /* Free memory, destroy cell dm */
2780: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
2781: PetscCall(DMDestroy(&dmc));
2782: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2783: PetscFunctionReturn(PETSC_SUCCESS);
2784: }
2786: /*@
2787: DMSwarmComputeMoments - Compute the first three particle moments for a given field
2789: Noncollective
2791: Input Parameters:
2792: + sw - the `DMSWARM`
2793: . coordinate - the coordinate field name
2794: - weight - the weight field name
2796: Output Parameter:
2797: . moments - the field moments
2799: Level: intermediate
2801: Notes:
2802: The `moments` array should be of length bs + 2, where bs is the block size of the coordinate field.
2804: The weight field must have blocksize 1 and data type `PETSC_REAL` or `PETSC_SCALAR`. For complex scalars,
2805: only the real part of each weight is used.
2807: .seealso: `DM`, `DMSWARM`, `DMPlexComputeMoments()`
2808: @*/
2809: PetscErrorCode DMSwarmComputeMoments(DM sw, const char coordinate[], const char weight[], PetscReal moments[])
2810: {
2811: const PetscReal *coords;
2812: void *weights;
2813: PetscDataType dtc, dtw;
2814: PetscInt bsc, bsw, Np;
2815: MPI_Comm comm;
2817: PetscFunctionBegin;
2819: PetscAssertPointer(coordinate, 2);
2820: PetscAssertPointer(weight, 3);
2821: PetscAssertPointer(moments, 4);
2822: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
2823: PetscCall(DMSwarmGetFieldInfo(sw, coordinate, &bsc, &dtc));
2824: PetscCall(DMSwarmGetFieldInfo(sw, weight, &bsw, &dtw));
2825: PetscCheck(dtc == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s must be real, not %s", coordinate, PetscDataTypes[dtc]);
2826: 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]);
2827: PetscCheck(bsw == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be a scalar, not blocksize %" PetscInt_FMT, weight, bsw);
2828: PetscCall(DMSwarmGetField(sw, coordinate, NULL, NULL, (void **)&coords));
2829: PetscCall(DMSwarmGetField(sw, weight, NULL, NULL, &weights));
2830: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2831: PetscCall(PetscArrayzero(moments, bsc + 2));
2832: for (PetscInt p = 0; p < Np; ++p) {
2833: const PetscReal *c = &coords[p * bsc];
2834: PetscReal wp;
2836: if (dtw == PETSC_REAL) wp = ((PetscReal *)weights)[p];
2837: else wp = PetscRealPart(((PetscScalar *)weights)[p]);
2838: moments[0] += wp;
2839: for (PetscInt d = 0; d < bsc; ++d) {
2840: moments[d + 1] += wp * c[d];
2841: moments[bsc + 1] += wp * PetscSqr(c[d]);
2842: }
2843: }
2844: PetscCall(DMSwarmRestoreField(sw, coordinate, NULL, NULL, (void **)&coords));
2845: PetscCall(DMSwarmRestoreField(sw, weight, NULL, NULL, &weights));
2846: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, moments, bsc + 2, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
2847: PetscFunctionReturn(PETSC_SUCCESS);
2848: }
2850: static PetscErrorCode DMSetFromOptions_Swarm(DM dm, PetscOptionItems PetscOptionsObject)
2851: {
2852: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2854: PetscFunctionBegin;
2855: PetscOptionsHeadBegin(PetscOptionsObject, "DMSwarm Options");
2856: PetscCall(PetscOptionsEnum("-dm_swarm_remap_type", "Remap algorithm", "DMSwarmSetRemapType", DMSwarmRemapTypeNames, (PetscEnum)swarm->remap_type, (PetscEnum *)&swarm->remap_type, NULL));
2857: PetscOptionsHeadEnd();
2858: PetscFunctionReturn(PETSC_SUCCESS);
2859: }
2861: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
2863: static PetscErrorCode DMInitialize_Swarm(DM sw)
2864: {
2865: PetscFunctionBegin;
2866: sw->ops->view = DMView_Swarm;
2867: sw->ops->load = NULL;
2868: sw->ops->setfromoptions = DMSetFromOptions_Swarm;
2869: sw->ops->clone = DMClone_Swarm;
2870: sw->ops->setup = DMSetup_Swarm;
2871: sw->ops->createlocalsection = NULL;
2872: sw->ops->createsectionpermutation = NULL;
2873: sw->ops->createdefaultconstraints = NULL;
2874: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
2875: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
2876: sw->ops->getlocaltoglobalmapping = NULL;
2877: sw->ops->createfieldis = NULL;
2878: sw->ops->createcoordinatedm = NULL;
2879: sw->ops->createcellcoordinatedm = NULL;
2880: sw->ops->getcoloring = NULL;
2881: sw->ops->creatematrix = DMCreateMatrix_Swarm;
2882: sw->ops->createinterpolation = NULL;
2883: sw->ops->createinjection = NULL;
2884: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
2885: sw->ops->creategradientmatrix = DMCreateGradientMatrix_Swarm;
2886: sw->ops->refine = NULL;
2887: sw->ops->coarsen = NULL;
2888: sw->ops->refinehierarchy = NULL;
2889: sw->ops->coarsenhierarchy = NULL;
2890: sw->ops->globaltolocalbegin = DMGlobalToLocalBegin_Swarm;
2891: sw->ops->globaltolocalend = DMGlobalToLocalEnd_Swarm;
2892: sw->ops->localtoglobalbegin = DMLocalToGlobalBegin_Swarm;
2893: sw->ops->localtoglobalend = DMLocalToGlobalEnd_Swarm;
2894: sw->ops->destroy = DMDestroy_Swarm;
2895: sw->ops->createsubdm = NULL;
2896: sw->ops->getdimpoints = NULL;
2897: sw->ops->locatepoints = NULL;
2898: sw->ops->projectfieldlocal = DMProjectFieldLocal_Swarm;
2899: PetscFunctionReturn(PETSC_SUCCESS);
2900: }
2902: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
2903: {
2904: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2906: PetscFunctionBegin;
2907: swarm->refct++;
2908: (*newdm)->data = swarm;
2909: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
2910: PetscCall(DMInitialize_Swarm(*newdm));
2911: (*newdm)->dim = dm->dim;
2912: PetscFunctionReturn(PETSC_SUCCESS);
2913: }
2915: /*MC
2916: DMSWARM = "swarm" - A `DM` object for particle methods, such as particle-in-cell (PIC), in which the underlying
2917: data is both (i) dynamic in length, (ii) and of arbitrary data type.
2919: Level: intermediate
2921: Notes:
2922: User data can be represented by `DMSWARM` through a registering "fields" which are to be stored on particles.
2923: To register a field, the user must provide;
2924: (a) a unique name;
2925: (b) the data type (or size in bytes);
2926: (c) the block size of the data.
2928: For example, suppose the application requires a unique id, energy, momentum and density to be stored
2929: on a set of particles. Then the following code could be used
2930: .vb
2931: DMSwarmInitializeFieldRegister(dm)
2932: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
2933: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
2934: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
2935: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
2936: DMSwarmFinalizeFieldRegister(dm)
2937: .ve
2939: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
2940: The only restriction imposed by `DMSWARM` is that all fields contain the same number of particles.
2942: To support particle methods, "migration" techniques are provided. These methods migrate data
2943: between ranks.
2945: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
2946: As a `DMSWARM` may internally define and store values of different data types,
2947: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
2948: fields should be used to define a `Vec` object via `DMSwarmVectorDefineField()`
2949: The specified field can be changed at any time - thereby permitting vectors
2950: compatible with different fields to be created.
2952: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via `DMSwarmCreateGlobalVectorFromField()`
2953: Here the data defining the field in the `DMSWARM` is shared with a `Vec`.
2954: This is inherently unsafe if you alter the size of the field at any time between
2955: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
2956: If the local size of the `DMSWARM` does not match the local size of the global vector
2957: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
2959: Additional high-level support is provided for Particle-In-Cell methods. Refer to `DMSwarmSetType()`.
2961: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`, `DMSwarmSetType()`, `DMSwarmType`, `DMSwarmCreateGlobalVectorFromField()`,
2962: `DMCreateGlobalVector()`, `DMCreateLocalVector()`
2963: M*/
2965: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
2966: {
2967: DM_Swarm *swarm;
2969: PetscFunctionBegin;
2971: PetscCall(PetscNew(&swarm));
2972: dm->data = swarm;
2973: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
2974: PetscCall(DMSwarmInitializeFieldRegister(dm));
2975: dm->dim = 0;
2976: swarm->refct = 1;
2977: swarm->issetup = PETSC_FALSE;
2978: swarm->swarm_type = DMSWARM_BASIC;
2979: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
2980: swarm->collect_type = DMSWARM_COLLECT_BASIC;
2981: swarm->migrate_error_on_missing_point = PETSC_FALSE;
2982: swarm->collect_view_active = PETSC_FALSE;
2983: swarm->collect_view_reset_nlocal = -1;
2984: PetscCall(DMInitialize_Swarm(dm));
2985: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
2986: PetscFunctionReturn(PETSC_SUCCESS);
2987: }
2989: /* Replace dm with the contents of ndm, and then destroy ndm
2990: - Share the DM_Swarm structure
2991: */
2992: /*@
2993: DMSwarmReplace - Replace the internal state of a `DMSWARM` with that of another `DMSWARM`, sharing the underlying particle data and destroying the source `DM`.
2995: Collective
2997: Input Parameters:
2998: + dm - the destination `DMSWARM`, whose current contents are discarded
2999: - ndm - pointer to the source `DMSWARM`; destroyed and set to `NULL` on return
3001: Level: developer
3003: Note:
3004: The dimension, periodicity, name, and shared reference to the underlying particle bucket are transferred from the source to the destination.
3006: .seealso: `DM`, `DMSWARM`, `DMSwarmDuplicate()`, `DMDestroy()`
3007: @*/
3008: PetscErrorCode DMSwarmReplace(DM dm, DM *ndm)
3009: {
3010: DM dmNew = *ndm;
3011: const PetscReal *maxCell, *Lstart, *L;
3012: PetscInt dim;
3014: PetscFunctionBegin;
3015: if (dm == dmNew) {
3016: PetscCall(DMDestroy(ndm));
3017: PetscFunctionReturn(PETSC_SUCCESS);
3018: }
3019: dm->setupcalled = dmNew->setupcalled;
3020: if (!dm->hdr.name) {
3021: const char *name;
3023: PetscCall(PetscObjectGetName((PetscObject)*ndm, &name));
3024: PetscCall(PetscObjectSetName((PetscObject)dm, name));
3025: }
3026: PetscCall(DMGetDimension(dmNew, &dim));
3027: PetscCall(DMSetDimension(dm, dim));
3028: PetscCall(DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L));
3029: PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
3030: PetscCall(DMDestroy_Swarm(dm));
3031: PetscCall(DMInitialize_Swarm(dm));
3032: dm->data = dmNew->data;
3033: ((DM_Swarm *)dmNew->data)->refct++;
3034: PetscCall(DMDestroy(ndm));
3035: PetscFunctionReturn(PETSC_SUCCESS);
3036: }
3038: /*@
3039: DMSwarmDuplicate - Creates a new `DMSWARM` with the same fields and cell `DM`s but no particles
3041: Collective
3043: Input Parameter:
3044: . sw - the `DMSWARM`
3046: Output Parameter:
3047: . nsw - the new `DMSWARM`
3049: Level: beginner
3051: .seealso: `DM`, `DMSWARM`, `DMSwarmCreate()`, `DMClone()`
3052: @*/
3053: PetscErrorCode DMSwarmDuplicate(DM sw, DM *nsw)
3054: {
3055: DM_Swarm *swarm = (DM_Swarm *)sw->data;
3056: DMSwarmDataField *fields;
3057: DMSwarmCellDM celldm, ncelldm;
3058: DMSwarmType stype;
3059: const char *name, **celldmnames;
3060: void *ctx;
3061: PetscInt dim, Nf, Ndm;
3062: PetscBool flg;
3064: PetscFunctionBegin;
3065: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), nsw));
3066: PetscCall(DMSetType(*nsw, DMSWARM));
3067: PetscCall(PetscObjectGetName((PetscObject)sw, &name));
3068: PetscCall(PetscObjectSetName((PetscObject)*nsw, name));
3069: PetscCall(DMGetDimension(sw, &dim));
3070: PetscCall(DMSetDimension(*nsw, dim));
3071: PetscCall(DMSwarmGetType(sw, &stype));
3072: PetscCall(DMSwarmSetType(*nsw, stype));
3073: PetscCall(DMGetApplicationContext(sw, &ctx));
3074: PetscCall(DMSetApplicationContext(*nsw, ctx));
3076: PetscCall(DMSwarmDataBucketGetDMSwarmDataFields(swarm->db, &Nf, &fields));
3077: for (PetscInt f = 0; f < Nf; ++f) {
3078: PetscCall(DMSwarmDataFieldStringInList(fields[f]->name, ((DM_Swarm *)(*nsw)->data)->db->nfields, (const DMSwarmDataField *)((DM_Swarm *)(*nsw)->data)->db->field, &flg));
3079: if (!flg) PetscCall(DMSwarmRegisterPetscDatatypeField(*nsw, fields[f]->name, fields[f]->bs, fields[f]->petsc_type));
3080: }
3082: PetscCall(DMSwarmGetCellDMNames(sw, &Ndm, &celldmnames));
3083: for (PetscInt c = 0; c < Ndm; ++c) {
3084: DM dm;
3085: PetscInt Ncf;
3086: const char **coordfields, **fields;
3088: PetscCall(DMSwarmGetCellDMByName(sw, celldmnames[c], &celldm));
3089: PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
3090: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncf, &coordfields));
3091: PetscCall(DMSwarmCellDMGetFields(celldm, &Nf, &fields));
3092: PetscCall(DMSwarmCellDMCreate(dm, Nf, fields, Ncf, coordfields, &ncelldm));
3093: PetscCall(DMSwarmAddCellDM(*nsw, ncelldm));
3094: PetscCall(DMSwarmCellDMDestroy(&ncelldm));
3095: }
3096: PetscCall(PetscFree(celldmnames));
3098: PetscCall(DMSetFromOptions(*nsw));
3099: PetscCall(DMSetUp(*nsw));
3100: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
3101: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
3102: PetscCall(DMSwarmSetCellDMActive(*nsw, name));
3103: PetscFunctionReturn(PETSC_SUCCESS);
3104: }
3106: PetscErrorCode DMLocalToGlobalBegin_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
3107: {
3108: PetscFunctionBegin;
3109: PetscFunctionReturn(PETSC_SUCCESS);
3110: }
3112: PetscErrorCode DMLocalToGlobalEnd_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
3113: {
3114: PetscFunctionBegin;
3115: switch (mode) {
3116: case INSERT_VALUES:
3117: PetscCall(VecCopy(l, g));
3118: break;
3119: case ADD_VALUES:
3120: PetscCall(VecAXPY(g, 1., l));
3121: break;
3122: default:
3123: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mode not supported: %d", mode);
3124: }
3125: PetscFunctionReturn(PETSC_SUCCESS);
3126: }
3128: PetscErrorCode DMGlobalToLocalBegin_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
3129: {
3130: PetscFunctionBegin;
3131: PetscFunctionReturn(PETSC_SUCCESS);
3132: }
3134: PetscErrorCode DMGlobalToLocalEnd_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
3135: {
3136: PetscFunctionBegin;
3137: PetscCall(DMLocalToGlobalEnd_Swarm(dm, g, mode, l));
3138: PetscFunctionReturn(PETSC_SUCCESS);
3139: }