Actual source code: swarm.c
1: #include "petscdm.h"
2: #include "petscdmswarm.h"
3: #include "petscstring.h"
4: #include "petscsys.h"
5: #define PETSCDM_DLL
6: #include <petsc/private/dmswarmimpl.h>
7: #include <petsc/private/hashsetij.h>
8: #include <petsc/private/petscfeimpl.h>
9: #include <petscviewer.h>
10: #include <petscdraw.h>
11: #include <petscdmplex.h>
12: #include <petscblaslapack.h>
13: #include "../src/dm/impls/swarm/data_bucket.h"
14: #include <petscdmlabel.h>
15: #include <petscsection.h>
17: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
18: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
19: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
21: const char *DMSwarmTypeNames[] = {"basic", "pic", NULL};
22: const char *DMSwarmMigrateTypeNames[] = {"basic", "dmcellnscatter", "dmcellexact", "user", NULL};
23: const char *DMSwarmCollectTypeNames[] = {"basic", "boundingbox", "general", "user", NULL};
24: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
26: const char DMSwarmField_pid[] = "DMSwarm_pid";
27: const char DMSwarmField_rank[] = "DMSwarm_rank";
28: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
30: PetscInt SwarmDataFieldId = -1;
32: #if defined(PETSC_HAVE_HDF5)
33: #include <petscviewerhdf5.h>
35: static PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
36: {
37: DM dm;
38: PetscReal seqval;
39: PetscInt seqnum, bs;
40: PetscBool isseq, ists;
42: PetscFunctionBegin;
43: PetscCall(VecGetDM(v, &dm));
44: PetscCall(VecGetBlockSize(v, &bs));
45: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particle_fields"));
46: PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
47: PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
48: PetscCall(PetscViewerHDF5IsTimestepping(viewer, &ists));
49: if (ists) PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
50: /* PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer)); */
51: PetscCall(VecViewNative(v, viewer));
52: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs));
53: PetscCall(PetscViewerHDF5PopGroup(viewer));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
58: {
59: DMSwarmCellDM celldm;
60: Vec coordinates;
61: PetscInt Np, Nfc;
62: PetscBool isseq;
63: const char **coordFields;
65: PetscFunctionBegin;
66: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
67: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
68: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
69: PetscCall(DMSwarmGetSize(dm, &Np));
70: PetscCall(DMSwarmCreateGlobalVectorFromField(dm, coordFields[0], &coordinates));
71: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
72: PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
73: PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
74: PetscCall(VecViewNative(coordinates, viewer));
75: PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
76: PetscCall(PetscViewerHDF5PopGroup(viewer));
77: PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, coordFields[0], &coordinates));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
80: #endif
82: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
83: {
84: DM dm;
85: #if defined(PETSC_HAVE_HDF5)
86: PetscBool ishdf5;
87: #endif
89: PetscFunctionBegin;
90: PetscCall(VecGetDM(v, &dm));
91: PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
92: #if defined(PETSC_HAVE_HDF5)
93: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
94: if (ishdf5) {
95: PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
98: #endif
99: PetscCall(VecViewNative(v, viewer));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: /*@C
104: DMSwarmVectorGetField - Gets the fields from which to define a `Vec` object
105: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
107: Not collective
109: Input Parameter:
110: . sw - a `DMSWARM`
112: Output Parameters:
113: + Nf - the number of fields
114: - fieldnames - the textual name given to each registered field, or NULL if it has not been set
116: Level: beginner
118: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
119: @*/
120: PetscErrorCode DMSwarmVectorGetField(DM sw, PetscInt *Nf, const char **fieldnames[])
121: {
122: DMSwarmCellDM celldm;
124: PetscFunctionBegin;
126: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
127: PetscCall(DMSwarmCellDMGetFields(celldm, Nf, fieldnames));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
133: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
135: Collective
137: Input Parameters:
138: + dm - a `DMSWARM`
139: - fieldname - the textual name given to each registered field
141: Level: beginner
143: Notes:
144: The field with name `fieldname` must be defined as having a data type of `PetscScalar`.
146: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
147: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
149: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineFields()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
150: @*/
151: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
152: {
153: PetscFunctionBegin;
154: PetscCall(DMSwarmVectorDefineFields(dm, 1, &fieldname));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: /*@C
159: DMSwarmVectorDefineFields - Sets the fields from which to define a `Vec` object
160: when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called
162: Collective, No Fortran support
164: Input Parameters:
165: + sw - a `DMSWARM`
166: . Nf - the number of fields
167: - fieldnames - the textual name given to each registered field
169: Level: beginner
171: Notes:
172: Each field with name in `fieldnames` must be defined as having a data type of `PetscScalar`.
174: This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
175: Multiple calls to `DMSwarmVectorDefineField()` are permitted.
177: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
178: @*/
179: PetscErrorCode DMSwarmVectorDefineFields(DM sw, PetscInt Nf, const char *fieldnames[])
180: {
181: DM_Swarm *swarm = (DM_Swarm *)sw->data;
182: DMSwarmCellDM celldm;
184: PetscFunctionBegin;
186: if (fieldnames) PetscAssertPointer(fieldnames, 3);
187: if (!swarm->issetup) PetscCall(DMSetUp(sw));
188: PetscCheck(Nf >= 0, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields must be non-negative, not %" PetscInt_FMT, Nf);
189: // Create a dummy cell DM if none has been specified (I think we should not support this mode)
190: if (!swarm->activeCellDM) {
191: DM dm;
192: DMSwarmCellDM celldm;
194: PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), &dm));
195: PetscCall(DMSetType(dm, DMSHELL));
196: PetscCall(PetscObjectSetName((PetscObject)dm, "dummy"));
197: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 0, NULL, &celldm));
198: PetscCall(DMDestroy(&dm));
199: PetscCall(DMSwarmAddCellDM(sw, celldm));
200: PetscCall(DMSwarmCellDMDestroy(&celldm));
201: PetscCall(DMSwarmSetCellDMActive(sw, "dummy"));
202: }
203: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
204: for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscFree(celldm->dmFields[f]));
205: PetscCall(PetscFree(celldm->dmFields));
207: celldm->Nf = Nf;
208: PetscCall(PetscMalloc1(Nf, &celldm->dmFields));
209: for (PetscInt f = 0; f < Nf; ++f) {
210: PetscDataType type;
212: // Check all fields are of type PETSC_REAL or PETSC_SCALAR
213: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], NULL, &type));
214: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
215: PetscCall(PetscStrallocpy(fieldnames[f], (char **)&celldm->dmFields[f]));
216: }
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /* requires DMSwarmDefineFieldVector has been called */
221: static PetscErrorCode DMCreateGlobalVector_Swarm(DM sw, Vec *vec)
222: {
223: DM_Swarm *swarm = (DM_Swarm *)sw->data;
224: DMSwarmCellDM celldm;
225: Vec x;
226: char name[PETSC_MAX_PATH_LEN];
227: PetscInt bs = 0, n;
229: PetscFunctionBegin;
230: if (!swarm->issetup) PetscCall(DMSetUp(sw));
231: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
232: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
233: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
235: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
236: for (PetscInt f = 0; f < celldm->Nf; ++f) {
237: PetscInt fbs;
238: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
239: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
240: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
241: bs += fbs;
242: }
243: PetscCall(VecCreate(PetscObjectComm((PetscObject)sw), &x));
244: PetscCall(PetscObjectSetName((PetscObject)x, name));
245: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
246: PetscCall(VecSetBlockSize(x, bs));
247: PetscCall(VecSetDM(x, sw));
248: PetscCall(VecSetFromOptions(x));
249: PetscCall(VecSetOperation(x, VECOP_VIEW, (void (*)(void))VecView_Swarm));
250: *vec = x;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: /* requires DMSwarmDefineFieldVector has been called */
255: static PetscErrorCode DMCreateLocalVector_Swarm(DM sw, Vec *vec)
256: {
257: DM_Swarm *swarm = (DM_Swarm *)sw->data;
258: DMSwarmCellDM celldm;
259: Vec x;
260: char name[PETSC_MAX_PATH_LEN];
261: PetscInt bs = 0, n;
263: PetscFunctionBegin;
264: if (!swarm->issetup) PetscCall(DMSetUp(sw));
265: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
266: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
267: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
269: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
270: for (PetscInt f = 0; f < celldm->Nf; ++f) {
271: PetscInt fbs;
272: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
273: PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
274: PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
275: bs += fbs;
276: }
277: PetscCall(VecCreate(PETSC_COMM_SELF, &x));
278: PetscCall(PetscObjectSetName((PetscObject)x, name));
279: PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
280: PetscCall(VecSetBlockSize(x, bs));
281: PetscCall(VecSetDM(x, sw));
282: PetscCall(VecSetFromOptions(x));
283: *vec = x;
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
288: {
289: DM_Swarm *swarm = (DM_Swarm *)dm->data;
290: DMSwarmDataField gfield;
291: PetscInt bs, nlocal, fid = -1, cfid = -2;
292: PetscBool flg;
294: PetscFunctionBegin;
295: /* check vector is an inplace array */
296: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
297: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
298: (void)flg; /* avoid compiler warning */
299: 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);
300: PetscCall(VecGetLocalSize(*vec, &nlocal));
301: PetscCall(VecGetBlockSize(*vec, &bs));
302: 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");
303: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
304: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
305: PetscCall(VecResetArray(*vec));
306: PetscCall(VecDestroy(vec));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
311: {
312: DM_Swarm *swarm = (DM_Swarm *)dm->data;
313: PetscDataType type;
314: PetscScalar *array;
315: PetscInt bs, n, fid;
316: char name[PETSC_MAX_PATH_LEN];
317: PetscMPIInt size;
318: PetscBool iscuda, iskokkos, iship;
320: PetscFunctionBegin;
321: if (!swarm->issetup) PetscCall(DMSetUp(dm));
322: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
323: PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
324: PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
326: PetscCallMPI(MPI_Comm_size(comm, &size));
327: PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
328: PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
329: PetscCall(PetscStrcmp(dm->vectype, VECHIP, &iship));
330: PetscCall(VecCreate(comm, vec));
331: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
332: PetscCall(VecSetBlockSize(*vec, bs));
333: if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
334: else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
335: else if (iship) PetscCall(VecSetType(*vec, VECHIP));
336: else PetscCall(VecSetType(*vec, VECSTANDARD));
337: PetscCall(VecPlaceArray(*vec, array));
339: PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
340: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
342: /* Set guard */
343: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
344: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));
346: PetscCall(VecSetDM(*vec, dm));
347: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: static PetscErrorCode DMSwarmDestroyVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], Vec *vec)
352: {
353: DM_Swarm *swarm = (DM_Swarm *)sw->data;
354: const PetscScalar *array;
355: PetscInt bs, n, id = 0, cid = -2;
356: PetscBool flg;
358: PetscFunctionBegin;
359: for (PetscInt f = 0; f < Nf; ++f) {
360: PetscInt fid;
362: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
363: id += fid;
364: }
365: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cid, flg));
366: (void)flg; /* avoid compiler warning */
367: 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);
368: PetscCall(VecGetLocalSize(*vec, &n));
369: PetscCall(VecGetBlockSize(*vec, &bs));
370: n /= bs;
371: PetscCheck(n == swarm->db->L, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
372: PetscCall(VecGetArrayRead(*vec, &array));
373: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
374: PetscScalar *farray;
375: PetscDataType ftype;
376: PetscInt fbs;
378: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
379: PetscCheck(off + fbs <= bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid blocksize %" PetscInt_FMT " < %" PetscInt_FMT, bs, off + fbs);
380: for (PetscInt i = 0; i < n; ++i) {
381: for (PetscInt b = 0; b < fbs; ++b) farray[i * fbs + b] = array[i * bs + off + b];
382: }
383: off += fbs;
384: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
385: }
386: PetscCall(VecRestoreArrayRead(*vec, &array));
387: PetscCall(VecDestroy(vec));
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: static PetscErrorCode DMSwarmCreateVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], MPI_Comm comm, Vec *vec)
392: {
393: DM_Swarm *swarm = (DM_Swarm *)sw->data;
394: PetscScalar *array;
395: PetscInt n, bs = 0, id = 0;
396: char name[PETSC_MAX_PATH_LEN];
398: PetscFunctionBegin;
399: if (!swarm->issetup) PetscCall(DMSetUp(sw));
400: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
401: for (PetscInt f = 0; f < Nf; ++f) {
402: PetscDataType ftype;
403: PetscInt fbs;
405: PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
406: PetscCheck(ftype == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
407: bs += fbs;
408: }
410: PetscCall(VecCreate(comm, vec));
411: PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
412: PetscCall(VecSetBlockSize(*vec, bs));
413: PetscCall(VecSetType(*vec, sw->vectype));
415: PetscCall(VecGetArrayWrite(*vec, &array));
416: for (PetscInt f = 0, off = 0; f < Nf; ++f) {
417: PetscScalar *farray;
418: PetscDataType ftype;
419: PetscInt fbs;
421: PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
422: for (PetscInt i = 0; i < n; ++i) {
423: for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = farray[i * fbs + b];
424: }
425: off += fbs;
426: PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
427: }
428: PetscCall(VecRestoreArrayWrite(*vec, &array));
430: PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
431: for (PetscInt f = 0; f < Nf; ++f) {
432: PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
433: PetscCall(PetscStrlcat(name, fieldnames[f], PETSC_MAX_PATH_LEN));
434: }
435: PetscCall(PetscObjectSetName((PetscObject)*vec, name));
437: for (PetscInt f = 0; f < Nf; ++f) {
438: PetscInt fid;
440: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
441: id += fid;
442: }
443: PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, id));
445: PetscCall(VecSetDM(*vec, sw));
446: PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm));
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
452: \hat f = \sum_i f_i \phi_i
454: and a particle function is given by
456: f = \sum_i w_i \delta(x - x_i)
458: then we want to require that
460: M \hat f = M_p f
462: where the particle mass matrix is given by
464: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
466: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
467: his integral. We allow this with the boolean flag.
468: */
469: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
470: {
471: const char *name = "Mass Matrix";
472: MPI_Comm comm;
473: DMSwarmCellDM celldm;
474: PetscDS prob;
475: PetscSection fsection, globalFSection;
476: PetscHSetIJ ht;
477: PetscLayout rLayout, colLayout;
478: PetscInt *dnz, *onz;
479: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
480: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
481: PetscScalar *elemMat;
482: PetscInt dim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0, totNc = 0;
483: const char **coordFields;
484: PetscReal **coordVals;
485: PetscInt *bs;
487: PetscFunctionBegin;
488: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
489: PetscCall(DMGetCoordinateDim(dmf, &dim));
490: PetscCall(DMGetDS(dmf, &prob));
491: PetscCall(PetscDSGetNumFields(prob, &Nf));
492: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
493: PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
494: PetscCall(DMGetLocalSection(dmf, &fsection));
495: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
496: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
497: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
499: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
500: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
501: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
503: PetscCall(PetscLayoutCreate(comm, &colLayout));
504: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
505: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
506: PetscCall(PetscLayoutSetUp(colLayout));
507: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
508: PetscCall(PetscLayoutDestroy(&colLayout));
510: PetscCall(PetscLayoutCreate(comm, &rLayout));
511: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
512: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
513: PetscCall(PetscLayoutSetUp(rLayout));
514: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
515: PetscCall(PetscLayoutDestroy(&rLayout));
517: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
518: PetscCall(PetscHSetIJCreate(&ht));
520: PetscCall(PetscSynchronizedFlush(comm, NULL));
521: for (PetscInt field = 0; field < Nf; ++field) {
522: PetscObject obj;
523: PetscClassId id;
524: PetscInt Nc;
526: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
527: PetscCall(PetscObjectGetClassId(obj, &id));
528: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
529: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
530: totNc += Nc;
531: }
532: /* count non-zeros */
533: PetscCall(DMSwarmSortGetAccess(dmc));
534: for (PetscInt field = 0; field < Nf; ++field) {
535: PetscObject obj;
536: PetscClassId id;
537: PetscInt Nc;
539: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
540: PetscCall(PetscObjectGetClassId(obj, &id));
541: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
542: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
544: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
545: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
546: PetscInt numFIndices, numCIndices;
548: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
549: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
550: maxC = PetscMax(maxC, numCIndices);
551: {
552: PetscHashIJKey key;
553: PetscBool missing;
554: for (PetscInt i = 0; i < numFIndices; ++i) {
555: key.j = findices[i]; /* global column (from Plex) */
556: if (key.j >= 0) {
557: /* Get indices for coarse elements */
558: for (PetscInt j = 0; j < numCIndices; ++j) {
559: for (PetscInt c = 0; c < Nc; ++c) {
560: // TODO Need field offset on particle here
561: key.i = cindices[j] * totNc + c + rStart; /* global cols (from Swarm) */
562: if (key.i < 0) continue;
563: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
564: if (missing) {
565: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
566: else ++onz[key.i - rStart];
567: } else SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
568: }
569: }
570: }
571: }
572: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
573: }
574: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
575: }
576: }
577: PetscCall(PetscHSetIJDestroy(&ht));
578: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
579: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
580: PetscCall(PetscFree2(dnz, onz));
581: PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
582: for (PetscInt field = 0; field < Nf; ++field) {
583: PetscTabulation Tcoarse;
584: PetscObject obj;
585: PetscClassId id;
586: PetscInt Nc;
588: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
589: PetscCall(PetscObjectGetClassId(obj, &id));
590: if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
591: else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
593: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
594: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
595: PetscInt *findices, *cindices;
596: PetscInt numFIndices, numCIndices;
598: /* TODO: Use DMField instead of assuming affine */
599: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
600: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
601: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
602: for (PetscInt j = 0; j < numCIndices; ++j) {
603: PetscReal xr[8];
604: PetscInt off = 0;
606: for (PetscInt i = 0; i < Nfc; ++i) {
607: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[j] * bs[i] + b];
608: }
609: 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);
610: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, xr, &xi[j * dim]);
611: }
612: if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
613: else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
614: /* Get elemMat entries by multiplying by weight */
615: PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
616: for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
617: for (PetscInt j = 0; j < numCIndices; ++j) {
618: for (PetscInt c = 0; c < Nc; ++c) {
619: // TODO Need field offset on particle and field here
620: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
621: elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
622: }
623: }
624: }
625: for (PetscInt j = 0; j < numCIndices; ++j)
626: // TODO Need field offset on particle here
627: for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
628: if (0) PetscCall(DMPrintCellMatrix(cell, name, numCIndices * Nc, numFIndices, elemMat));
629: PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
630: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
631: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
632: PetscCall(PetscTabulationDestroy(&Tcoarse));
633: }
634: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
635: }
636: PetscCall(PetscFree3(elemMat, rowIDXs, xi));
637: PetscCall(DMSwarmSortRestoreAccess(dmc));
638: PetscCall(PetscFree3(v0, J, invJ));
639: PetscCall(PetscFree2(coordVals, bs));
640: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
641: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: /* Returns empty matrix for use with SNES FD */
646: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
647: {
648: Vec field;
649: PetscInt size;
651: PetscFunctionBegin;
652: PetscCall(DMGetGlobalVector(sw, &field));
653: PetscCall(VecGetLocalSize(field, &size));
654: PetscCall(DMRestoreGlobalVector(sw, &field));
655: PetscCall(MatCreate(PETSC_COMM_WORLD, m));
656: PetscCall(MatSetFromOptions(*m));
657: PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
658: PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
659: PetscCall(MatZeroEntries(*m));
660: PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
661: PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
662: PetscCall(MatShift(*m, 1.0));
663: PetscCall(MatSetDM(*m, sw));
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: /* FEM cols, Particle rows */
668: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
669: {
670: DMSwarmCellDM celldm;
671: PetscSection gsf;
672: PetscInt m, n, Np, bs;
673: void *ctx;
675: PetscFunctionBegin;
676: PetscCall(DMSwarmGetCellDMActive(dmCoarse, &celldm));
677: PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)dmCoarse), PETSC_ERR_USER, "Active cell DM does not define any fields");
678: PetscCall(DMGetGlobalSection(dmFine, &gsf));
679: PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
680: PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
681: PetscCall(DMSwarmCellDMGetBlockSize(celldm, dmCoarse, &bs));
682: n = Np * bs;
683: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
684: PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
685: PetscCall(MatSetType(*mass, dmCoarse->mattype));
686: PetscCall(DMGetApplicationContext(dmFine, &ctx));
688: PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
689: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
694: {
695: const char *name = "Mass Matrix Square";
696: MPI_Comm comm;
697: DMSwarmCellDM celldm;
698: PetscDS prob;
699: PetscSection fsection, globalFSection;
700: PetscHSetIJ ht;
701: PetscLayout rLayout, colLayout;
702: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
703: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
704: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
705: PetscScalar *elemMat, *elemMatSq;
706: PetscInt cdim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0;
707: const char **coordFields;
708: PetscReal **coordVals;
709: PetscInt *bs;
711: PetscFunctionBegin;
712: PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
713: PetscCall(DMGetCoordinateDim(dmf, &cdim));
714: PetscCall(DMGetDS(dmf, &prob));
715: PetscCall(PetscDSGetNumFields(prob, &Nf));
716: PetscCall(PetscDSGetTotalDimension(prob, &totDim));
717: PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
718: PetscCall(DMGetLocalSection(dmf, &fsection));
719: PetscCall(DMGetGlobalSection(dmf, &globalFSection));
720: PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
721: PetscCall(MatGetLocalSize(mass, &locRows, &locCols));
723: PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
724: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
725: PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));
727: PetscCall(PetscLayoutCreate(comm, &colLayout));
728: PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
729: PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
730: PetscCall(PetscLayoutSetUp(colLayout));
731: PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
732: PetscCall(PetscLayoutDestroy(&colLayout));
734: PetscCall(PetscLayoutCreate(comm, &rLayout));
735: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
736: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
737: PetscCall(PetscLayoutSetUp(rLayout));
738: PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
739: PetscCall(PetscLayoutDestroy(&rLayout));
741: PetscCall(DMPlexGetDepth(dmf, &depth));
742: PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
743: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
744: PetscCall(PetscMalloc1(maxAdjSize, &adj));
746: PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
747: PetscCall(PetscHSetIJCreate(&ht));
748: /* Count nonzeros
749: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
750: */
751: PetscCall(DMSwarmSortGetAccess(dmc));
752: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
753: PetscInt *cindices;
754: PetscInt numCIndices;
755: #if 0
756: PetscInt adjSize = maxAdjSize, a, j;
757: #endif
759: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
760: maxC = PetscMax(maxC, numCIndices);
761: /* Diagonal block */
762: for (PetscInt i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
763: #if 0
764: /* Off-diagonal blocks */
765: PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
766: for (a = 0; a < adjSize; ++a) {
767: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
768: const PetscInt ncell = adj[a];
769: PetscInt *ncindices;
770: PetscInt numNCIndices;
772: PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
773: {
774: PetscHashIJKey key;
775: PetscBool missing;
777: for (i = 0; i < numCIndices; ++i) {
778: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
779: if (key.i < 0) continue;
780: for (j = 0; j < numNCIndices; ++j) {
781: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
782: if (key.j < 0) continue;
783: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
784: if (missing) {
785: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
786: else ++onz[key.i - rStart];
787: }
788: }
789: }
790: }
791: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
792: }
793: }
794: #endif
795: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
796: }
797: PetscCall(PetscHSetIJDestroy(&ht));
798: PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
799: PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
800: PetscCall(PetscFree2(dnz, onz));
801: PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
802: /* Fill in values
803: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
804: Start just by producing block diagonal
805: Could loop over adjacent cells
806: Produce neighboring element matrix
807: TODO Determine which columns and rows correspond to shared dual vector
808: Do MatMatMult with rectangular matrices
809: Insert block
810: */
811: for (PetscInt field = 0; field < Nf; ++field) {
812: PetscTabulation Tcoarse;
813: PetscObject obj;
814: PetscInt Nc;
816: PetscCall(PetscDSGetDiscretization(prob, field, &obj));
817: PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
818: PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
819: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
820: for (PetscInt cell = cStart; cell < cEnd; ++cell) {
821: PetscInt *findices, *cindices;
822: PetscInt numFIndices, numCIndices;
824: /* TODO: Use DMField instead of assuming affine */
825: PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
826: PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
827: PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
828: for (PetscInt p = 0; p < numCIndices; ++p) {
829: PetscReal xr[8];
830: PetscInt off = 0;
832: for (PetscInt i = 0; i < Nfc; ++i) {
833: for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[p] * bs[i] + b];
834: }
835: 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);
836: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[p * cdim]);
837: }
838: PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
839: /* Get elemMat entries by multiplying by weight */
840: PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
841: for (PetscInt i = 0; i < numFIndices; ++i) {
842: for (PetscInt p = 0; p < numCIndices; ++p) {
843: for (PetscInt c = 0; c < Nc; ++c) {
844: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
845: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
846: }
847: }
848: }
849: PetscCall(PetscTabulationDestroy(&Tcoarse));
850: for (PetscInt p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
851: if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
852: /* Block diagonal */
853: if (numCIndices) {
854: PetscBLASInt blasn, blask;
855: PetscScalar one = 1.0, zero = 0.0;
857: PetscCall(PetscBLASIntCast(numCIndices, &blasn));
858: PetscCall(PetscBLASIntCast(numFIndices, &blask));
859: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
860: }
861: PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
862: /* TODO off-diagonal */
863: PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
864: PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
865: }
866: for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
867: }
868: PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
869: PetscCall(PetscFree(adj));
870: PetscCall(DMSwarmSortRestoreAccess(dmc));
871: PetscCall(PetscFree3(v0, J, invJ));
872: PetscCall(PetscFree2(coordVals, bs));
873: PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
874: PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
875: PetscFunctionReturn(PETSC_SUCCESS);
876: }
878: /*@
879: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
881: Collective
883: Input Parameters:
884: + dmCoarse - a `DMSWARM`
885: - dmFine - a `DMPLEX`
887: Output Parameter:
888: . mass - the square of the particle mass matrix
890: Level: advanced
892: Note:
893: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
894: future to compute the full normal equations.
896: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
897: @*/
898: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
899: {
900: PetscInt n;
901: void *ctx;
903: PetscFunctionBegin;
904: PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
905: PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
906: PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
907: PetscCall(MatSetType(*mass, dmCoarse->mattype));
908: PetscCall(DMGetApplicationContext(dmFine, &ctx));
910: PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
911: PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: /*@
916: DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
918: Collective
920: Input Parameters:
921: + dm - a `DMSWARM`
922: - fieldname - the textual name given to a registered field
924: Output Parameter:
925: . vec - the vector
927: Level: beginner
929: Note:
930: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.
932: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
933: @*/
934: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
935: {
936: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
938: PetscFunctionBegin;
940: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: /*@
945: DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
947: Collective
949: Input Parameters:
950: + dm - a `DMSWARM`
951: - fieldname - the textual name given to a registered field
953: Output Parameter:
954: . vec - the vector
956: Level: beginner
958: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
959: @*/
960: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
961: {
962: PetscFunctionBegin;
964: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: /*@
969: DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field
971: Collective
973: Input Parameters:
974: + dm - a `DMSWARM`
975: - fieldname - the textual name given to a registered field
977: Output Parameter:
978: . vec - the vector
980: Level: beginner
982: Note:
983: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
985: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
986: @*/
987: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
988: {
989: MPI_Comm comm = PETSC_COMM_SELF;
991: PetscFunctionBegin;
992: PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
993: PetscFunctionReturn(PETSC_SUCCESS);
994: }
996: /*@
997: DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field
999: Collective
1001: Input Parameters:
1002: + dm - a `DMSWARM`
1003: - fieldname - the textual name given to a registered field
1005: Output Parameter:
1006: . vec - the vector
1008: Level: beginner
1010: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
1011: @*/
1012: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1013: {
1014: PetscFunctionBegin;
1016: PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }
1020: /*@
1021: DMSwarmCreateGlobalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1023: Collective
1025: Input Parameters:
1026: + dm - a `DMSWARM`
1027: . Nf - the number of fields
1028: - fieldnames - the textual names given to the registered fields
1030: Output Parameter:
1031: . vec - the vector
1033: Level: beginner
1035: Notes:
1036: The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromFields()`.
1038: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateGlobalVectorFromField()`
1040: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromFields()`
1041: @*/
1042: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1043: {
1044: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
1046: PetscFunctionBegin;
1048: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1049: PetscFunctionReturn(PETSC_SUCCESS);
1050: }
1052: /*@
1053: DMSwarmDestroyGlobalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1055: Collective
1057: Input Parameters:
1058: + dm - a `DMSWARM`
1059: . Nf - the number of fields
1060: - fieldnames - the textual names given to the registered fields
1062: Output Parameter:
1063: . vec - the vector
1065: Level: beginner
1067: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1068: @*/
1069: PetscErrorCode DMSwarmDestroyGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1070: {
1071: PetscFunctionBegin;
1073: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: /*@
1078: DMSwarmCreateLocalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set
1080: Collective
1082: Input Parameters:
1083: + dm - a `DMSWARM`
1084: . Nf - the number of fields
1085: - fieldnames - the textual names given to the registered fields
1087: Output Parameter:
1088: . vec - the vector
1090: Level: beginner
1092: Notes:
1093: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
1095: This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateLocalVectorFromField()`
1097: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1098: @*/
1099: PetscErrorCode DMSwarmCreateLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1100: {
1101: MPI_Comm comm = PETSC_COMM_SELF;
1103: PetscFunctionBegin;
1104: PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: /*@
1109: DMSwarmDestroyLocalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set
1111: Collective
1113: Input Parameters:
1114: + dm - a `DMSWARM`
1115: . Nf - the number of fields
1116: - fieldnames - the textual names given to the registered fields
1118: Output Parameter:
1119: . vec - the vector
1121: Level: beginner
1123: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromFields()`
1124: @*/
1125: PetscErrorCode DMSwarmDestroyLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1126: {
1127: PetscFunctionBegin;
1129: PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }
1133: /*@
1134: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`
1136: Collective
1138: Input Parameter:
1139: . dm - a `DMSWARM`
1141: Level: beginner
1143: Note:
1144: After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.
1146: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1147: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1148: @*/
1149: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
1150: {
1151: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1153: PetscFunctionBegin;
1154: if (!swarm->field_registration_initialized) {
1155: swarm->field_registration_initialized = PETSC_TRUE;
1156: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
1157: PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT)); /* used for communication */
1158: }
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: /*@
1163: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`
1165: Collective
1167: Input Parameter:
1168: . dm - a `DMSWARM`
1170: Level: beginner
1172: Note:
1173: After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.
1175: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1176: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1177: @*/
1178: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
1179: {
1180: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1182: PetscFunctionBegin;
1183: if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
1184: swarm->field_registration_finalized = PETSC_TRUE;
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: /*@
1189: DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`
1191: Not Collective
1193: Input Parameters:
1194: + dm - a `DMSWARM`
1195: . nlocal - the length of each registered field
1196: - buffer - the length of the buffer used to efficient dynamic re-sizing
1198: Level: beginner
1200: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
1201: @*/
1202: PetscErrorCode DMSwarmSetLocalSizes(DM dm, PetscInt nlocal, PetscInt buffer)
1203: {
1204: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1206: PetscFunctionBegin;
1207: PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
1208: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
1209: PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: /*@
1214: DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`
1216: Collective
1218: Input Parameters:
1219: + sw - a `DMSWARM`
1220: - dm - the `DM` to attach to the `DMSWARM`
1222: Level: beginner
1224: Note:
1225: The attached `DM` (dm) will be queried for point location and
1226: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1228: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
1229: @*/
1230: PetscErrorCode DMSwarmSetCellDM(DM sw, DM dm)
1231: {
1232: DMSwarmCellDM celldm;
1233: const char *name;
1234: char *coordName;
1236: PetscFunctionBegin;
1239: PetscCall(PetscStrallocpy(DMSwarmPICField_coor, &coordName));
1240: PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 1, (const char **)&coordName, &celldm));
1241: PetscCall(PetscFree(coordName));
1242: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1243: PetscCall(DMSwarmAddCellDM(sw, celldm));
1244: PetscCall(DMSwarmCellDMDestroy(&celldm));
1245: PetscCall(DMSwarmSetCellDMActive(sw, name));
1246: PetscFunctionReturn(PETSC_SUCCESS);
1247: }
1249: /*@
1250: DMSwarmGetCellDM - Fetches the active cell `DM`
1252: Collective
1254: Input Parameter:
1255: . sw - a `DMSWARM`
1257: Output Parameter:
1258: . dm - the active `DM` for the `DMSWARM`
1260: Level: beginner
1262: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1263: @*/
1264: PetscErrorCode DMSwarmGetCellDM(DM sw, DM *dm)
1265: {
1266: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1267: DMSwarmCellDM celldm;
1269: PetscFunctionBegin;
1270: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)&celldm));
1271: PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "There is no cell DM named %s in this Swarm", swarm->activeCellDM);
1272: PetscCall(DMSwarmCellDMGetDM(celldm, dm));
1273: PetscFunctionReturn(PETSC_SUCCESS);
1274: }
1276: /*@
1277: DMSwarmSetCellDMActive - Activates a cell `DM` for a `DMSWARM`
1279: Collective
1281: Input Parameters:
1282: + sw - a `DMSWARM`
1283: - name - name of the cell `DM` to active for the `DMSWARM`
1285: Level: beginner
1287: Note:
1288: The attached `DM` (dmcell) will be queried for point location and
1289: neighbor MPI-rank information if `DMSwarmMigrate()` is called.
1291: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1292: @*/
1293: PetscErrorCode DMSwarmSetCellDMActive(DM sw, const char name[])
1294: {
1295: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1296: DMSwarmCellDM celldm;
1298: PetscFunctionBegin;
1300: PetscCall(PetscInfo(sw, "Setting cell DM to %s\n", name));
1301: PetscCall(PetscFree(swarm->activeCellDM));
1302: PetscCall(PetscStrallocpy(name, (char **)&swarm->activeCellDM));
1303: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1304: PetscFunctionReturn(PETSC_SUCCESS);
1305: }
1307: /*@
1308: DMSwarmGetCellDMActive - Returns the active cell `DM` for a `DMSWARM`
1310: Collective
1312: Input Parameter:
1313: . sw - a `DMSWARM`
1315: Output Parameter:
1316: . celldm - the active `DMSwarmCellDM`
1318: Level: beginner
1320: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1321: @*/
1322: PetscErrorCode DMSwarmGetCellDMActive(DM sw, DMSwarmCellDM *celldm)
1323: {
1324: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1326: PetscFunctionBegin;
1328: PetscCheck(swarm->activeCellDM, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no active cell DM");
1329: PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)celldm));
1330: PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has invalid cell DM for %s", swarm->activeCellDM);
1331: PetscFunctionReturn(PETSC_SUCCESS);
1332: }
1334: /*@
1335: DMSwarmAddCellDM - Adds a cell `DM` to the `DMSWARM`
1337: Collective
1339: Input Parameters:
1340: + sw - a `DMSWARM`
1341: - celldm - the `DMSwarmCellDM`
1343: Level: beginner
1345: Note:
1346: Cell DMs with the same name will share the cellid field
1348: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPushCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1349: @*/
1350: PetscErrorCode DMSwarmAddCellDM(DM sw, DMSwarmCellDM celldm)
1351: {
1352: DM_Swarm *swarm = (DM_Swarm *)sw->data;
1353: const char *name;
1354: PetscInt dim;
1355: PetscBool flg;
1356: MPI_Comm comm;
1358: PetscFunctionBegin;
1360: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1362: PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1363: PetscCall(PetscObjectListAdd(&swarm->cellDMs, name, (PetscObject)celldm));
1364: PetscCall(DMGetDimension(sw, &dim));
1365: for (PetscInt f = 0; f < celldm->Nfc; ++f) {
1366: PetscCall(DMSwarmDataFieldStringInList(celldm->coordFields[f], swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1367: if (!flg) {
1368: PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->coordFields[f], dim, PETSC_DOUBLE));
1369: } else {
1370: PetscDataType dt;
1371: PetscInt bs;
1373: PetscCall(DMSwarmGetFieldInfo(sw, celldm->coordFields[f], &bs, &dt));
1374: PetscCheck(bs == dim, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has blocksize %" PetscInt_FMT " != %" PetscInt_FMT " spatial dimension", celldm->coordFields[f], bs, dim);
1375: PetscCheck(dt == PETSC_DOUBLE, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has datatype %s != PETSC_DOUBLE", celldm->coordFields[f], PetscDataTypes[dt]);
1376: }
1377: }
1378: // Assume that DMs with the same name share the cellid field
1379: PetscCall(DMSwarmDataFieldStringInList(celldm->cellid, swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1380: if (!flg) {
1381: PetscBool isShell, isDummy;
1382: const char *name;
1384: // Allow dummy DMSHELL (I don't think we should support this mode)
1385: PetscCall(PetscObjectTypeCompare((PetscObject)celldm->dm, DMSHELL, &isShell));
1386: PetscCall(PetscObjectGetName((PetscObject)celldm->dm, &name));
1387: PetscCall(PetscStrcmp(name, "dummy", &isDummy));
1388: if (!isShell || !isDummy) PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->cellid, 1, PETSC_INT));
1389: }
1390: PetscCall(DMSwarmSetCellDMActive(sw, name));
1391: PetscFunctionReturn(PETSC_SUCCESS);
1392: }
1394: /*@
1395: DMSwarmGetLocalSize - Retrieves the local length of fields registered
1397: Not Collective
1399: Input Parameter:
1400: . dm - a `DMSWARM`
1402: Output Parameter:
1403: . nlocal - the length of each registered field
1405: Level: beginner
1407: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
1408: @*/
1409: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
1410: {
1411: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1413: PetscFunctionBegin;
1414: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
1415: PetscFunctionReturn(PETSC_SUCCESS);
1416: }
1418: /*@
1419: DMSwarmGetSize - Retrieves the total length of fields registered
1421: Collective
1423: Input Parameter:
1424: . dm - a `DMSWARM`
1426: Output Parameter:
1427: . n - the total length of each registered field
1429: Level: beginner
1431: Note:
1432: This calls `MPI_Allreduce()` upon each call (inefficient but safe)
1434: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
1435: @*/
1436: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
1437: {
1438: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1439: PetscInt nlocal;
1441: PetscFunctionBegin;
1442: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1443: PetscCallMPI(MPIU_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1444: PetscFunctionReturn(PETSC_SUCCESS);
1445: }
1447: /*@
1448: DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type
1450: Collective
1452: Input Parameters:
1453: + dm - a `DMSWARM`
1454: . fieldname - the textual name to identify this field
1455: . blocksize - the number of each data type
1456: - type - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)
1458: Level: beginner
1460: Notes:
1461: The textual name for each registered field must be unique.
1463: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1464: @*/
1465: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1466: {
1467: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1468: size_t size;
1470: PetscFunctionBegin;
1471: PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1472: PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1474: PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1475: PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1476: PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1477: PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1478: PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1480: PetscCall(PetscDataTypeGetSize(type, &size));
1481: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1482: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1483: {
1484: DMSwarmDataField gfield;
1486: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1487: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1488: }
1489: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1490: PetscFunctionReturn(PETSC_SUCCESS);
1491: }
1493: /*@C
1494: DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`
1496: Collective
1498: Input Parameters:
1499: + dm - a `DMSWARM`
1500: . fieldname - the textual name to identify this field
1501: - size - the size in bytes of the user struct of each data type
1503: Level: beginner
1505: Note:
1506: The textual name for each registered field must be unique.
1508: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1509: @*/
1510: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1511: {
1512: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1514: PetscFunctionBegin;
1515: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1516: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1517: PetscFunctionReturn(PETSC_SUCCESS);
1518: }
1520: /*@C
1521: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`
1523: Collective
1525: Input Parameters:
1526: + dm - a `DMSWARM`
1527: . fieldname - the textual name to identify this field
1528: . size - the size in bytes of the user data type
1529: - blocksize - the number of each data type
1531: Level: beginner
1533: Note:
1534: The textual name for each registered field must be unique.
1536: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1537: @*/
1538: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1539: {
1540: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1542: PetscFunctionBegin;
1543: PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1544: {
1545: DMSwarmDataField gfield;
1547: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1548: PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1549: }
1550: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1551: PetscFunctionReturn(PETSC_SUCCESS);
1552: }
1554: /*@C
1555: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1557: Not Collective, No Fortran Support
1559: Input Parameters:
1560: + dm - a `DMSWARM`
1561: - fieldname - the textual name to identify this field
1563: Output Parameters:
1564: + blocksize - the number of each data type
1565: . type - the data type
1566: - data - pointer to raw array
1568: Level: beginner
1570: Notes:
1571: The array must be returned using a matching call to `DMSwarmRestoreField()`.
1573: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1574: @*/
1575: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1576: {
1577: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1578: DMSwarmDataField gfield;
1580: PetscFunctionBegin;
1582: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1583: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1584: PetscCall(DMSwarmDataFieldGetAccess(gfield));
1585: PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1586: if (blocksize) *blocksize = gfield->bs;
1587: if (type) *type = gfield->petsc_type;
1588: PetscFunctionReturn(PETSC_SUCCESS);
1589: }
1591: /*@C
1592: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1594: Not Collective, No Fortran Support
1596: Input Parameters:
1597: + dm - a `DMSWARM`
1598: - fieldname - the textual name to identify this field
1600: Output Parameters:
1601: + blocksize - the number of each data type
1602: . type - the data type
1603: - data - pointer to raw array
1605: Level: beginner
1607: Notes:
1608: The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.
1610: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1611: @*/
1612: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1613: {
1614: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1615: DMSwarmDataField gfield;
1617: PetscFunctionBegin;
1619: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1620: PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1621: if (data) *data = NULL;
1622: PetscFunctionReturn(PETSC_SUCCESS);
1623: }
1625: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
1626: {
1627: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1628: DMSwarmDataField gfield;
1630: PetscFunctionBegin;
1632: PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1633: if (blocksize) *blocksize = gfield->bs;
1634: if (type) *type = gfield->petsc_type;
1635: PetscFunctionReturn(PETSC_SUCCESS);
1636: }
1638: /*@
1639: DMSwarmAddPoint - Add space for one new point in the `DMSWARM`
1641: Not Collective
1643: Input Parameter:
1644: . dm - a `DMSWARM`
1646: Level: beginner
1648: Notes:
1649: The new point will have all fields initialized to zero.
1651: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
1652: @*/
1653: PetscErrorCode DMSwarmAddPoint(DM dm)
1654: {
1655: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1657: PetscFunctionBegin;
1658: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1659: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1660: PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
1661: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: /*@
1666: DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`
1668: Not Collective
1670: Input Parameters:
1671: + dm - a `DMSWARM`
1672: - npoints - the number of new points to add
1674: Level: beginner
1676: Notes:
1677: The new point will have all fields initialized to zero.
1679: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
1680: @*/
1681: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1682: {
1683: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1684: PetscInt nlocal;
1686: PetscFunctionBegin;
1687: PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1688: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1689: nlocal = nlocal + npoints;
1690: PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
1691: PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1692: PetscFunctionReturn(PETSC_SUCCESS);
1693: }
1695: /*@
1696: DMSwarmRemovePoint - Remove the last point from the `DMSWARM`
1698: Not Collective
1700: Input Parameter:
1701: . dm - a `DMSWARM`
1703: Level: beginner
1705: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
1706: @*/
1707: PetscErrorCode DMSwarmRemovePoint(DM dm)
1708: {
1709: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1711: PetscFunctionBegin;
1712: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1713: PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
1714: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1715: PetscFunctionReturn(PETSC_SUCCESS);
1716: }
1718: /*@
1719: DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`
1721: Not Collective
1723: Input Parameters:
1724: + dm - a `DMSWARM`
1725: - idx - index of point to remove
1727: Level: beginner
1729: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1730: @*/
1731: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
1732: {
1733: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1735: PetscFunctionBegin;
1736: PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1737: PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
1738: PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: }
1742: /*@
1743: DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`
1745: Not Collective
1747: Input Parameters:
1748: + dm - a `DMSWARM`
1749: . pi - the index of the point to copy
1750: - pj - the point index where the copy should be located
1752: Level: beginner
1754: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
1755: @*/
1756: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
1757: {
1758: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1760: PetscFunctionBegin;
1761: if (!swarm->issetup) PetscCall(DMSetUp(dm));
1762: PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
1763: PetscFunctionReturn(PETSC_SUCCESS);
1764: }
1766: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
1767: {
1768: PetscFunctionBegin;
1769: PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
1770: PetscFunctionReturn(PETSC_SUCCESS);
1771: }
1773: /*@
1774: DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks
1776: Collective
1778: Input Parameters:
1779: + dm - the `DMSWARM`
1780: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1782: Level: advanced
1784: Notes:
1785: The `DM` will be modified to accommodate received points.
1786: If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
1787: Different styles of migration are supported. See `DMSwarmSetMigrateType()`.
1789: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
1790: @*/
1791: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
1792: {
1793: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1795: PetscFunctionBegin;
1796: PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
1797: switch (swarm->migrate_type) {
1798: case DMSWARM_MIGRATE_BASIC:
1799: PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
1800: break;
1801: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1802: PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
1803: break;
1804: case DMSWARM_MIGRATE_DMCELLEXACT:
1805: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1806: case DMSWARM_MIGRATE_USER:
1807: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
1808: default:
1809: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
1810: }
1811: PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
1812: PetscCall(DMClearGlobalVectors(dm));
1813: PetscFunctionReturn(PETSC_SUCCESS);
1814: }
1816: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
1818: /*
1819: DMSwarmCollectViewCreate
1821: * Applies a collection method and gathers point neighbour points into dm
1823: Notes:
1824: Users should call DMSwarmCollectViewDestroy() after
1825: they have finished computations associated with the collected points
1826: */
1828: /*@
1829: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1830: in neighbour ranks into the `DMSWARM`
1832: Collective
1834: Input Parameter:
1835: . dm - the `DMSWARM`
1837: Level: advanced
1839: Notes:
1840: Users should call `DMSwarmCollectViewDestroy()` after
1841: they have finished computations associated with the collected points
1843: Different collect methods are supported. See `DMSwarmSetCollectType()`.
1845: Developer Note:
1846: Create and Destroy routines create new objects that can get destroyed, they do not change the state
1847: of the current object.
1849: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
1850: @*/
1851: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1852: {
1853: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1854: PetscInt ng;
1856: PetscFunctionBegin;
1857: PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
1858: PetscCall(DMSwarmGetLocalSize(dm, &ng));
1859: switch (swarm->collect_type) {
1860: case DMSWARM_COLLECT_BASIC:
1861: PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
1862: break;
1863: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1864: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1865: case DMSWARM_COLLECT_GENERAL:
1866: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
1867: default:
1868: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
1869: }
1870: swarm->collect_view_active = PETSC_TRUE;
1871: swarm->collect_view_reset_nlocal = ng;
1872: PetscFunctionReturn(PETSC_SUCCESS);
1873: }
1875: /*@
1876: DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`
1878: Collective
1880: Input Parameters:
1881: . dm - the `DMSWARM`
1883: Notes:
1884: Users should call `DMSwarmCollectViewCreate()` before this function is called.
1886: Level: advanced
1888: Developer Note:
1889: Create and Destroy routines create new objects that can get destroyed, they do not change the state
1890: of the current object.
1892: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
1893: @*/
1894: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1895: {
1896: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1898: PetscFunctionBegin;
1899: PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
1900: PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
1901: swarm->collect_view_active = PETSC_FALSE;
1902: PetscFunctionReturn(PETSC_SUCCESS);
1903: }
1905: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
1906: {
1907: PetscInt dim;
1909: PetscFunctionBegin;
1910: PetscCall(DMSwarmSetNumSpecies(dm, 1));
1911: PetscCall(DMGetDimension(dm, &dim));
1912: PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1913: PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
1914: PetscFunctionReturn(PETSC_SUCCESS);
1915: }
1917: /*@
1918: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
1920: Collective
1922: Input Parameters:
1923: + dm - the `DMSWARM`
1924: - Npc - The number of particles per cell in the cell `DM`
1926: Level: intermediate
1928: Notes:
1929: The user must use `DMSwarmSetCellDM()` to set the cell `DM` first. The particles are placed randomly inside each cell. If only
1930: one particle is in each cell, it is placed at the centroid.
1932: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1933: @*/
1934: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1935: {
1936: DM cdm;
1937: DMSwarmCellDM celldm;
1938: PetscRandom rnd;
1939: DMPolytopeType ct;
1940: PetscBool simplex;
1941: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1942: PetscInt dim, d, cStart, cEnd, c, p, Nfc;
1943: const char **coordFields;
1945: PetscFunctionBeginUser;
1946: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
1947: PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
1948: PetscCall(PetscRandomSetType(rnd, PETSCRAND48));
1950: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
1951: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1952: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
1953: PetscCall(DMSwarmGetCellDM(dm, &cdm));
1954: PetscCall(DMGetDimension(cdm, &dim));
1955: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
1956: PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
1957: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1959: PetscCall(PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
1960: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1961: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&coords));
1962: for (c = cStart; c < cEnd; ++c) {
1963: if (Npc == 1) {
1964: PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
1965: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
1966: } else {
1967: PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
1968: for (p = 0; p < Npc; ++p) {
1969: const PetscInt n = c * Npc + p;
1970: PetscReal sum = 0.0, refcoords[3];
1972: for (d = 0; d < dim; ++d) {
1973: PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
1974: sum += refcoords[d];
1975: }
1976: if (simplex && sum > 0.0)
1977: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
1978: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
1979: }
1980: }
1981: }
1982: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&coords));
1983: PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
1984: PetscCall(PetscRandomDestroy(&rnd));
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: /*@
1989: DMSwarmSetType - Set particular flavor of `DMSWARM`
1991: Collective
1993: Input Parameters:
1994: + dm - the `DMSWARM`
1995: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)
1997: Level: advanced
1999: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2000: @*/
2001: PetscErrorCode DMSwarmSetType(DM dm, DMSwarmType stype)
2002: {
2003: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2005: PetscFunctionBegin;
2006: swarm->swarm_type = stype;
2007: if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(dm));
2008: PetscFunctionReturn(PETSC_SUCCESS);
2009: }
2011: static PetscErrorCode DMSetup_Swarm(DM dm)
2012: {
2013: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2014: PetscMPIInt rank;
2015: PetscInt p, npoints, *rankval;
2017: PetscFunctionBegin;
2018: if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
2019: swarm->issetup = PETSC_TRUE;
2021: if (swarm->swarm_type == DMSWARM_PIC) {
2022: DMSwarmCellDM celldm;
2024: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2025: PetscCheck(celldm, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "No active cell DM. DMSWARM_PIC requires you call DMSwarmSetCellDM() or DMSwarmAddCellDM()");
2026: if (celldm->dm->ops->locatepointssubdomain) {
2027: /* check methods exists for exact ownership identificiation */
2028: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
2029: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
2030: } else {
2031: /* check methods exist for point location AND rank neighbor identification */
2032: if (celldm->dm->ops->locatepoints) {
2033: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));
2034: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
2036: if (celldm->dm->ops->getneighbors) {
2037: PetscCall(PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));
2038: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
2040: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
2041: }
2042: }
2044: PetscCall(DMSwarmFinalizeFieldRegister(dm));
2046: /* check some fields were registered */
2047: PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
2049: /* check local sizes were set */
2050: PetscCheck(swarm->db->L != -1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Local sizes must be set via DMSwarmSetLocalSizes()");
2052: /* initialize values in pid and rank placeholders */
2053: /* TODO: [pid - use MPI_Scan] */
2054: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
2055: PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL));
2056: PetscCall(DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
2057: for (p = 0; p < npoints; p++) rankval[p] = rank;
2058: PetscCall(DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
2059: PetscFunctionReturn(PETSC_SUCCESS);
2060: }
2062: static PetscErrorCode DMDestroy_Swarm(DM dm)
2063: {
2064: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2066: PetscFunctionBegin;
2067: if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2068: PetscCall(PetscObjectListDestroy(&swarm->cellDMs));
2069: PetscCall(PetscFree(swarm->activeCellDM));
2070: PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
2071: PetscCall(PetscFree(swarm));
2072: PetscFunctionReturn(PETSC_SUCCESS);
2073: }
2075: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
2076: {
2077: DM cdm;
2078: DMSwarmCellDM celldm;
2079: PetscDraw draw;
2080: PetscReal *coords, oldPause, radius = 0.01;
2081: PetscInt Np, p, bs, Nfc;
2082: const char **coordFields;
2084: PetscFunctionBegin;
2085: PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
2086: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2087: PetscCall(DMSwarmGetCellDM(dm, &cdm));
2088: PetscCall(PetscDrawGetPause(draw, &oldPause));
2089: PetscCall(PetscDrawSetPause(draw, 0.0));
2090: PetscCall(DMView(cdm, viewer));
2091: PetscCall(PetscDrawSetPause(draw, oldPause));
2093: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2094: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2095: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2096: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2097: PetscCall(DMSwarmGetField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2098: for (p = 0; p < Np; ++p) {
2099: const PetscInt i = p * bs;
2101: PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
2102: }
2103: PetscCall(DMSwarmRestoreField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2104: PetscCall(PetscDrawFlush(draw));
2105: PetscCall(PetscDrawPause(draw));
2106: PetscCall(PetscDrawSave(draw));
2107: PetscFunctionReturn(PETSC_SUCCESS);
2108: }
2110: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
2111: {
2112: PetscViewerFormat format;
2113: PetscInt *sizes;
2114: PetscInt dim, Np, maxSize = 17;
2115: MPI_Comm comm;
2116: PetscMPIInt rank, size, p;
2117: const char *name, *cellid;
2119: PetscFunctionBegin;
2120: PetscCall(PetscViewerGetFormat(viewer, &format));
2121: PetscCall(DMGetDimension(dm, &dim));
2122: PetscCall(DMSwarmGetLocalSize(dm, &Np));
2123: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2124: PetscCallMPI(MPI_Comm_rank(comm, &rank));
2125: PetscCallMPI(MPI_Comm_size(comm, &size));
2126: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
2127: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
2128: else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
2129: if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
2130: else PetscCall(PetscCalloc1(3, &sizes));
2131: if (size < maxSize) {
2132: PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
2133: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of particles per rank:"));
2134: for (p = 0; p < size; ++p) {
2135: if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
2136: }
2137: } else {
2138: PetscInt locMinMax[2] = {Np, Np};
2140: PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
2141: PetscCall(PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
2142: }
2143: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2144: PetscCall(PetscFree(sizes));
2145: if (format == PETSC_VIEWER_ASCII_INFO) {
2146: DMSwarmCellDM celldm;
2147: PetscInt *cell;
2149: PetscCall(PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n"));
2150: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2151: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2152: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
2153: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&cell));
2154: for (p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " p%d: %" PetscInt_FMT "\n", p, cell[p]));
2155: PetscCall(PetscViewerFlush(viewer));
2156: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2157: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&cell));
2158: }
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
2163: {
2164: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2165: PetscBool iascii, ibinary, isvtk, isdraw, ispython;
2166: #if defined(PETSC_HAVE_HDF5)
2167: PetscBool ishdf5;
2168: #endif
2170: PetscFunctionBegin;
2173: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2174: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
2175: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
2176: #if defined(PETSC_HAVE_HDF5)
2177: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2178: #endif
2179: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2180: PetscCall(PetscObjectHasFunction((PetscObject)viewer, "PetscViewerPythonViewObject_C", &ispython));
2181: if (iascii) {
2182: PetscViewerFormat format;
2184: PetscCall(PetscViewerGetFormat(viewer, &format));
2185: switch (format) {
2186: case PETSC_VIEWER_ASCII_INFO_DETAIL:
2187: PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
2188: break;
2189: default:
2190: PetscCall(DMView_Swarm_Ascii(dm, viewer));
2191: }
2192: } else {
2193: #if defined(PETSC_HAVE_HDF5)
2194: if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
2195: #endif
2196: if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
2197: if (ispython) PetscCall(PetscViewerPythonViewObject(viewer, (PetscObject)dm));
2198: }
2199: PetscFunctionReturn(PETSC_SUCCESS);
2200: }
2202: /*@
2203: DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
2204: 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.
2206: Noncollective
2208: Input Parameters:
2209: + sw - the `DMSWARM`
2210: . cellID - the integer id of the cell to be extracted and filtered
2211: - cellswarm - The `DMSWARM` to receive the cell
2213: Level: beginner
2215: Notes:
2216: This presently only supports `DMSWARM_PIC` type
2218: Should be restored with `DMSwarmRestoreCellSwarm()`
2220: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
2222: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
2223: @*/
2224: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2225: {
2226: DM_Swarm *original = (DM_Swarm *)sw->data;
2227: DMLabel label;
2228: DM dmc, subdmc;
2229: PetscInt *pids, particles, dim;
2230: const char *name;
2232: PetscFunctionBegin;
2233: /* Configure new swarm */
2234: PetscCall(DMSetType(cellswarm, DMSWARM));
2235: PetscCall(DMGetDimension(sw, &dim));
2236: PetscCall(DMSetDimension(cellswarm, dim));
2237: PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
2238: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
2239: PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
2240: PetscCall(DMSwarmSortGetAccess(sw));
2241: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
2242: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2243: PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
2244: PetscCall(DMSwarmSortRestoreAccess(sw));
2245: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2246: PetscCall(DMSwarmGetCellDM(sw, &dmc));
2247: PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
2248: PetscCall(DMAddLabel(dmc, label));
2249: PetscCall(DMLabelSetValue(label, cellID, 1));
2250: PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdmc));
2251: PetscCall(PetscObjectGetName((PetscObject)dmc, &name));
2252: PetscCall(PetscObjectSetName((PetscObject)subdmc, name));
2253: PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
2254: PetscCall(DMLabelDestroy(&label));
2255: PetscFunctionReturn(PETSC_SUCCESS);
2256: }
2258: /*@
2259: DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.
2261: Noncollective
2263: Input Parameters:
2264: + sw - the parent `DMSWARM`
2265: . cellID - the integer id of the cell to be copied back into the parent swarm
2266: - cellswarm - the cell swarm object
2268: Level: beginner
2270: Note:
2271: This only supports `DMSWARM_PIC` types of `DMSWARM`s
2273: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
2274: @*/
2275: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2276: {
2277: DM dmc;
2278: PetscInt *pids, particles, p;
2280: PetscFunctionBegin;
2281: PetscCall(DMSwarmSortGetAccess(sw));
2282: PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2283: PetscCall(DMSwarmSortRestoreAccess(sw));
2284: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
2285: for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
2286: /* Free memory, destroy cell dm */
2287: PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
2288: PetscCall(DMDestroy(&dmc));
2289: PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2290: PetscFunctionReturn(PETSC_SUCCESS);
2291: }
2293: /*@
2294: DMSwarmComputeMoments - Compute the first three particle moments for a given field
2296: Noncollective
2298: Input Parameters:
2299: + sw - the `DMSWARM`
2300: . coordinate - the coordinate field name
2301: - weight - the weight field name
2303: Output Parameter:
2304: . moments - the field moments
2306: Level: intermediate
2308: Notes:
2309: The `moments` array should be of length bs + 2, where bs is the block size of the coordinate field.
2311: The weight field must be a scalar, having blocksize 1.
2313: .seealso: `DM`, `DMSWARM`, `DMPlexComputeMoments()`
2314: @*/
2315: PetscErrorCode DMSwarmComputeMoments(DM sw, const char coordinate[], const char weight[], PetscReal moments[])
2316: {
2317: const PetscReal *coords;
2318: const PetscReal *w;
2319: PetscReal *mom;
2320: PetscDataType dtc, dtw;
2321: PetscInt bsc, bsw, Np;
2322: MPI_Comm comm;
2324: PetscFunctionBegin;
2326: PetscAssertPointer(coordinate, 2);
2327: PetscAssertPointer(weight, 3);
2328: PetscAssertPointer(moments, 4);
2329: PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
2330: PetscCall(DMSwarmGetField(sw, coordinate, &bsc, &dtc, (void **)&coords));
2331: PetscCall(DMSwarmGetField(sw, weight, &bsw, &dtw, (void **)&w));
2332: PetscCheck(dtc == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s must be real, not %s", coordinate, PetscDataTypes[dtc]);
2333: PetscCheck(dtw == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be real, not %s", weight, PetscDataTypes[dtw]);
2334: PetscCheck(bsw == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be a scalar, not blocksize %" PetscInt_FMT, weight, bsw);
2335: PetscCall(DMSwarmGetLocalSize(sw, &Np));
2336: PetscCall(DMGetWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2337: PetscCall(PetscArrayzero(mom, bsc + 2));
2338: for (PetscInt p = 0; p < Np; ++p) {
2339: const PetscReal *c = &coords[p * bsc];
2340: const PetscReal wp = w[p];
2342: mom[0] += wp;
2343: for (PetscInt d = 0; d < bsc; ++d) {
2344: mom[d + 1] += wp * c[d];
2345: mom[d + bsc + 1] += wp * PetscSqr(c[d]);
2346: }
2347: }
2348: PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&coords));
2349: PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
2350: PetscCallMPI(MPIU_Allreduce(mom, moments, bsc + 2, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
2351: PetscCall(DMRestoreWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2352: PetscFunctionReturn(PETSC_SUCCESS);
2353: }
2355: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
2357: static PetscErrorCode DMInitialize_Swarm(DM sw)
2358: {
2359: PetscFunctionBegin;
2360: sw->ops->view = DMView_Swarm;
2361: sw->ops->load = NULL;
2362: sw->ops->setfromoptions = NULL;
2363: sw->ops->clone = DMClone_Swarm;
2364: sw->ops->setup = DMSetup_Swarm;
2365: sw->ops->createlocalsection = NULL;
2366: sw->ops->createsectionpermutation = NULL;
2367: sw->ops->createdefaultconstraints = NULL;
2368: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
2369: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
2370: sw->ops->getlocaltoglobalmapping = NULL;
2371: sw->ops->createfieldis = NULL;
2372: sw->ops->createcoordinatedm = NULL;
2373: sw->ops->getcoloring = NULL;
2374: sw->ops->creatematrix = DMCreateMatrix_Swarm;
2375: sw->ops->createinterpolation = NULL;
2376: sw->ops->createinjection = NULL;
2377: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
2378: sw->ops->refine = NULL;
2379: sw->ops->coarsen = NULL;
2380: sw->ops->refinehierarchy = NULL;
2381: sw->ops->coarsenhierarchy = NULL;
2382: sw->ops->globaltolocalbegin = NULL;
2383: sw->ops->globaltolocalend = NULL;
2384: sw->ops->localtoglobalbegin = NULL;
2385: sw->ops->localtoglobalend = NULL;
2386: sw->ops->destroy = DMDestroy_Swarm;
2387: sw->ops->createsubdm = NULL;
2388: sw->ops->getdimpoints = NULL;
2389: sw->ops->locatepoints = NULL;
2390: PetscFunctionReturn(PETSC_SUCCESS);
2391: }
2393: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
2394: {
2395: DM_Swarm *swarm = (DM_Swarm *)dm->data;
2397: PetscFunctionBegin;
2398: swarm->refct++;
2399: (*newdm)->data = swarm;
2400: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
2401: PetscCall(DMInitialize_Swarm(*newdm));
2402: (*newdm)->dim = dm->dim;
2403: PetscFunctionReturn(PETSC_SUCCESS);
2404: }
2406: /*MC
2407: DMSWARM = "swarm" - A `DM` object used to represent arrays of data (fields) of arbitrary data type.
2408: This implementation was designed for particle methods in which the underlying
2409: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
2411: Level: intermediate
2413: Notes:
2414: User data can be represented by `DMSWARM` through a registering "fields".
2415: To register a field, the user must provide;
2416: (a) a unique name;
2417: (b) the data type (or size in bytes);
2418: (c) the block size of the data.
2420: For example, suppose the application requires a unique id, energy, momentum and density to be stored
2421: on a set of particles. Then the following code could be used
2422: .vb
2423: DMSwarmInitializeFieldRegister(dm)
2424: DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
2425: DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
2426: DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
2427: DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
2428: DMSwarmFinalizeFieldRegister(dm)
2429: .ve
2431: The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
2432: The only restriction imposed by `DMSWARM` is that all fields contain the same number of points.
2434: To support particle methods, "migration" techniques are provided. These methods migrate data
2435: between ranks.
2437: `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
2438: As a `DMSWARM` may internally define and store values of different data types,
2439: before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
2440: fields should be used to define a `Vec` object via
2441: `DMSwarmVectorDefineField()`
2442: The specified field can be changed at any time - thereby permitting vectors
2443: compatible with different fields to be created.
2445: A dual representation of fields in the `DMSWARM` and a Vec object is permitted via
2446: `DMSwarmCreateGlobalVectorFromField()`
2447: Here the data defining the field in the `DMSWARM` is shared with a Vec.
2448: This is inherently unsafe if you alter the size of the field at any time between
2449: calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
2450: If the local size of the `DMSWARM` does not match the local size of the global vector
2451: when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.
2453: Additional high-level support is provided for Particle-In-Cell methods.
2454: Please refer to `DMSwarmSetType()`.
2456: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`
2457: M*/
2459: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
2460: {
2461: DM_Swarm *swarm;
2463: PetscFunctionBegin;
2465: PetscCall(PetscNew(&swarm));
2466: dm->data = swarm;
2467: PetscCall(DMSwarmDataBucketCreate(&swarm->db));
2468: PetscCall(DMSwarmInitializeFieldRegister(dm));
2469: dm->dim = 0;
2470: swarm->refct = 1;
2471: swarm->issetup = PETSC_FALSE;
2472: swarm->swarm_type = DMSWARM_BASIC;
2473: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
2474: swarm->collect_type = DMSWARM_COLLECT_BASIC;
2475: swarm->migrate_error_on_missing_point = PETSC_FALSE;
2476: swarm->collect_view_active = PETSC_FALSE;
2477: swarm->collect_view_reset_nlocal = -1;
2478: PetscCall(DMInitialize_Swarm(dm));
2479: if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
2480: PetscFunctionReturn(PETSC_SUCCESS);
2481: }
2483: /* Replace dm with the contents of ndm, and then destroy ndm
2484: - Share the DM_Swarm structure
2485: */
2486: PetscErrorCode DMSwarmReplace(DM dm, DM *ndm)
2487: {
2488: DM dmNew = *ndm;
2489: const PetscReal *maxCell, *Lstart, *L;
2490: PetscInt dim;
2492: PetscFunctionBegin;
2493: if (dm == dmNew) {
2494: PetscCall(DMDestroy(ndm));
2495: PetscFunctionReturn(PETSC_SUCCESS);
2496: }
2497: dm->setupcalled = dmNew->setupcalled;
2498: if (!dm->hdr.name) {
2499: const char *name;
2501: PetscCall(PetscObjectGetName((PetscObject)*ndm, &name));
2502: PetscCall(PetscObjectSetName((PetscObject)dm, name));
2503: }
2504: PetscCall(DMGetDimension(dmNew, &dim));
2505: PetscCall(DMSetDimension(dm, dim));
2506: PetscCall(DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L));
2507: PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
2508: PetscCall(DMDestroy_Swarm(dm));
2509: PetscCall(DMInitialize_Swarm(dm));
2510: dm->data = dmNew->data;
2511: ((DM_Swarm *)dmNew->data)->refct++;
2512: PetscCall(DMDestroy(ndm));
2513: PetscFunctionReturn(PETSC_SUCCESS);
2514: }