Actual source code: swarm.c
1: #define PETSCDM_DLL
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 *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};
22: const char DMSwarmField_pid[] = "DMSwarm_pid";
23: const char DMSwarmField_rank[] = "DMSwarm_rank";
24: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
25: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
27: PetscInt SwarmDataFieldId = -1;
29: #if defined(PETSC_HAVE_HDF5)
30: #include <petscviewerhdf5.h>
32: PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
33: {
34: DM dm;
35: PetscReal seqval;
36: PetscInt seqnum, bs;
37: PetscBool isseq;
39: VecGetDM(v, &dm);
40: VecGetBlockSize(v, &bs);
41: PetscViewerHDF5PushGroup(viewer, "/particle_fields");
42: PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq);
43: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
44: PetscViewerHDF5SetTimestep(viewer, seqnum);
45: /* DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer); */
46: VecViewNative(v, viewer);
47: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs);
48: PetscViewerHDF5PopGroup(viewer);
49: return 0;
50: }
52: PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
53: {
54: Vec coordinates;
55: PetscInt Np;
56: PetscBool isseq;
58: DMSwarmGetSize(dm, &Np);
59: DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
60: PetscObjectSetName((PetscObject)coordinates, "coordinates");
61: PetscViewerHDF5PushGroup(viewer, "/particles");
62: PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq);
63: VecViewNative(coordinates, viewer);
64: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np);
65: PetscViewerHDF5PopGroup(viewer);
66: DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
67: return 0;
68: }
69: #endif
71: PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
72: {
73: DM dm;
74: #if defined(PETSC_HAVE_HDF5)
75: PetscBool ishdf5;
76: #endif
78: VecGetDM(v, &dm);
80: #if defined(PETSC_HAVE_HDF5)
81: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
82: if (ishdf5) {
83: VecView_Swarm_HDF5_Internal(v, viewer);
84: return 0;
85: }
86: #endif
87: VecViewNative(v, viewer);
88: return 0;
89: }
91: /*@C
92: DMSwarmVectorDefineField - Sets the field from which to define a Vec object
93: when DMCreateLocalVector(), or DMCreateGlobalVector() is called
95: Collective on dm
97: Input parameters:
98: + dm - a DMSwarm
99: - fieldname - the textual name given to a registered field
101: Level: beginner
103: Notes:
105: The field with name fieldname must be defined as having a data type of PetscScalar.
107: This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
108: Mutiple calls to DMSwarmVectorDefineField() are permitted.
110: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
111: @*/
112: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
113: {
114: DM_Swarm *swarm = (DM_Swarm *)dm->data;
115: PetscInt bs, n;
116: PetscScalar *array;
117: PetscDataType type;
119: if (!swarm->issetup) DMSetUp(dm);
120: DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
121: DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array);
123: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
125: PetscSNPrintf(swarm->vec_field_name, PETSC_MAX_PATH_LEN - 1, "%s", fieldname);
126: swarm->vec_field_set = PETSC_TRUE;
127: swarm->vec_field_bs = bs;
128: swarm->vec_field_nlocal = n;
129: DMSwarmRestoreField(dm, fieldname, &bs, &type, (void **)&array);
130: return 0;
131: }
133: /* requires DMSwarmDefineFieldVector has been called */
134: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm, Vec *vec)
135: {
136: DM_Swarm *swarm = (DM_Swarm *)dm->data;
137: Vec x;
138: char name[PETSC_MAX_PATH_LEN];
140: if (!swarm->issetup) DMSetUp(dm);
144: PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmField_%s", swarm->vec_field_name);
145: VecCreate(PetscObjectComm((PetscObject)dm), &x);
146: PetscObjectSetName((PetscObject)x, name);
147: VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE);
148: VecSetBlockSize(x, swarm->vec_field_bs);
149: VecSetDM(x, dm);
150: VecSetFromOptions(x);
151: VecSetDM(x, dm);
152: VecSetOperation(x, VECOP_VIEW, (void (*)(void))VecView_Swarm);
153: *vec = x;
154: return 0;
155: }
157: /* requires DMSwarmDefineFieldVector has been called */
158: PetscErrorCode DMCreateLocalVector_Swarm(DM dm, Vec *vec)
159: {
160: DM_Swarm *swarm = (DM_Swarm *)dm->data;
161: Vec x;
162: char name[PETSC_MAX_PATH_LEN];
164: if (!swarm->issetup) DMSetUp(dm);
168: PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmField_%s", swarm->vec_field_name);
169: VecCreate(PETSC_COMM_SELF, &x);
170: PetscObjectSetName((PetscObject)x, name);
171: VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE);
172: VecSetBlockSize(x, swarm->vec_field_bs);
173: VecSetDM(x, dm);
174: VecSetFromOptions(x);
175: *vec = x;
176: return 0;
177: }
179: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
180: {
181: DM_Swarm *swarm = (DM_Swarm *)dm->data;
182: DMSwarmDataField gfield;
183: PetscInt bs, nlocal, fid = -1, cfid = -2;
184: PetscBool flg;
186: /* check vector is an inplace array */
187: DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid);
188: PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg);
190: VecGetLocalSize(*vec, &nlocal);
191: VecGetBlockSize(*vec, &bs);
193: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
194: DMSwarmDataFieldRestoreAccess(gfield);
195: VecResetArray(*vec);
196: VecDestroy(vec);
197: return 0;
198: }
200: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
201: {
202: DM_Swarm *swarm = (DM_Swarm *)dm->data;
203: PetscDataType type;
204: PetscScalar *array;
205: PetscInt bs, n, fid;
206: char name[PETSC_MAX_PATH_LEN];
207: PetscMPIInt size;
208: PetscBool iscuda, iskokkos;
210: if (!swarm->issetup) DMSetUp(dm);
211: DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
212: DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array);
215: MPI_Comm_size(comm, &size);
216: PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos);
217: PetscStrcmp(dm->vectype, VECCUDA, &iscuda);
218: VecCreate(comm, vec);
219: VecSetSizes(*vec, n * bs, PETSC_DETERMINE);
220: VecSetBlockSize(*vec, bs);
221: if (iskokkos) VecSetType(*vec, VECKOKKOS);
222: else if (iscuda) VecSetType(*vec, VECCUDA);
223: else VecSetType(*vec, VECSTANDARD);
224: VecPlaceArray(*vec, array);
226: PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname);
227: PetscObjectSetName((PetscObject)*vec, name);
229: /* Set guard */
230: DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid);
231: PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid);
233: VecSetDM(*vec, dm);
234: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm);
235: return 0;
236: }
238: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
240: \hat f = \sum_i f_i \phi_i
242: and a particle function is given by
244: f = \sum_i w_i \delta(x - x_i)
246: then we want to require that
248: M \hat f = M_p f
250: where the particle mass matrix is given by
252: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
254: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
255: his integral. We allow this with the boolean flag.
256: */
257: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
258: {
259: const char *name = "Mass Matrix";
260: MPI_Comm comm;
261: PetscDS prob;
262: PetscSection fsection, globalFSection;
263: PetscHSetIJ ht;
264: PetscLayout rLayout, colLayout;
265: PetscInt *dnz, *onz;
266: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
267: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
268: PetscScalar *elemMat;
269: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
271: PetscObjectGetComm((PetscObject)mass, &comm);
272: DMGetCoordinateDim(dmf, &dim);
273: DMGetDS(dmf, &prob);
274: PetscDSGetNumFields(prob, &Nf);
275: PetscDSGetTotalDimension(prob, &totDim);
276: PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ);
277: DMGetLocalSection(dmf, &fsection);
278: DMGetGlobalSection(dmf, &globalFSection);
279: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
280: MatGetLocalSize(mass, &locRows, &locCols);
282: PetscLayoutCreate(comm, &colLayout);
283: PetscLayoutSetLocalSize(colLayout, locCols);
284: PetscLayoutSetBlockSize(colLayout, 1);
285: PetscLayoutSetUp(colLayout);
286: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
287: PetscLayoutDestroy(&colLayout);
289: PetscLayoutCreate(comm, &rLayout);
290: PetscLayoutSetLocalSize(rLayout, locRows);
291: PetscLayoutSetBlockSize(rLayout, 1);
292: PetscLayoutSetUp(rLayout);
293: PetscLayoutGetRange(rLayout, &rStart, NULL);
294: PetscLayoutDestroy(&rLayout);
296: PetscCalloc2(locRows, &dnz, locRows, &onz);
297: PetscHSetIJCreate(&ht);
299: PetscSynchronizedFlush(comm, NULL);
300: /* count non-zeros */
301: DMSwarmSortGetAccess(dmc);
302: for (field = 0; field < Nf; ++field) {
303: for (cell = cStart; cell < cEnd; ++cell) {
304: PetscInt c, i;
305: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
306: PetscInt numFIndices, numCIndices;
308: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
309: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
310: maxC = PetscMax(maxC, numCIndices);
311: {
312: PetscHashIJKey key;
313: PetscBool missing;
314: for (i = 0; i < numFIndices; ++i) {
315: key.j = findices[i]; /* global column (from Plex) */
316: if (key.j >= 0) {
317: /* Get indices for coarse elements */
318: for (c = 0; c < numCIndices; ++c) {
319: key.i = cindices[c] + rStart; /* global cols (from Swarm) */
320: if (key.i < 0) continue;
321: PetscHSetIJQueryAdd(ht, key, &missing);
322: if (missing) {
323: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
324: else ++onz[key.i - rStart];
325: } else SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
326: }
327: }
328: }
329: PetscFree(cindices);
330: }
331: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
332: }
333: }
334: PetscHSetIJDestroy(&ht);
335: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
336: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
337: PetscFree2(dnz, onz);
338: PetscMalloc3(maxC * totDim, &elemMat, maxC, &rowIDXs, maxC * dim, &xi);
339: for (field = 0; field < Nf; ++field) {
340: PetscTabulation Tcoarse;
341: PetscObject obj;
342: PetscReal *coords;
343: PetscInt Nc, i;
345: PetscDSGetDiscretization(prob, field, &obj);
346: PetscFEGetNumComponents((PetscFE)obj, &Nc);
348: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
349: for (cell = cStart; cell < cEnd; ++cell) {
350: PetscInt *findices, *cindices;
351: PetscInt numFIndices, numCIndices;
352: PetscInt p, c;
354: /* TODO: Use DMField instead of assuming affine */
355: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
356: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
357: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
358: for (p = 0; p < numCIndices; ++p) CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p] * dim], &xi[p * dim]);
359: PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse);
360: /* Get elemMat entries by multiplying by weight */
361: PetscArrayzero(elemMat, numCIndices * totDim);
362: for (i = 0; i < numFIndices; ++i) {
363: for (p = 0; p < numCIndices; ++p) {
364: for (c = 0; c < Nc; ++c) {
365: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
366: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
367: }
368: }
369: }
370: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
371: if (0) DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);
372: MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);
373: PetscFree(cindices);
374: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
375: PetscTabulationDestroy(&Tcoarse);
376: }
377: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
378: }
379: PetscFree3(elemMat, rowIDXs, xi);
380: DMSwarmSortRestoreAccess(dmc);
381: PetscFree3(v0, J, invJ);
382: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
383: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
384: return 0;
385: }
387: /* Returns empty matrix for use with SNES FD */
388: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
389: {
390: Vec field;
391: PetscInt size;
393: DMGetGlobalVector(sw, &field);
394: VecGetLocalSize(field, &size);
395: DMRestoreGlobalVector(sw, &field);
396: MatCreate(PETSC_COMM_WORLD, m);
397: MatSetFromOptions(*m);
398: MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size);
399: MatSeqAIJSetPreallocation(*m, 1, NULL);
400: MatZeroEntries(*m);
401: MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY);
402: MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY);
403: MatShift(*m, 1.0);
404: MatSetDM(*m, sw);
405: return 0;
406: }
408: /* FEM cols, Particle rows */
409: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
410: {
411: PetscSection gsf;
412: PetscInt m, n;
413: void *ctx;
415: DMGetGlobalSection(dmFine, &gsf);
416: PetscSectionGetConstrainedStorageSize(gsf, &m);
417: DMSwarmGetLocalSize(dmCoarse, &n);
418: MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass);
419: MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
420: MatSetType(*mass, dmCoarse->mattype);
421: DMGetApplicationContext(dmFine, &ctx);
423: DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
424: MatViewFromOptions(*mass, NULL, "-mass_mat_view");
425: return 0;
426: }
428: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
429: {
430: const char *name = "Mass Matrix Square";
431: MPI_Comm comm;
432: PetscDS prob;
433: PetscSection fsection, globalFSection;
434: PetscHSetIJ ht;
435: PetscLayout rLayout, colLayout;
436: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
437: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
438: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
439: PetscScalar *elemMat, *elemMatSq;
440: PetscInt cdim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
442: PetscObjectGetComm((PetscObject)mass, &comm);
443: DMGetCoordinateDim(dmf, &cdim);
444: DMGetDS(dmf, &prob);
445: PetscDSGetNumFields(prob, &Nf);
446: PetscDSGetTotalDimension(prob, &totDim);
447: PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ);
448: DMGetLocalSection(dmf, &fsection);
449: DMGetGlobalSection(dmf, &globalFSection);
450: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
451: MatGetLocalSize(mass, &locRows, &locCols);
453: PetscLayoutCreate(comm, &colLayout);
454: PetscLayoutSetLocalSize(colLayout, locCols);
455: PetscLayoutSetBlockSize(colLayout, 1);
456: PetscLayoutSetUp(colLayout);
457: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
458: PetscLayoutDestroy(&colLayout);
460: PetscLayoutCreate(comm, &rLayout);
461: PetscLayoutSetLocalSize(rLayout, locRows);
462: PetscLayoutSetBlockSize(rLayout, 1);
463: PetscLayoutSetUp(rLayout);
464: PetscLayoutGetRange(rLayout, &rStart, NULL);
465: PetscLayoutDestroy(&rLayout);
467: DMPlexGetDepth(dmf, &depth);
468: DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize);
469: maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
470: PetscMalloc1(maxAdjSize, &adj);
472: PetscCalloc2(locRows, &dnz, locRows, &onz);
473: PetscHSetIJCreate(&ht);
474: /* Count nonzeros
475: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
476: */
477: DMSwarmSortGetAccess(dmc);
478: for (cell = cStart; cell < cEnd; ++cell) {
479: PetscInt i;
480: PetscInt *cindices;
481: PetscInt numCIndices;
482: #if 0
483: PetscInt adjSize = maxAdjSize, a, j;
484: #endif
486: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
487: maxC = PetscMax(maxC, numCIndices);
488: /* Diagonal block */
489: for (i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
490: #if 0
491: /* Off-diagonal blocks */
492: DMPlexGetAdjacency(dmf, cell, &adjSize, &adj);
493: for (a = 0; a < adjSize; ++a) {
494: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
495: const PetscInt ncell = adj[a];
496: PetscInt *ncindices;
497: PetscInt numNCIndices;
499: DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices);
500: {
501: PetscHashIJKey key;
502: PetscBool missing;
504: for (i = 0; i < numCIndices; ++i) {
505: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
506: if (key.i < 0) continue;
507: for (j = 0; j < numNCIndices; ++j) {
508: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
509: if (key.j < 0) continue;
510: PetscHSetIJQueryAdd(ht, key, &missing);
511: if (missing) {
512: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
513: else ++onz[key.i - rStart];
514: }
515: }
516: }
517: }
518: PetscFree(ncindices);
519: }
520: }
521: #endif
522: PetscFree(cindices);
523: }
524: PetscHSetIJDestroy(&ht);
525: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
526: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
527: PetscFree2(dnz, onz);
528: PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi);
529: /* Fill in values
530: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
531: Start just by producing block diagonal
532: Could loop over adjacent cells
533: Produce neighboring element matrix
534: TODO Determine which columns and rows correspond to shared dual vector
535: Do MatMatMult with rectangular matrices
536: Insert block
537: */
538: for (field = 0; field < Nf; ++field) {
539: PetscTabulation Tcoarse;
540: PetscObject obj;
541: PetscReal *coords;
542: PetscInt Nc, i;
544: PetscDSGetDiscretization(prob, field, &obj);
545: PetscFEGetNumComponents((PetscFE)obj, &Nc);
547: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
548: for (cell = cStart; cell < cEnd; ++cell) {
549: PetscInt *findices, *cindices;
550: PetscInt numFIndices, numCIndices;
551: PetscInt p, c;
553: /* TODO: Use DMField instead of assuming affine */
554: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
555: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
556: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
557: for (p = 0; p < numCIndices; ++p) CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p] * cdim], &xi[p * cdim]);
558: PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse);
559: /* Get elemMat entries by multiplying by weight */
560: PetscArrayzero(elemMat, numCIndices * totDim);
561: for (i = 0; i < numFIndices; ++i) {
562: for (p = 0; p < numCIndices; ++p) {
563: for (c = 0; c < Nc; ++c) {
564: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
565: elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
566: }
567: }
568: }
569: PetscTabulationDestroy(&Tcoarse);
570: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
571: if (0) DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);
572: /* Block diagonal */
573: if (numCIndices) {
574: PetscBLASInt blasn, blask;
575: PetscScalar one = 1.0, zero = 0.0;
577: PetscBLASIntCast(numCIndices, &blasn);
578: PetscBLASIntCast(numFIndices, &blask);
579: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
580: }
581: MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES);
582: /* TODO Off-diagonal */
583: PetscFree(cindices);
584: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
585: }
586: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
587: }
588: PetscFree4(elemMat, elemMatSq, rowIDXs, xi);
589: PetscFree(adj);
590: DMSwarmSortRestoreAccess(dmc);
591: PetscFree3(v0, J, invJ);
592: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
593: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
594: return 0;
595: }
597: /*@C
598: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
600: Collective on dmCoarse
602: Input parameters:
603: + dmCoarse - a DMSwarm
604: - dmFine - a DMPlex
606: Output parameter:
607: . mass - the square of the particle mass matrix
609: Level: advanced
611: Notes:
612: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
613: future to compute the full normal equations.
615: .seealso: `DMCreateMassMatrix()`
616: @*/
617: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
618: {
619: PetscInt n;
620: void *ctx;
622: DMSwarmGetLocalSize(dmCoarse, &n);
623: MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass);
624: MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
625: MatSetType(*mass, dmCoarse->mattype);
626: DMGetApplicationContext(dmFine, &ctx);
628: DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
629: MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view");
630: return 0;
631: }
633: /*@C
634: DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
636: Collective on dm
638: Input parameters:
639: + dm - a DMSwarm
640: - fieldname - the textual name given to a registered field
642: Output parameter:
643: . vec - the vector
645: Level: beginner
647: Notes:
648: The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
650: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
651: @*/
652: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
653: {
654: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
656: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
657: return 0;
658: }
660: /*@C
661: DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
663: Collective on dm
665: Input parameters:
666: + dm - a DMSwarm
667: - fieldname - the textual name given to a registered field
669: Output parameter:
670: . vec - the vector
672: Level: beginner
674: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
675: @*/
676: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
677: {
678: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
679: return 0;
680: }
682: /*@C
683: DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
685: Collective on dm
687: Input parameters:
688: + dm - a DMSwarm
689: - fieldname - the textual name given to a registered field
691: Output parameter:
692: . vec - the vector
694: Level: beginner
696: Notes:
697: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
699: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
700: @*/
701: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
702: {
703: MPI_Comm comm = PETSC_COMM_SELF;
705: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
706: return 0;
707: }
709: /*@C
710: DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
712: Collective on dm
714: Input parameters:
715: + dm - a DMSwarm
716: - fieldname - the textual name given to a registered field
718: Output parameter:
719: . vec - the vector
721: Level: beginner
723: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
724: @*/
725: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
726: {
727: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
728: return 0;
729: }
731: /*@C
732: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
734: Collective on dm
736: Input parameter:
737: . dm - a DMSwarm
739: Level: beginner
741: Notes:
742: After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
744: .seealso: `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
745: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
746: @*/
747: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
748: {
749: DM_Swarm *swarm = (DM_Swarm *)dm->data;
751: if (!swarm->field_registration_initialized) {
752: swarm->field_registration_initialized = PETSC_TRUE;
753: DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64); /* unique identifer */
754: DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT); /* used for communication */
755: }
756: return 0;
757: }
759: /*@
760: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
762: Collective on dm
764: Input parameter:
765: . dm - a DMSwarm
767: Level: beginner
769: Notes:
770: After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
772: .seealso: `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
773: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
774: @*/
775: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
776: {
777: DM_Swarm *swarm = (DM_Swarm *)dm->data;
779: if (!swarm->field_registration_finalized) DMSwarmDataBucketFinalize(swarm->db);
780: swarm->field_registration_finalized = PETSC_TRUE;
781: return 0;
782: }
784: /*@
785: DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
787: Not collective
789: Input parameters:
790: + dm - a DMSwarm
791: . nlocal - the length of each registered field
792: - buffer - the length of the buffer used to efficient dynamic re-sizing
794: Level: beginner
796: .seealso: `DMSwarmGetLocalSize()`
797: @*/
798: PetscErrorCode DMSwarmSetLocalSizes(DM dm, PetscInt nlocal, PetscInt buffer)
799: {
800: DM_Swarm *swarm = (DM_Swarm *)dm->data;
802: PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0);
803: DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer);
804: PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0);
805: return 0;
806: }
808: /*@
809: DMSwarmSetCellDM - Attachs a DM to a DMSwarm
811: Collective on dm
813: Input parameters:
814: + dm - a DMSwarm
815: - dmcell - the DM to attach to the DMSwarm
817: Level: beginner
819: Notes:
820: The attached DM (dmcell) will be queried for point location and
821: neighbor MPI-rank information if DMSwarmMigrate() is called.
823: .seealso: `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
824: @*/
825: PetscErrorCode DMSwarmSetCellDM(DM dm, DM dmcell)
826: {
827: DM_Swarm *swarm = (DM_Swarm *)dm->data;
829: swarm->dmcell = dmcell;
830: return 0;
831: }
833: /*@
834: DMSwarmGetCellDM - Fetches the attached cell DM
836: Collective on dm
838: Input parameter:
839: . dm - a DMSwarm
841: Output parameter:
842: . dmcell - the DM which was attached to the DMSwarm
844: Level: beginner
846: .seealso: `DMSwarmSetCellDM()`
847: @*/
848: PetscErrorCode DMSwarmGetCellDM(DM dm, DM *dmcell)
849: {
850: DM_Swarm *swarm = (DM_Swarm *)dm->data;
852: *dmcell = swarm->dmcell;
853: return 0;
854: }
856: /*@
857: DMSwarmGetLocalSize - Retrives the local length of fields registered
859: Not collective
861: Input parameter:
862: . dm - a DMSwarm
864: Output parameter:
865: . nlocal - the length of each registered field
867: Level: beginner
869: .seealso: `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
870: @*/
871: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
872: {
873: DM_Swarm *swarm = (DM_Swarm *)dm->data;
875: DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL);
876: return 0;
877: }
879: /*@
880: DMSwarmGetSize - Retrives the total length of fields registered
882: Collective on dm
884: Input parameter:
885: . dm - a DMSwarm
887: Output parameter:
888: . n - the total length of each registered field
890: Level: beginner
892: Note:
893: This calls MPI_Allreduce upon each call (inefficient but safe)
895: .seealso: `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
896: @*/
897: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
898: {
899: DM_Swarm *swarm = (DM_Swarm *)dm->data;
900: PetscInt nlocal;
902: DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL);
903: MPI_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
904: return 0;
905: }
907: /*@C
908: DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
910: Collective on dm
912: Input parameters:
913: + dm - a DMSwarm
914: . fieldname - the textual name to identify this field
915: . blocksize - the number of each data type
916: - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
918: Level: beginner
920: Notes:
921: The textual name for each registered field must be unique.
923: .seealso: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
924: @*/
925: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
926: {
927: DM_Swarm *swarm = (DM_Swarm *)dm->data;
928: size_t size;
939: PetscDataTypeGetSize(type, &size);
940: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
941: DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL);
942: {
943: DMSwarmDataField gfield;
945: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
946: DMSwarmDataFieldSetBlockSize(gfield, blocksize);
947: }
948: swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
949: return 0;
950: }
952: /*@C
953: DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
955: Collective on dm
957: Input parameters:
958: + dm - a DMSwarm
959: . fieldname - the textual name to identify this field
960: - size - the size in bytes of the user struct of each data type
962: Level: beginner
964: Notes:
965: The textual name for each registered field must be unique.
967: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
968: @*/
969: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
970: {
971: DM_Swarm *swarm = (DM_Swarm *)dm->data;
973: DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL);
974: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
975: return 0;
976: }
978: /*@C
979: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
981: Collective on dm
983: Input parameters:
984: + dm - a DMSwarm
985: . fieldname - the textual name to identify this field
986: . size - the size in bytes of the user data type
987: - blocksize - the number of each data type
989: Level: beginner
991: Notes:
992: The textual name for each registered field must be unique.
994: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
995: @*/
996: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
997: {
998: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1000: DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL);
1001: {
1002: DMSwarmDataField gfield;
1004: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
1005: DMSwarmDataFieldSetBlockSize(gfield, blocksize);
1006: }
1007: swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1008: return 0;
1009: }
1011: /*@C
1012: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1014: Not collective
1016: Input parameters:
1017: + dm - a DMSwarm
1018: - fieldname - the textual name to identify this field
1020: Output parameters:
1021: + blocksize - the number of each data type
1022: . type - the data type
1023: - data - pointer to raw array
1025: Level: beginner
1027: Notes:
1028: The array must be returned using a matching call to DMSwarmRestoreField().
1030: .seealso: `DMSwarmRestoreField()`
1031: @*/
1032: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1033: {
1034: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1035: DMSwarmDataField gfield;
1037: if (!swarm->issetup) DMSetUp(dm);
1038: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
1039: DMSwarmDataFieldGetAccess(gfield);
1040: DMSwarmDataFieldGetEntries(gfield, data);
1041: if (blocksize) *blocksize = gfield->bs;
1042: if (type) *type = gfield->petsc_type;
1043: return 0;
1044: }
1046: /*@C
1047: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1049: Not collective
1051: Input parameters:
1052: + dm - a DMSwarm
1053: - fieldname - the textual name to identify this field
1055: Output parameters:
1056: + blocksize - the number of each data type
1057: . type - the data type
1058: - data - pointer to raw array
1060: Level: beginner
1062: Notes:
1063: The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
1065: .seealso: `DMSwarmGetField()`
1066: @*/
1067: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1068: {
1069: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1070: DMSwarmDataField gfield;
1072: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
1073: DMSwarmDataFieldRestoreAccess(gfield);
1074: if (data) *data = NULL;
1075: return 0;
1076: }
1078: /*@
1079: DMSwarmAddPoint - Add space for one new point in the DMSwarm
1081: Not collective
1083: Input parameter:
1084: . dm - a DMSwarm
1086: Level: beginner
1088: Notes:
1089: The new point will have all fields initialized to zero.
1091: .seealso: `DMSwarmAddNPoints()`
1092: @*/
1093: PetscErrorCode DMSwarmAddPoint(DM dm)
1094: {
1095: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1097: if (!swarm->issetup) DMSetUp(dm);
1098: PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0);
1099: DMSwarmDataBucketAddPoint(swarm->db);
1100: PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0);
1101: return 0;
1102: }
1104: /*@
1105: DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
1107: Not collective
1109: Input parameters:
1110: + dm - a DMSwarm
1111: - npoints - the number of new points to add
1113: Level: beginner
1115: Notes:
1116: The new point will have all fields initialized to zero.
1118: .seealso: `DMSwarmAddPoint()`
1119: @*/
1120: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1121: {
1122: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1123: PetscInt nlocal;
1125: PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0);
1126: DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL);
1127: nlocal = nlocal + npoints;
1128: DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
1129: PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0);
1130: return 0;
1131: }
1133: /*@
1134: DMSwarmRemovePoint - Remove the last point from the DMSwarm
1136: Not collective
1138: Input parameter:
1139: . dm - a DMSwarm
1141: Level: beginner
1143: .seealso: `DMSwarmRemovePointAtIndex()`
1144: @*/
1145: PetscErrorCode DMSwarmRemovePoint(DM dm)
1146: {
1147: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1149: PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0);
1150: DMSwarmDataBucketRemovePoint(swarm->db);
1151: PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0);
1152: return 0;
1153: }
1155: /*@
1156: DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
1158: Not collective
1160: Input parameters:
1161: + dm - a DMSwarm
1162: - idx - index of point to remove
1164: Level: beginner
1166: .seealso: `DMSwarmRemovePoint()`
1167: @*/
1168: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
1169: {
1170: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1172: PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0);
1173: DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx);
1174: PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0);
1175: return 0;
1176: }
1178: /*@
1179: DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
1181: Not collective
1183: Input parameters:
1184: + dm - a DMSwarm
1185: . pi - the index of the point to copy
1186: - pj - the point index where the copy should be located
1188: Level: beginner
1190: .seealso: `DMSwarmRemovePoint()`
1191: @*/
1192: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
1193: {
1194: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1196: if (!swarm->issetup) DMSetUp(dm);
1197: DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj);
1198: return 0;
1199: }
1201: PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
1202: {
1203: DMSwarmMigrate_Push_Basic(dm, remove_sent_points);
1204: return 0;
1205: }
1207: /*@
1208: DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
1210: Collective on dm
1212: Input parameters:
1213: + dm - the DMSwarm
1214: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1216: Notes:
1217: The DM will be modified to accommodate received points.
1218: If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
1219: Different styles of migration are supported. See DMSwarmSetMigrateType().
1221: Level: advanced
1223: .seealso: `DMSwarmSetMigrateType()`
1224: @*/
1225: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
1226: {
1227: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1229: PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0);
1230: switch (swarm->migrate_type) {
1231: case DMSWARM_MIGRATE_BASIC:
1232: DMSwarmMigrate_Basic(dm, remove_sent_points);
1233: break;
1234: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1235: DMSwarmMigrate_CellDMScatter(dm, remove_sent_points);
1236: break;
1237: case DMSWARM_MIGRATE_DMCELLEXACT:
1238: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1239: case DMSWARM_MIGRATE_USER:
1240: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
1241: default:
1242: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
1243: }
1244: PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0);
1245: DMClearGlobalVectors(dm);
1246: return 0;
1247: }
1249: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm, PetscInt *globalsize);
1251: /*
1252: DMSwarmCollectViewCreate
1254: * Applies a collection method and gathers point neighbour points into dm
1256: Notes:
1257: Users should call DMSwarmCollectViewDestroy() after
1258: they have finished computations associated with the collected points
1259: */
1261: /*@
1262: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1263: in neighbour ranks into the DMSwarm
1265: Collective on dm
1267: Input parameter:
1268: . dm - the DMSwarm
1270: Notes:
1271: Users should call DMSwarmCollectViewDestroy() after
1272: they have finished computations associated with the collected points
1273: Different collect methods are supported. See DMSwarmSetCollectType().
1275: Level: advanced
1277: .seealso: `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
1278: @*/
1279: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1280: {
1281: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1282: PetscInt ng;
1285: DMSwarmGetLocalSize(dm, &ng);
1286: switch (swarm->collect_type) {
1287: case DMSWARM_COLLECT_BASIC:
1288: DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng);
1289: break;
1290: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1291: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1292: case DMSWARM_COLLECT_GENERAL:
1293: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
1294: default:
1295: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
1296: }
1297: swarm->collect_view_active = PETSC_TRUE;
1298: swarm->collect_view_reset_nlocal = ng;
1299: return 0;
1300: }
1302: /*@
1303: DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
1305: Collective on dm
1307: Input parameters:
1308: . dm - the DMSwarm
1310: Notes:
1311: Users should call DMSwarmCollectViewCreate() before this function is called.
1313: Level: advanced
1315: .seealso: `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
1316: @*/
1317: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1318: {
1319: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1322: DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1);
1323: swarm->collect_view_active = PETSC_FALSE;
1324: return 0;
1325: }
1327: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1328: {
1329: PetscInt dim;
1331: DMSwarmSetNumSpecies(dm, 1);
1332: DMGetDimension(dm, &dim);
1335: DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_coor, dim, PETSC_DOUBLE);
1336: DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_cellid, 1, PETSC_INT);
1337: return 0;
1338: }
1340: /*@
1341: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
1343: Collective on dm
1345: Input parameters:
1346: + dm - the DMSwarm
1347: - Npc - The number of particles per cell in the cell DM
1349: Notes:
1350: The user must use DMSwarmSetCellDM() to set the cell DM first. The particles are placed randomly inside each cell. If only
1351: one particle is in each cell, it is placed at the centroid.
1353: Level: intermediate
1355: .seealso: `DMSwarmSetCellDM()`
1356: @*/
1357: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1358: {
1359: DM cdm;
1360: PetscRandom rnd;
1361: DMPolytopeType ct;
1362: PetscBool simplex;
1363: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1364: PetscInt dim, d, cStart, cEnd, c, p;
1367: PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd);
1368: PetscRandomSetInterval(rnd, -1.0, 1.0);
1369: PetscRandomSetType(rnd, PETSCRAND48);
1371: DMSwarmGetCellDM(dm, &cdm);
1372: DMGetDimension(cdm, &dim);
1373: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
1374: DMPlexGetCellType(cdm, cStart, &ct);
1375: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1377: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ);
1378: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1379: DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
1380: for (c = cStart; c < cEnd; ++c) {
1381: if (Npc == 1) {
1382: DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL);
1383: for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
1384: } else {
1385: DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ); /* affine */
1386: for (p = 0; p < Npc; ++p) {
1387: const PetscInt n = c * Npc + p;
1388: PetscReal sum = 0.0, refcoords[3];
1390: for (d = 0; d < dim; ++d) {
1391: PetscRandomGetValueReal(rnd, &refcoords[d]);
1392: sum += refcoords[d];
1393: }
1394: if (simplex && sum > 0.0)
1395: for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
1396: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
1397: }
1398: }
1399: }
1400: DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
1401: PetscFree5(centroid, xi0, v0, J, invJ);
1402: PetscRandomDestroy(&rnd);
1403: return 0;
1404: }
1406: /*@
1407: DMSwarmSetType - Set particular flavor of DMSwarm
1409: Collective on dm
1411: Input parameters:
1412: + dm - the DMSwarm
1413: - stype - the DMSwarm type (e.g. DMSWARM_PIC)
1415: Level: advanced
1417: .seealso: `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
1418: @*/
1419: PetscErrorCode DMSwarmSetType(DM dm, DMSwarmType stype)
1420: {
1421: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1423: swarm->swarm_type = stype;
1424: if (swarm->swarm_type == DMSWARM_PIC) DMSwarmSetUpPIC(dm);
1425: return 0;
1426: }
1428: PetscErrorCode DMSetup_Swarm(DM dm)
1429: {
1430: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1431: PetscMPIInt rank;
1432: PetscInt p, npoints, *rankval;
1434: if (swarm->issetup) return 0;
1435: swarm->issetup = PETSC_TRUE;
1437: if (swarm->swarm_type == DMSWARM_PIC) {
1438: /* check dmcell exists */
1441: if (swarm->dmcell->ops->locatepointssubdomain) {
1442: /* check methods exists for exact ownership identificiation */
1443: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1444: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1445: } else {
1446: /* check methods exist for point location AND rank neighbor identification */
1447: if (swarm->dmcell->ops->locatepoints) {
1448: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1449: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1451: if (swarm->dmcell->ops->getneighbors) {
1452: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1453: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1455: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1456: }
1457: }
1459: DMSwarmFinalizeFieldRegister(dm);
1461: /* check some fields were registered */
1464: /* check local sizes were set */
1467: /* initialize values in pid and rank placeholders */
1468: /* TODO: [pid - use MPI_Scan] */
1469: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1470: DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
1471: DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
1472: for (p = 0; p < npoints; p++) rankval[p] = (PetscInt)rank;
1473: DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
1474: return 0;
1475: }
1477: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1479: PetscErrorCode DMDestroy_Swarm(DM dm)
1480: {
1481: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1483: if (--swarm->refct > 0) return 0;
1484: DMSwarmDataBucketDestroy(&swarm->db);
1485: if (swarm->sort_context) DMSwarmSortDestroy(&swarm->sort_context);
1486: PetscFree(swarm);
1487: return 0;
1488: }
1490: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1491: {
1492: DM cdm;
1493: PetscDraw draw;
1494: PetscReal *coords, oldPause, radius = 0.01;
1495: PetscInt Np, p, bs;
1497: PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL);
1498: PetscViewerDrawGetDraw(viewer, 0, &draw);
1499: DMSwarmGetCellDM(dm, &cdm);
1500: PetscDrawGetPause(draw, &oldPause);
1501: PetscDrawSetPause(draw, 0.0);
1502: DMView(cdm, viewer);
1503: PetscDrawSetPause(draw, oldPause);
1505: DMSwarmGetLocalSize(dm, &Np);
1506: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords);
1507: for (p = 0; p < Np; ++p) {
1508: const PetscInt i = p * bs;
1510: PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE);
1511: }
1512: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords);
1513: PetscDrawFlush(draw);
1514: PetscDrawPause(draw);
1515: PetscDrawSave(draw);
1516: return 0;
1517: }
1519: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
1520: {
1521: PetscViewerFormat format;
1522: PetscInt *sizes;
1523: PetscInt dim, Np, maxSize = 17;
1524: MPI_Comm comm;
1525: PetscMPIInt rank, size, p;
1526: const char *name;
1528: PetscViewerGetFormat(viewer, &format);
1529: DMGetDimension(dm, &dim);
1530: DMSwarmGetLocalSize(dm, &Np);
1531: PetscObjectGetComm((PetscObject)dm, &comm);
1532: MPI_Comm_rank(comm, &rank);
1533: MPI_Comm_size(comm, &size);
1534: PetscObjectGetName((PetscObject)dm, &name);
1535: if (name) PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s");
1536: else PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s");
1537: if (size < maxSize) PetscCalloc1(size, &sizes);
1538: else PetscCalloc1(3, &sizes);
1539: if (size < maxSize) {
1540: MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
1541: PetscViewerASCIIPrintf(viewer, " Number of particles per rank:");
1542: for (p = 0; p < size; ++p) {
1543: if (rank == 0) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]);
1544: }
1545: } else {
1546: PetscInt locMinMax[2] = {Np, Np};
1548: PetscGlobalMinMaxInt(comm, locMinMax, sizes);
1549: PetscViewerASCIIPrintf(viewer, " Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]);
1550: }
1551: PetscViewerASCIIPrintf(viewer, "\n");
1552: PetscFree(sizes);
1553: if (format == PETSC_VIEWER_ASCII_INFO) {
1554: PetscInt *cell;
1556: PetscViewerASCIIPrintf(viewer, " Cells containing each particle:\n");
1557: PetscViewerASCIIPushSynchronized(viewer);
1558: DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell);
1559: for (p = 0; p < Np; ++p) PetscViewerASCIISynchronizedPrintf(viewer, " p%d: %" PetscInt_FMT "\n", p, cell[p]);
1560: PetscViewerFlush(viewer);
1561: PetscViewerASCIIPopSynchronized(viewer);
1562: DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell);
1563: }
1564: return 0;
1565: }
1567: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1568: {
1569: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1570: PetscBool iascii, ibinary, isvtk, isdraw;
1571: #if defined(PETSC_HAVE_HDF5)
1572: PetscBool ishdf5;
1573: #endif
1577: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1578: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary);
1579: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1580: #if defined(PETSC_HAVE_HDF5)
1581: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1582: #endif
1583: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1584: if (iascii) {
1585: PetscViewerFormat format;
1587: PetscViewerGetFormat(viewer, &format);
1588: switch (format) {
1589: case PETSC_VIEWER_ASCII_INFO_DETAIL:
1590: DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT);
1591: break;
1592: default:
1593: DMView_Swarm_Ascii(dm, viewer);
1594: }
1595: } else {
1596: #if defined(PETSC_HAVE_HDF5)
1597: if (ishdf5) DMSwarmView_HDF5(dm, viewer);
1598: #endif
1599: if (isdraw) DMSwarmView_Draw(dm, viewer);
1600: }
1601: return 0;
1602: }
1604: /*@C
1605: DMSwarmGetCellSwarm - Extracts a single cell from the DMSwarm object, returns it as a single cell DMSwarm.
1606: 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.
1608: Important: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
1610: Noncollective
1612: Input parameters:
1613: + sw - the DMSwarm
1614: . cellID - the integer id of the cell to be extracted and filtered
1615: - cellswarm - The DMSwarm to receive the cell
1617: Level: beginner
1619: Notes:
1620: This presently only supports DMSWARM_PIC type
1622: Should be restored with DMSwarmRestoreCellSwarm()
1624: .seealso: `DMSwarmRestoreCellSwarm()`
1625: @*/
1626: PETSC_EXTERN PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1627: {
1628: DM_Swarm *original = (DM_Swarm *)sw->data;
1629: DMLabel label;
1630: DM dmc, subdmc;
1631: PetscInt *pids, particles, dim;
1633: /* Configure new swarm */
1634: DMSetType(cellswarm, DMSWARM);
1635: DMGetDimension(sw, &dim);
1636: DMSetDimension(cellswarm, dim);
1637: DMSwarmSetType(cellswarm, DMSWARM_PIC);
1638: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
1639: DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db);
1640: DMSwarmSortGetAccess(sw);
1641: DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles);
1642: DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids);
1643: DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db);
1644: DMSwarmSortRestoreAccess(sw);
1645: PetscFree(pids);
1646: DMSwarmGetCellDM(sw, &dmc);
1647: DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label);
1648: DMAddLabel(dmc, label);
1649: DMLabelSetValue(label, cellID, 1);
1650: DMPlexFilter(dmc, label, 1, &subdmc);
1651: DMSwarmSetCellDM(cellswarm, subdmc);
1652: DMLabelDestroy(&label);
1653: return 0;
1654: }
1656: /*@C
1657: DMSwarmRestoreCellSwarm - Restores a DMSwarm object obtained with DMSwarmGetCellSwarm(). All fields are copied back into the parent swarm.
1659: Noncollective
1661: Input parameters:
1662: + sw - the parent DMSwarm
1663: . cellID - the integer id of the cell to be copied back into the parent swarm
1664: - cellswarm - the cell swarm object
1666: Level: beginner
1668: Note:
1669: This only supports DMSWARM_PIC types of DMSwarms
1671: .seealso: `DMSwarmGetCellSwarm()`
1672: @*/
1673: PETSC_EXTERN PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1674: {
1675: DM dmc;
1676: PetscInt *pids, particles, p;
1678: DMSwarmSortGetAccess(sw);
1679: DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids);
1680: DMSwarmSortRestoreAccess(sw);
1681: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
1682: for (p = 0; p < particles; ++p) DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]);
1683: /* Free memory, destroy cell dm */
1684: DMSwarmGetCellDM(cellswarm, &dmc);
1685: DMDestroy(&dmc);
1686: PetscFree(pids);
1687: return 0;
1688: }
1690: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
1692: static PetscErrorCode DMInitialize_Swarm(DM sw)
1693: {
1694: sw->dim = 0;
1695: sw->ops->view = DMView_Swarm;
1696: sw->ops->load = NULL;
1697: sw->ops->setfromoptions = NULL;
1698: sw->ops->clone = DMClone_Swarm;
1699: sw->ops->setup = DMSetup_Swarm;
1700: sw->ops->createlocalsection = NULL;
1701: sw->ops->createdefaultconstraints = NULL;
1702: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1703: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
1704: sw->ops->getlocaltoglobalmapping = NULL;
1705: sw->ops->createfieldis = NULL;
1706: sw->ops->createcoordinatedm = NULL;
1707: sw->ops->getcoloring = NULL;
1708: sw->ops->creatematrix = DMCreateMatrix_Swarm;
1709: sw->ops->createinterpolation = NULL;
1710: sw->ops->createinjection = NULL;
1711: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1712: sw->ops->refine = NULL;
1713: sw->ops->coarsen = NULL;
1714: sw->ops->refinehierarchy = NULL;
1715: sw->ops->coarsenhierarchy = NULL;
1716: sw->ops->globaltolocalbegin = NULL;
1717: sw->ops->globaltolocalend = NULL;
1718: sw->ops->localtoglobalbegin = NULL;
1719: sw->ops->localtoglobalend = NULL;
1720: sw->ops->destroy = DMDestroy_Swarm;
1721: sw->ops->createsubdm = NULL;
1722: sw->ops->getdimpoints = NULL;
1723: sw->ops->locatepoints = NULL;
1724: return 0;
1725: }
1727: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
1728: {
1729: DM_Swarm *swarm = (DM_Swarm *)dm->data;
1731: swarm->refct++;
1732: (*newdm)->data = swarm;
1733: PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM);
1734: DMInitialize_Swarm(*newdm);
1735: (*newdm)->dim = dm->dim;
1736: return 0;
1737: }
1739: /*MC
1741: DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1742: This implementation was designed for particle methods in which the underlying
1743: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1745: User data can be represented by DMSwarm through a registering "fields".
1746: To register a field, the user must provide;
1747: (a) a unique name;
1748: (b) the data type (or size in bytes);
1749: (c) the block size of the data.
1751: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1752: on a set of particles. Then the following code could be used
1754: $ DMSwarmInitializeFieldRegister(dm)
1755: $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1756: $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1757: $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1758: $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1759: $ DMSwarmFinalizeFieldRegister(dm)
1761: The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1762: The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1764: To support particle methods, "migration" techniques are provided. These methods migrate data
1765: between ranks.
1767: DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1768: As a DMSwarm may internally define and store values of different data types,
1769: before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1770: fields should be used to define a Vec object via
1771: DMSwarmVectorDefineField()
1772: The specified field can be changed at any time - thereby permitting vectors
1773: compatible with different fields to be created.
1775: A dual representation of fields in the DMSwarm and a Vec object is permitted via
1776: DMSwarmCreateGlobalVectorFromField()
1777: Here the data defining the field in the DMSwarm is shared with a Vec.
1778: This is inherently unsafe if you alter the size of the field at any time between
1779: calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1780: If the local size of the DMSwarm does not match the local size of the global vector
1781: when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1783: Additional high-level support is provided for Particle-In-Cell methods.
1784: Please refer to the man page for DMSwarmSetType().
1786: Level: beginner
1788: .seealso: `DMType`, `DMCreate()`, `DMSetType()`
1789: M*/
1790: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1791: {
1792: DM_Swarm *swarm;
1795: PetscNew(&swarm);
1796: dm->data = swarm;
1797: DMSwarmDataBucketCreate(&swarm->db);
1798: DMSwarmInitializeFieldRegister(dm);
1799: swarm->refct = 1;
1800: swarm->vec_field_set = PETSC_FALSE;
1801: swarm->issetup = PETSC_FALSE;
1802: swarm->swarm_type = DMSWARM_BASIC;
1803: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1804: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1805: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1806: swarm->dmcell = NULL;
1807: swarm->collect_view_active = PETSC_FALSE;
1808: swarm->collect_view_reset_nlocal = -1;
1809: DMInitialize_Swarm(dm);
1810: if (SwarmDataFieldId == -1) PetscObjectComposedDataRegister(&SwarmDataFieldId);
1811: return 0;
1812: }