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