1 #define PETSCDM_DLL
2 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/
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
11 PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
12 PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
13 PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
14
15 const char* DMSwarmTypeNames[] = { "basic", "pic", NULL };
16 const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", NULL };
17 const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", NULL };
18 const char* DMSwarmPICLayoutTypeNames[] = { "regular", "gauss", "subdivision", NULL };
19
20 const char DMSwarmField_pid[] = "DMSwarm_pid";
21 const char DMSwarmField_rank[] = "DMSwarm_rank";
22 const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
23 const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
24
25 PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
26 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
27
28 #if defined(PETSC_HAVE_HDF5)
29 #include <petscviewerhdf5.h>
30
VecView_Swarm_HDF5_Internal(Vec v,PetscViewer viewer)31 PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
32 {
33 DM dm;
34 PetscReal seqval;
35 PetscInt seqnum, bs;
36 PetscBool isseq;
37 PetscErrorCode ierr;
38
39 PetscFunctionBegin;
40 ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
41 ierr = VecGetBlockSize(v, &bs);CHKERRQ(ierr);
42 ierr = PetscViewerHDF5PushGroup(viewer, "/particle_fields");CHKERRQ(ierr);
43 ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr);
44 ierr = DMGetOutputSequenceNumber(dm, &seqnum, &seqval);CHKERRQ(ierr);
45 ierr = PetscViewerHDF5SetTimestep(viewer, seqnum);CHKERRQ(ierr);
46 /* ierr = DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);CHKERRQ(ierr); */
47 if (isseq) {ierr = VecView_Seq(v, viewer);CHKERRQ(ierr);}
48 else {ierr = VecView_MPI(v, viewer);CHKERRQ(ierr);}
49 ierr = PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) v, "Nc", PETSC_INT, (void *) &bs);CHKERRQ(ierr);
50 ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
51 PetscFunctionReturn(0);
52 }
53
DMSwarmView_HDF5(DM dm,PetscViewer viewer)54 PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
55 {
56 Vec coordinates;
57 PetscInt Np;
58 PetscBool isseq;
59 PetscErrorCode ierr;
60
61 PetscFunctionBegin;
62 ierr = DMSwarmGetSize(dm, &Np);CHKERRQ(ierr);
63 ierr = DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);CHKERRQ(ierr);
64 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
65 ierr = PetscViewerHDF5PushGroup(viewer, "/particles");CHKERRQ(ierr);
66 ierr = PetscObjectTypeCompare((PetscObject) coordinates, VECSEQ, &isseq);CHKERRQ(ierr);
67 if (isseq) {ierr = VecView_Seq(coordinates, viewer);CHKERRQ(ierr);}
68 else {ierr = VecView_MPI(coordinates, viewer);CHKERRQ(ierr);}
69 ierr = PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) coordinates, "Np", PETSC_INT, (void *) &Np);CHKERRQ(ierr);
70 ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
71 ierr = DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);CHKERRQ(ierr);
72 PetscFunctionReturn(0);
73 }
74 #endif
75
VecView_Swarm(Vec v,PetscViewer viewer)76 PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
77 {
78 DM dm;
79 PetscBool ishdf5;
80 PetscErrorCode ierr;
81
82 PetscFunctionBegin;
83 ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
84 if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
85 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr);
86 if (ishdf5) {
87 #if defined(PETSC_HAVE_HDF5)
88 ierr = VecView_Swarm_HDF5_Internal(v, viewer);CHKERRQ(ierr);
89 #else
90 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
91 #endif
92 } else {
93 PetscBool isseq;
94
95 ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr);
96 if (isseq) {ierr = VecView_Seq(v, viewer);CHKERRQ(ierr);}
97 else {ierr = VecView_MPI(v, viewer);CHKERRQ(ierr);}
98 }
99 PetscFunctionReturn(0);
100 }
101
102 /*@C
103 DMSwarmVectorDefineField - Sets the field from which to define a Vec object
104 when DMCreateLocalVector(), or DMCreateGlobalVector() is called
105
106 Collective on dm
107
108 Input parameters:
109 + dm - a DMSwarm
110 - fieldname - the textual name given to a registered field
111
112 Level: beginner
113
114 Notes:
115
116 The field with name fieldname must be defined as having a data type of PetscScalar.
117
118 This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
119 Mutiple calls to DMSwarmVectorDefineField() are permitted.
120
121 .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
122 @*/
DMSwarmVectorDefineField(DM dm,const char fieldname[])123 PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
124 {
125 DM_Swarm *swarm = (DM_Swarm*)dm->data;
126 PetscErrorCode ierr;
127 PetscInt bs,n;
128 PetscScalar *array;
129 PetscDataType type;
130
131 PetscFunctionBegin;
132 if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
133 ierr = DMSwarmDataBucketGetSizes(swarm->db,&n,NULL,NULL);CHKERRQ(ierr);
134 ierr = DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr);
135
136 /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
137 if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
138 ierr = PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);CHKERRQ(ierr);
139 swarm->vec_field_set = PETSC_TRUE;
140 swarm->vec_field_bs = bs;
141 swarm->vec_field_nlocal = n;
142 ierr = DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr);
143 PetscFunctionReturn(0);
144 }
145
146 /* requires DMSwarmDefineFieldVector has been called */
DMCreateGlobalVector_Swarm(DM dm,Vec * vec)147 PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
148 {
149 DM_Swarm *swarm = (DM_Swarm*)dm->data;
150 PetscErrorCode ierr;
151 Vec x;
152 char name[PETSC_MAX_PATH_LEN];
153
154 PetscFunctionBegin;
155 if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
156 if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
157 if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
158
159 ierr = PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);CHKERRQ(ierr);
160 ierr = VecCreate(PetscObjectComm((PetscObject)dm),&x);CHKERRQ(ierr);
161 ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr);
162 ierr = VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);CHKERRQ(ierr);
163 ierr = VecSetBlockSize(x,swarm->vec_field_bs);CHKERRQ(ierr);
164 ierr = VecSetDM(x,dm);CHKERRQ(ierr);
165 ierr = VecSetFromOptions(x);CHKERRQ(ierr);
166 ierr = VecSetDM(x, dm);CHKERRQ(ierr);
167 ierr = VecSetOperation(x, VECOP_VIEW, (void (*)(void)) VecView_Swarm);CHKERRQ(ierr);
168 *vec = x;
169 PetscFunctionReturn(0);
170 }
171
172 /* requires DMSwarmDefineFieldVector has been called */
DMCreateLocalVector_Swarm(DM dm,Vec * vec)173 PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
174 {
175 DM_Swarm *swarm = (DM_Swarm*)dm->data;
176 PetscErrorCode ierr;
177 Vec x;
178 char name[PETSC_MAX_PATH_LEN];
179
180 PetscFunctionBegin;
181 if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
182 if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
183 if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
184
185 ierr = PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);CHKERRQ(ierr);
186 ierr = VecCreate(PETSC_COMM_SELF,&x);CHKERRQ(ierr);
187 ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr);
188 ierr = VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);CHKERRQ(ierr);
189 ierr = VecSetBlockSize(x,swarm->vec_field_bs);CHKERRQ(ierr);
190 ierr = VecSetDM(x,dm);CHKERRQ(ierr);
191 ierr = VecSetFromOptions(x);CHKERRQ(ierr);
192 *vec = x;
193 PetscFunctionReturn(0);
194 }
195
DMSwarmDestroyVectorFromField_Private(DM dm,const char fieldname[],Vec * vec)196 static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
197 {
198 DM_Swarm *swarm = (DM_Swarm *) dm->data;
199 DMSwarmDataField gfield;
200 void (*fptr)(void);
201 PetscInt bs, nlocal;
202 char name[PETSC_MAX_PATH_LEN];
203 PetscErrorCode ierr;
204
205 PetscFunctionBegin;
206 ierr = VecGetLocalSize(*vec, &nlocal);CHKERRQ(ierr);
207 ierr = VecGetBlockSize(*vec, &bs);CHKERRQ(ierr);
208 if (nlocal/bs != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid"); /* Stale data */
209 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);CHKERRQ(ierr);
210 /* check vector is an inplace array */
211 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);CHKERRQ(ierr);
212 ierr = PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);CHKERRQ(ierr);
213 if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
214 ierr = DMSwarmDataFieldRestoreAccess(gfield);CHKERRQ(ierr);
215 ierr = VecDestroy(vec);CHKERRQ(ierr);
216 PetscFunctionReturn(0);
217 }
218
DMSwarmCreateVectorFromField_Private(DM dm,const char fieldname[],MPI_Comm comm,Vec * vec)219 static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
220 {
221 DM_Swarm *swarm = (DM_Swarm *) dm->data;
222 PetscDataType type;
223 PetscScalar *array;
224 PetscInt bs, n;
225 char name[PETSC_MAX_PATH_LEN];
226 PetscMPIInt size;
227 PetscErrorCode ierr;
228
229 PetscFunctionBegin;
230 if (!swarm->issetup) {ierr = DMSetUp(dm);CHKERRQ(ierr);}
231 ierr = DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);CHKERRQ(ierr);
232 ierr = DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);CHKERRQ(ierr);
233 if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
234
235 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
236 if (size == 1) {
237 ierr = VecCreateSeqWithArray(comm, bs, n*bs, array, vec);CHKERRQ(ierr);
238 } else {
239 ierr = VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);CHKERRQ(ierr);
240 }
241 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);CHKERRQ(ierr);
242 ierr = PetscObjectSetName((PetscObject) *vec, name);CHKERRQ(ierr);
243
244 /* Set guard */
245 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);CHKERRQ(ierr);
246 ierr = PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);CHKERRQ(ierr);
247
248 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr);
249 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Swarm);CHKERRQ(ierr);
250 PetscFunctionReturn(0);
251 }
252
253 /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
254
255 \hat f = \sum_i f_i \phi_i
256
257 and a particle function is given by
258
259 f = \sum_i w_i \delta(x - x_i)
260
261 then we want to require that
262
263 M \hat f = M_p f
264
265 where the particle mass matrix is given by
266
267 (M_p)_{ij} = \int \phi_i \delta(x - x_j)
268
269 The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
270 his integral. We allow this with the boolean flag.
271 */
DMSwarmComputeMassMatrix_Private(DM dmc,DM dmf,Mat mass,PetscBool useDeltaFunction,void * ctx)272 static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
273 {
274 const char *name = "Mass Matrix";
275 MPI_Comm comm;
276 PetscDS prob;
277 PetscSection fsection, globalFSection;
278 PetscHSetIJ ht;
279 PetscLayout rLayout, colLayout;
280 PetscInt *dnz, *onz;
281 PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
282 PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
283 PetscScalar *elemMat;
284 PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
285 PetscErrorCode ierr;
286
287 PetscFunctionBegin;
288 ierr = PetscObjectGetComm((PetscObject) mass, &comm);CHKERRQ(ierr);
289 ierr = DMGetCoordinateDim(dmf, &dim);CHKERRQ(ierr);
290 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr);
291 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr);
292 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
293 ierr = PetscMalloc3(dim, &v0, dim*dim, &J, dim*dim,&invJ);CHKERRQ(ierr);
294 ierr = DMGetLocalSection(dmf, &fsection);CHKERRQ(ierr);
295 ierr = DMGetGlobalSection(dmf, &globalFSection);CHKERRQ(ierr);
296 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr);
297 ierr = MatGetLocalSize(mass, &locRows, &locCols);CHKERRQ(ierr);
298
299 ierr = PetscLayoutCreate(comm, &colLayout);CHKERRQ(ierr);
300 ierr = PetscLayoutSetLocalSize(colLayout, locCols);CHKERRQ(ierr);
301 ierr = PetscLayoutSetBlockSize(colLayout, 1);CHKERRQ(ierr);
302 ierr = PetscLayoutSetUp(colLayout);CHKERRQ(ierr);
303 ierr = PetscLayoutGetRange(colLayout, &colStart, &colEnd);CHKERRQ(ierr);
304 ierr = PetscLayoutDestroy(&colLayout);CHKERRQ(ierr);
305
306 ierr = PetscLayoutCreate(comm, &rLayout);CHKERRQ(ierr);
307 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
308 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
309 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
310 ierr = PetscLayoutGetRange(rLayout, &rStart, NULL);CHKERRQ(ierr);
311 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
312
313 ierr = PetscCalloc2(locRows, &dnz, locRows, &onz);CHKERRQ(ierr);
314 ierr = PetscHSetIJCreate(&ht);CHKERRQ(ierr);
315
316 ierr = PetscSynchronizedFlush(comm, NULL);CHKERRQ(ierr);
317 /* count non-zeros */
318 ierr = DMSwarmSortGetAccess(dmc);CHKERRQ(ierr);
319 for (field = 0; field < Nf; ++field) {
320 for (cell = cStart; cell < cEnd; ++cell) {
321 PetscInt c, i;
322 PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
323 PetscInt numFIndices, numCIndices;
324
325 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr);
326 ierr = DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);CHKERRQ(ierr);
327 maxC = PetscMax(maxC, numCIndices);
328 {
329 PetscHashIJKey key;
330 PetscBool missing;
331 for (i = 0; i < numFIndices; ++i) {
332 key.j = findices[i]; /* global column (from Plex) */
333 if (key.j >= 0) {
334 /* Get indices for coarse elements */
335 for (c = 0; c < numCIndices; ++c) {
336 key.i = cindices[c] + rStart; /* global cols (from Swarm) */
337 if (key.i < 0) continue;
338 ierr = PetscHSetIJQueryAdd(ht, key, &missing);CHKERRQ(ierr);
339 if (missing) {
340 if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
341 else ++onz[key.i - rStart];
342 } else SETERRQ2(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Set new value at %D,%D", key.i, key.j);
343 }
344 }
345 }
346 ierr = PetscFree(cindices);CHKERRQ(ierr);
347 }
348 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr);
349 }
350 }
351 ierr = PetscHSetIJDestroy(&ht);CHKERRQ(ierr);
352 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr);
353 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);CHKERRQ(ierr);
354 ierr = PetscFree2(dnz, onz);CHKERRQ(ierr);
355 ierr = PetscMalloc3(maxC*totDim, &elemMat, maxC, &rowIDXs, maxC*dim, &xi);CHKERRQ(ierr);
356 for (field = 0; field < Nf; ++field) {
357 PetscTabulation Tcoarse;
358 PetscObject obj;
359 PetscReal *coords;
360 PetscInt Nc, i;
361
362 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr);
363 ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);
364 if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
365 ierr = DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr);
366 for (cell = cStart; cell < cEnd; ++cell) {
367 PetscInt *findices , *cindices;
368 PetscInt numFIndices, numCIndices;
369 PetscInt p, c;
370
371 /* TODO: Use DMField instead of assuming affine */
372 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
373 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr);
374 ierr = DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);CHKERRQ(ierr);
375 for (p = 0; p < numCIndices; ++p) {
376 CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p]*dim], &xi[p*dim]);
377 }
378 ierr = PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);CHKERRQ(ierr);
379 /* Get elemMat entries by multiplying by weight */
380 ierr = PetscArrayzero(elemMat, numCIndices*totDim);CHKERRQ(ierr);
381 for (i = 0; i < numFIndices; ++i) {
382 for (p = 0; p < numCIndices; ++p) {
383 for (c = 0; c < Nc; ++c) {
384 /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
385 elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
386 }
387 }
388 }
389 for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
390 if (0) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);}
391 ierr = MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);CHKERRQ(ierr);
392 ierr = PetscFree(cindices);CHKERRQ(ierr);
393 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr);
394 ierr = PetscTabulationDestroy(&Tcoarse);CHKERRQ(ierr);
395 }
396 ierr = DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr);
397 }
398 ierr = PetscFree3(elemMat, rowIDXs, xi);CHKERRQ(ierr);
399 ierr = DMSwarmSortRestoreAccess(dmc);CHKERRQ(ierr);
400 ierr = PetscFree3(v0, J, invJ);CHKERRQ(ierr);
401 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
402 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
403 PetscFunctionReturn(0);
404 }
405
406 /* FEM cols, Particle rows */
DMCreateMassMatrix_Swarm(DM dmCoarse,DM dmFine,Mat * mass)407 static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
408 {
409 PetscSection gsf;
410 PetscInt m, n;
411 void *ctx;
412 PetscErrorCode ierr;
413
414 PetscFunctionBegin;
415 ierr = DMGetGlobalSection(dmFine, &gsf);CHKERRQ(ierr);
416 ierr = PetscSectionGetConstrainedStorageSize(gsf, &m);CHKERRQ(ierr);
417 ierr = DMSwarmGetLocalSize(dmCoarse, &n);CHKERRQ(ierr);
418 ierr = MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);CHKERRQ(ierr);
419 ierr = MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
420 ierr = MatSetType(*mass, dmCoarse->mattype);CHKERRQ(ierr);
421 ierr = DMGetApplicationContext(dmFine, &ctx);CHKERRQ(ierr);
422
423 ierr = DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);CHKERRQ(ierr);
424 ierr = MatViewFromOptions(*mass, NULL, "-mass_mat_view");CHKERRQ(ierr);
425 PetscFunctionReturn(0);
426 }
427
DMSwarmComputeMassMatrixSquare_Private(DM dmc,DM dmf,Mat mass,PetscBool useDeltaFunction,void * ctx)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;
441 PetscErrorCode ierr;
442
443 PetscFunctionBegin;
444 ierr = PetscObjectGetComm((PetscObject) mass, &comm);CHKERRQ(ierr);
445 ierr = DMGetCoordinateDim(dmf, &cdim);CHKERRQ(ierr);
446 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr);
447 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr);
448 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
449 ierr = PetscMalloc3(cdim, &v0, cdim*cdim, &J, cdim*cdim,&invJ);CHKERRQ(ierr);
450 ierr = DMGetLocalSection(dmf, &fsection);CHKERRQ(ierr);
451 ierr = DMGetGlobalSection(dmf, &globalFSection);CHKERRQ(ierr);
452 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr);
453 ierr = MatGetLocalSize(mass, &locRows, &locCols);CHKERRQ(ierr);
454
455 ierr = PetscLayoutCreate(comm, &colLayout);CHKERRQ(ierr);
456 ierr = PetscLayoutSetLocalSize(colLayout, locCols);CHKERRQ(ierr);
457 ierr = PetscLayoutSetBlockSize(colLayout, 1);CHKERRQ(ierr);
458 ierr = PetscLayoutSetUp(colLayout);CHKERRQ(ierr);
459 ierr = PetscLayoutGetRange(colLayout, &colStart, &colEnd);CHKERRQ(ierr);
460 ierr = PetscLayoutDestroy(&colLayout);CHKERRQ(ierr);
461
462 ierr = PetscLayoutCreate(comm, &rLayout);CHKERRQ(ierr);
463 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
464 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
465 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
466 ierr = PetscLayoutGetRange(rLayout, &rStart, NULL);CHKERRQ(ierr);
467 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
468
469 ierr = DMPlexGetDepth(dmf, &depth);CHKERRQ(ierr);
470 ierr = DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
471 maxAdjSize = PetscPowInt(maxConeSize*maxSupportSize, depth);
472 ierr = PetscMalloc1(maxAdjSize, &adj);CHKERRQ(ierr);
473
474 ierr = PetscCalloc2(locRows, &dnz, locRows, &onz);CHKERRQ(ierr);
475 ierr = PetscHSetIJCreate(&ht);CHKERRQ(ierr);
476 /* Count nonzeros
477 This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
478 */
479 ierr = DMSwarmSortGetAccess(dmc);CHKERRQ(ierr);
480 for (cell = cStart; cell < cEnd; ++cell) {
481 PetscInt i;
482 PetscInt *cindices;
483 PetscInt numCIndices;
484 #if 0
485 PetscInt adjSize = maxAdjSize, a, j;
486 #endif
487
488 ierr = DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);CHKERRQ(ierr);
489 maxC = PetscMax(maxC, numCIndices);
490 /* Diagonal block */
491 for (i = 0; i < numCIndices; ++i) {dnz[cindices[i]] += numCIndices;}
492 #if 0
493 /* Off-diagonal blocks */
494 ierr = DMPlexGetAdjacency(dmf, cell, &adjSize, &adj);CHKERRQ(ierr);
495 for (a = 0; a < adjSize; ++a) {
496 if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
497 const PetscInt ncell = adj[a];
498 PetscInt *ncindices;
499 PetscInt numNCIndices;
500
501 ierr = DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices);CHKERRQ(ierr);
502 {
503 PetscHashIJKey key;
504 PetscBool missing;
505
506 for (i = 0; i < numCIndices; ++i) {
507 key.i = cindices[i] + rStart; /* global rows (from Swarm) */
508 if (key.i < 0) continue;
509 for (j = 0; j < numNCIndices; ++j) {
510 key.j = ncindices[j] + rStart; /* global column (from Swarm) */
511 if (key.j < 0) continue;
512 ierr = PetscHSetIJQueryAdd(ht, key, &missing);CHKERRQ(ierr);
513 if (missing) {
514 if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
515 else ++onz[key.i - rStart];
516 }
517 }
518 }
519 }
520 ierr = PetscFree(ncindices);CHKERRQ(ierr);
521 }
522 }
523 #endif
524 ierr = PetscFree(cindices);CHKERRQ(ierr);
525 }
526 ierr = PetscHSetIJDestroy(&ht);CHKERRQ(ierr);
527 ierr = MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr);
528 ierr = MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);CHKERRQ(ierr);
529 ierr = PetscFree2(dnz, onz);CHKERRQ(ierr);
530 ierr = PetscMalloc4(maxC*totDim, &elemMat, maxC*maxC, &elemMatSq, maxC, &rowIDXs, maxC*cdim, &xi);CHKERRQ(ierr);
531 /* Fill in values
532 Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
533 Start just by producing block diagonal
534 Could loop over adjacent cells
535 Produce neighboring element matrix
536 TODO Determine which columns and rows correspond to shared dual vector
537 Do MatMatMult with rectangular matrices
538 Insert block
539 */
540 for (field = 0; field < Nf; ++field) {
541 PetscTabulation Tcoarse;
542 PetscObject obj;
543 PetscReal *coords;
544 PetscInt Nc, i;
545
546 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr);
547 ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);
548 if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
549 ierr = DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr);
550 for (cell = cStart; cell < cEnd; ++cell) {
551 PetscInt *findices , *cindices;
552 PetscInt numFIndices, numCIndices;
553 PetscInt p, c;
554
555 /* TODO: Use DMField instead of assuming affine */
556 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
557 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr);
558 ierr = DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);CHKERRQ(ierr);
559 for (p = 0; p < numCIndices; ++p) {
560 CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p]*cdim], &xi[p*cdim]);
561 }
562 ierr = PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);CHKERRQ(ierr);
563 /* Get elemMat entries by multiplying by weight */
564 ierr = PetscArrayzero(elemMat, numCIndices*totDim);CHKERRQ(ierr);
565 for (i = 0; i < numFIndices; ++i) {
566 for (p = 0; p < numCIndices; ++p) {
567 for (c = 0; c < Nc; ++c) {
568 /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
569 elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
570 }
571 }
572 }
573 ierr = PetscTabulationDestroy(&Tcoarse);CHKERRQ(ierr);
574 for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
575 if (0) {ierr = DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);CHKERRQ(ierr);}
576 /* Block diagonal */
577 {
578 PetscBLASInt blasn, blask;
579 PetscScalar one = 1.0, zero = 0.0;
580
581 ierr = PetscBLASIntCast(numCIndices, &blasn);CHKERRQ(ierr);
582 ierr = PetscBLASIntCast(numFIndices, &blask);CHKERRQ(ierr);
583 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&blasn,&blasn,&blask,&one,elemMat,&blask,elemMat,&blask,&zero,elemMatSq,&blasn));
584 }
585 ierr = MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES);CHKERRQ(ierr);
586 /* TODO Off-diagonal */
587 ierr = PetscFree(cindices);CHKERRQ(ierr);
588 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);CHKERRQ(ierr);
589 }
590 ierr = DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr);
591 }
592 ierr = PetscFree4(elemMat, elemMatSq, rowIDXs, xi);CHKERRQ(ierr);
593 ierr = PetscFree(adj);CHKERRQ(ierr);
594 ierr = DMSwarmSortRestoreAccess(dmc);CHKERRQ(ierr);
595 ierr = PetscFree3(v0, J, invJ);CHKERRQ(ierr);
596 ierr = MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
597 ierr = MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
598 PetscFunctionReturn(0);
599 }
600
601 /*@C
602 DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
603
604 Collective on dmCoarse
605
606 Input parameters:
607 + dmCoarse - a DMSwarm
608 - dmFine - a DMPlex
609
610 Output parameter:
611 . mass - the square of the particle mass matrix
612
613 Level: advanced
614
615 Notes:
616 We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
617 future to compute the full normal equations.
618
619 .seealso: DMCreateMassMatrix()
620 @*/
DMSwarmCreateMassMatrixSquare(DM dmCoarse,DM dmFine,Mat * mass)621 PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
622 {
623 PetscInt n;
624 void *ctx;
625 PetscErrorCode ierr;
626
627 PetscFunctionBegin;
628 ierr = DMSwarmGetLocalSize(dmCoarse, &n);CHKERRQ(ierr);
629 ierr = MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);CHKERRQ(ierr);
630 ierr = MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
631 ierr = MatSetType(*mass, dmCoarse->mattype);CHKERRQ(ierr);
632 ierr = DMGetApplicationContext(dmFine, &ctx);CHKERRQ(ierr);
633
634 ierr = DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);CHKERRQ(ierr);
635 ierr = MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view");CHKERRQ(ierr);
636 PetscFunctionReturn(0);
637 }
638
639
640 /*@C
641 DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
642
643 Collective on dm
644
645 Input parameters:
646 + dm - a DMSwarm
647 - fieldname - the textual name given to a registered field
648
649 Output parameter:
650 . vec - the vector
651
652 Level: beginner
653
654 Notes:
655 The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
656
657 .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
658 @*/
DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec * vec)659 PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
660 {
661 MPI_Comm comm = PetscObjectComm((PetscObject) dm);
662 PetscErrorCode ierr;
663
664 PetscFunctionBegin;
665 ierr = DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);CHKERRQ(ierr);
666 PetscFunctionReturn(0);
667 }
668
669 /*@C
670 DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
671
672 Collective on dm
673
674 Input parameters:
675 + dm - a DMSwarm
676 - fieldname - the textual name given to a registered field
677
678 Output parameter:
679 . vec - the vector
680
681 Level: beginner
682
683 .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
684 @*/
DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec * vec)685 PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
686 {
687 PetscErrorCode ierr;
688
689 PetscFunctionBegin;
690 ierr = DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);CHKERRQ(ierr);
691 PetscFunctionReturn(0);
692 }
693
694 /*@C
695 DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
696
697 Collective on dm
698
699 Input parameters:
700 + dm - a DMSwarm
701 - fieldname - the textual name given to a registered field
702
703 Output parameter:
704 . vec - the vector
705
706 Level: beginner
707
708 Notes:
709 The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
710
711 .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
712 @*/
DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec * vec)713 PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
714 {
715 MPI_Comm comm = PETSC_COMM_SELF;
716 PetscErrorCode ierr;
717
718 PetscFunctionBegin;
719 ierr = DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);CHKERRQ(ierr);
720 PetscFunctionReturn(0);
721 }
722
723 /*@C
724 DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
725
726 Collective on dm
727
728 Input parameters:
729 + dm - a DMSwarm
730 - fieldname - the textual name given to a registered field
731
732 Output parameter:
733 . vec - the vector
734
735 Level: beginner
736
737 .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
738 @*/
DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec * vec)739 PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
740 {
741 PetscErrorCode ierr;
742
743 PetscFunctionBegin;
744 ierr = DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);CHKERRQ(ierr);
745 PetscFunctionReturn(0);
746 }
747
748 /*
749 PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
750 {
751 PetscFunctionReturn(0);
752 }
753
754 PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
755 {
756 PetscFunctionReturn(0);
757 }
758 */
759
760 /*@C
761 DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
762
763 Collective on dm
764
765 Input parameter:
766 . dm - a DMSwarm
767
768 Level: beginner
769
770 Notes:
771 After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
772
773 .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
774 DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
775 @*/
DMSwarmInitializeFieldRegister(DM dm)776 PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
777 {
778 DM_Swarm *swarm = (DM_Swarm *) dm->data;
779 PetscErrorCode ierr;
780
781 PetscFunctionBegin;
782 if (!swarm->field_registration_initialized) {
783 swarm->field_registration_initialized = PETSC_TRUE;
784 ierr = DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_INT64);CHKERRQ(ierr); /* unique identifer */
785 ierr = DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT);CHKERRQ(ierr); /* used for communication */
786 }
787 PetscFunctionReturn(0);
788 }
789
790 /*@
791 DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
792
793 Collective on dm
794
795 Input parameter:
796 . dm - a DMSwarm
797
798 Level: beginner
799
800 Notes:
801 After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
802
803 .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
804 DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
805 @*/
DMSwarmFinalizeFieldRegister(DM dm)806 PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
807 {
808 DM_Swarm *swarm = (DM_Swarm*)dm->data;
809 PetscErrorCode ierr;
810
811 PetscFunctionBegin;
812 if (!swarm->field_registration_finalized) {
813 ierr = DMSwarmDataBucketFinalize(swarm->db);CHKERRQ(ierr);
814 }
815 swarm->field_registration_finalized = PETSC_TRUE;
816 PetscFunctionReturn(0);
817 }
818
819 /*@
820 DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
821
822 Not collective
823
824 Input parameters:
825 + dm - a DMSwarm
826 . nlocal - the length of each registered field
827 - buffer - the length of the buffer used to efficient dynamic re-sizing
828
829 Level: beginner
830
831 .seealso: DMSwarmGetLocalSize()
832 @*/
DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)833 PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
834 {
835 DM_Swarm *swarm = (DM_Swarm*)dm->data;
836 PetscErrorCode ierr;
837
838 PetscFunctionBegin;
839 ierr = PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);CHKERRQ(ierr);
840 ierr = DMSwarmDataBucketSetSizes(swarm->db,nlocal,buffer);CHKERRQ(ierr);
841 ierr = PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);CHKERRQ(ierr);
842 PetscFunctionReturn(0);
843 }
844
845 /*@
846 DMSwarmSetCellDM - Attachs a DM to a DMSwarm
847
848 Collective on dm
849
850 Input parameters:
851 + dm - a DMSwarm
852 - dmcell - the DM to attach to the DMSwarm
853
854 Level: beginner
855
856 Notes:
857 The attached DM (dmcell) will be queried for point location and
858 neighbor MPI-rank information if DMSwarmMigrate() is called.
859
860 .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
861 @*/
DMSwarmSetCellDM(DM dm,DM dmcell)862 PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
863 {
864 DM_Swarm *swarm = (DM_Swarm*)dm->data;
865
866 PetscFunctionBegin;
867 swarm->dmcell = dmcell;
868 PetscFunctionReturn(0);
869 }
870
871 /*@
872 DMSwarmGetCellDM - Fetches the attached cell DM
873
874 Collective on dm
875
876 Input parameter:
877 . dm - a DMSwarm
878
879 Output parameter:
880 . dmcell - the DM which was attached to the DMSwarm
881
882 Level: beginner
883
884 .seealso: DMSwarmSetCellDM()
885 @*/
DMSwarmGetCellDM(DM dm,DM * dmcell)886 PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
887 {
888 DM_Swarm *swarm = (DM_Swarm*)dm->data;
889
890 PetscFunctionBegin;
891 *dmcell = swarm->dmcell;
892 PetscFunctionReturn(0);
893 }
894
895 /*@
896 DMSwarmGetLocalSize - Retrives the local length of fields registered
897
898 Not collective
899
900 Input parameter:
901 . dm - a DMSwarm
902
903 Output parameter:
904 . nlocal - the length of each registered field
905
906 Level: beginner
907
908 .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
909 @*/
DMSwarmGetLocalSize(DM dm,PetscInt * nlocal)910 PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
911 {
912 DM_Swarm *swarm = (DM_Swarm*)dm->data;
913 PetscErrorCode ierr;
914
915 PetscFunctionBegin;
916 if (nlocal) {ierr = DMSwarmDataBucketGetSizes(swarm->db,nlocal,NULL,NULL);CHKERRQ(ierr);}
917 PetscFunctionReturn(0);
918 }
919
920 /*@
921 DMSwarmGetSize - Retrives the total length of fields registered
922
923 Collective on dm
924
925 Input parameter:
926 . dm - a DMSwarm
927
928 Output parameter:
929 . n - the total length of each registered field
930
931 Level: beginner
932
933 Note:
934 This calls MPI_Allreduce upon each call (inefficient but safe)
935
936 .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
937 @*/
DMSwarmGetSize(DM dm,PetscInt * n)938 PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
939 {
940 DM_Swarm *swarm = (DM_Swarm*)dm->data;
941 PetscErrorCode ierr;
942 PetscInt nlocal,ng;
943
944 PetscFunctionBegin;
945 ierr = DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);CHKERRQ(ierr);
946 ierr = MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
947 if (n) { *n = ng; }
948 PetscFunctionReturn(0);
949 }
950
951 /*@C
952 DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
953
954 Collective on dm
955
956 Input parameters:
957 + dm - a DMSwarm
958 . fieldname - the textual name to identify this field
959 . blocksize - the number of each data type
960 - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
961
962 Level: beginner
963
964 Notes:
965 The textual name for each registered field must be unique.
966
967 .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
968 @*/
DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)969 PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
970 {
971 PetscErrorCode ierr;
972 DM_Swarm *swarm = (DM_Swarm*)dm->data;
973 size_t size;
974
975 PetscFunctionBegin;
976 if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
977 if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
978
979 if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
980 if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
981 if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
982 if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
983 if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
984
985 ierr = PetscDataTypeGetSize(type, &size);CHKERRQ(ierr);
986 /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
987 ierr = DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);CHKERRQ(ierr);
988 {
989 DMSwarmDataField gfield;
990
991 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
992 ierr = DMSwarmDataFieldSetBlockSize(gfield,blocksize);CHKERRQ(ierr);
993 }
994 swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
995 PetscFunctionReturn(0);
996 }
997
998 /*@C
999 DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
1000
1001 Collective on dm
1002
1003 Input parameters:
1004 + dm - a DMSwarm
1005 . fieldname - the textual name to identify this field
1006 - size - the size in bytes of the user struct of each data type
1007
1008 Level: beginner
1009
1010 Notes:
1011 The textual name for each registered field must be unique.
1012
1013 .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
1014 @*/
DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)1015 PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
1016 {
1017 PetscErrorCode ierr;
1018 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1019
1020 PetscFunctionBegin;
1021 ierr = DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);CHKERRQ(ierr);
1022 swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
1023 PetscFunctionReturn(0);
1024 }
1025
1026 /*@C
1027 DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
1028
1029 Collective on dm
1030
1031 Input parameters:
1032 + dm - a DMSwarm
1033 . fieldname - the textual name to identify this field
1034 . size - the size in bytes of the user data type
1035 - blocksize - the number of each data type
1036
1037 Level: beginner
1038
1039 Notes:
1040 The textual name for each registered field must be unique.
1041
1042 .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
1043 @*/
DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)1044 PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
1045 {
1046 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1047 PetscErrorCode ierr;
1048
1049 PetscFunctionBegin;
1050 ierr = DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);CHKERRQ(ierr);
1051 {
1052 DMSwarmDataField gfield;
1053
1054 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
1055 ierr = DMSwarmDataFieldSetBlockSize(gfield,blocksize);CHKERRQ(ierr);
1056 }
1057 swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1058 PetscFunctionReturn(0);
1059 }
1060
1061 /*@C
1062 DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1063
1064 Not collective
1065
1066 Input parameters:
1067 + dm - a DMSwarm
1068 - fieldname - the textual name to identify this field
1069
1070 Output parameters:
1071 + blocksize - the number of each data type
1072 . type - the data type
1073 - data - pointer to raw array
1074
1075 Level: beginner
1076
1077 Notes:
1078 The array must be returned using a matching call to DMSwarmRestoreField().
1079
1080 .seealso: DMSwarmRestoreField()
1081 @*/
DMSwarmGetField(DM dm,const char fieldname[],PetscInt * blocksize,PetscDataType * type,void ** data)1082 PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
1083 {
1084 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1085 DMSwarmDataField gfield;
1086 PetscErrorCode ierr;
1087
1088 PetscFunctionBegin;
1089 if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
1090 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
1091 ierr = DMSwarmDataFieldGetAccess(gfield);CHKERRQ(ierr);
1092 ierr = DMSwarmDataFieldGetEntries(gfield,data);CHKERRQ(ierr);
1093 if (blocksize) {*blocksize = gfield->bs; }
1094 if (type) { *type = gfield->petsc_type; }
1095 PetscFunctionReturn(0);
1096 }
1097
1098 /*@C
1099 DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1100
1101 Not collective
1102
1103 Input parameters:
1104 + dm - a DMSwarm
1105 - fieldname - the textual name to identify this field
1106
1107 Output parameters:
1108 + blocksize - the number of each data type
1109 . type - the data type
1110 - data - pointer to raw array
1111
1112 Level: beginner
1113
1114 Notes:
1115 The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
1116
1117 .seealso: DMSwarmGetField()
1118 @*/
DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt * blocksize,PetscDataType * type,void ** data)1119 PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
1120 {
1121 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1122 DMSwarmDataField gfield;
1123 PetscErrorCode ierr;
1124
1125 PetscFunctionBegin;
1126 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
1127 ierr = DMSwarmDataFieldRestoreAccess(gfield);CHKERRQ(ierr);
1128 if (data) *data = NULL;
1129 PetscFunctionReturn(0);
1130 }
1131
1132 /*@
1133 DMSwarmAddPoint - Add space for one new point in the DMSwarm
1134
1135 Not collective
1136
1137 Input parameter:
1138 . dm - a DMSwarm
1139
1140 Level: beginner
1141
1142 Notes:
1143 The new point will have all fields initialized to zero.
1144
1145 .seealso: DMSwarmAddNPoints()
1146 @*/
DMSwarmAddPoint(DM dm)1147 PetscErrorCode DMSwarmAddPoint(DM dm)
1148 {
1149 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1150 PetscErrorCode ierr;
1151
1152 PetscFunctionBegin;
1153 if (!swarm->issetup) {ierr = DMSetUp(dm);CHKERRQ(ierr);}
1154 ierr = PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);CHKERRQ(ierr);
1155 ierr = DMSwarmDataBucketAddPoint(swarm->db);CHKERRQ(ierr);
1156 ierr = PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);CHKERRQ(ierr);
1157 PetscFunctionReturn(0);
1158 }
1159
1160 /*@
1161 DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
1162
1163 Not collective
1164
1165 Input parameters:
1166 + dm - a DMSwarm
1167 - npoints - the number of new points to add
1168
1169 Level: beginner
1170
1171 Notes:
1172 The new point will have all fields initialized to zero.
1173
1174 .seealso: DMSwarmAddPoint()
1175 @*/
DMSwarmAddNPoints(DM dm,PetscInt npoints)1176 PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
1177 {
1178 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1179 PetscErrorCode ierr;
1180 PetscInt nlocal;
1181
1182 PetscFunctionBegin;
1183 ierr = PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);CHKERRQ(ierr);
1184 ierr = DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);CHKERRQ(ierr);
1185 nlocal = nlocal + npoints;
1186 ierr = DMSwarmDataBucketSetSizes(swarm->db,nlocal,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr);
1187 ierr = PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);CHKERRQ(ierr);
1188 PetscFunctionReturn(0);
1189 }
1190
1191 /*@
1192 DMSwarmRemovePoint - Remove the last point from the DMSwarm
1193
1194 Not collective
1195
1196 Input parameter:
1197 . dm - a DMSwarm
1198
1199 Level: beginner
1200
1201 .seealso: DMSwarmRemovePointAtIndex()
1202 @*/
DMSwarmRemovePoint(DM dm)1203 PetscErrorCode DMSwarmRemovePoint(DM dm)
1204 {
1205 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1206 PetscErrorCode ierr;
1207
1208 PetscFunctionBegin;
1209 ierr = PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);CHKERRQ(ierr);
1210 ierr = DMSwarmDataBucketRemovePoint(swarm->db);CHKERRQ(ierr);
1211 ierr = PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);CHKERRQ(ierr);
1212 PetscFunctionReturn(0);
1213 }
1214
1215 /*@
1216 DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
1217
1218 Not collective
1219
1220 Input parameters:
1221 + dm - a DMSwarm
1222 - idx - index of point to remove
1223
1224 Level: beginner
1225
1226 .seealso: DMSwarmRemovePoint()
1227 @*/
DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)1228 PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
1229 {
1230 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1231 PetscErrorCode ierr;
1232
1233 PetscFunctionBegin;
1234 ierr = PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);CHKERRQ(ierr);
1235 ierr = DMSwarmDataBucketRemovePointAtIndex(swarm->db,idx);CHKERRQ(ierr);
1236 ierr = PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);CHKERRQ(ierr);
1237 PetscFunctionReturn(0);
1238 }
1239
1240 /*@
1241 DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
1242
1243 Not collective
1244
1245 Input parameters:
1246 + dm - a DMSwarm
1247 . pi - the index of the point to copy
1248 - pj - the point index where the copy should be located
1249
1250 Level: beginner
1251
1252 .seealso: DMSwarmRemovePoint()
1253 @*/
DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)1254 PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
1255 {
1256 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1257 PetscErrorCode ierr;
1258
1259 PetscFunctionBegin;
1260 if (!swarm->issetup) {ierr = DMSetUp(dm);CHKERRQ(ierr);}
1261 ierr = DMSwarmDataBucketCopyPoint(swarm->db,pi,swarm->db,pj);CHKERRQ(ierr);
1262 PetscFunctionReturn(0);
1263 }
1264
DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)1265 PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
1266 {
1267 PetscErrorCode ierr;
1268
1269 PetscFunctionBegin;
1270 ierr = DMSwarmMigrate_Push_Basic(dm,remove_sent_points);CHKERRQ(ierr);
1271 PetscFunctionReturn(0);
1272 }
1273
1274 /*@
1275 DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
1276
1277 Collective on dm
1278
1279 Input parameters:
1280 + dm - the DMSwarm
1281 - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1282
1283 Notes:
1284 The DM will be modified to accomodate received points.
1285 If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
1286 Different styles of migration are supported. See DMSwarmSetMigrateType().
1287
1288 Level: advanced
1289
1290 .seealso: DMSwarmSetMigrateType()
1291 @*/
DMSwarmMigrate(DM dm,PetscBool remove_sent_points)1292 PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
1293 {
1294 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1295 PetscErrorCode ierr;
1296
1297 PetscFunctionBegin;
1298 ierr = PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);CHKERRQ(ierr);
1299 switch (swarm->migrate_type) {
1300 case DMSWARM_MIGRATE_BASIC:
1301 ierr = DMSwarmMigrate_Basic(dm,remove_sent_points);CHKERRQ(ierr);
1302 break;
1303 case DMSWARM_MIGRATE_DMCELLNSCATTER:
1304 ierr = DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);CHKERRQ(ierr);
1305 break;
1306 case DMSWARM_MIGRATE_DMCELLEXACT:
1307 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1308 break;
1309 case DMSWARM_MIGRATE_USER:
1310 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
1311 break;
1312 default:
1313 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
1314 break;
1315 }
1316 ierr = PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);CHKERRQ(ierr);
1317 PetscFunctionReturn(0);
1318 }
1319
1320 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
1321
1322 /*
1323 DMSwarmCollectViewCreate
1324
1325 * Applies a collection method and gathers point neighbour points into dm
1326
1327 Notes:
1328 Users should call DMSwarmCollectViewDestroy() after
1329 they have finished computations associated with the collected points
1330 */
1331
1332 /*@
1333 DMSwarmCollectViewCreate - Applies a collection method and gathers points
1334 in neighbour MPI-ranks into the DMSwarm
1335
1336 Collective on dm
1337
1338 Input parameter:
1339 . dm - the DMSwarm
1340
1341 Notes:
1342 Users should call DMSwarmCollectViewDestroy() after
1343 they have finished computations associated with the collected points
1344 Different collect methods are supported. See DMSwarmSetCollectType().
1345
1346 Level: advanced
1347
1348 .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
1349 @*/
DMSwarmCollectViewCreate(DM dm)1350 PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1351 {
1352 PetscErrorCode ierr;
1353 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1354 PetscInt ng;
1355
1356 PetscFunctionBegin;
1357 if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
1358 ierr = DMSwarmGetLocalSize(dm,&ng);CHKERRQ(ierr);
1359 switch (swarm->collect_type) {
1360
1361 case DMSWARM_COLLECT_BASIC:
1362 ierr = DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);CHKERRQ(ierr);
1363 break;
1364 case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1365 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1366 break;
1367 case DMSWARM_COLLECT_GENERAL:
1368 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
1369 break;
1370 default:
1371 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
1372 break;
1373 }
1374 swarm->collect_view_active = PETSC_TRUE;
1375 swarm->collect_view_reset_nlocal = ng;
1376 PetscFunctionReturn(0);
1377 }
1378
1379 /*@
1380 DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
1381
1382 Collective on dm
1383
1384 Input parameters:
1385 . dm - the DMSwarm
1386
1387 Notes:
1388 Users should call DMSwarmCollectViewCreate() before this function is called.
1389
1390 Level: advanced
1391
1392 .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
1393 @*/
DMSwarmCollectViewDestroy(DM dm)1394 PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1395 {
1396 PetscErrorCode ierr;
1397 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1398
1399 PetscFunctionBegin;
1400 if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
1401 ierr = DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);CHKERRQ(ierr);
1402 swarm->collect_view_active = PETSC_FALSE;
1403 PetscFunctionReturn(0);
1404 }
1405
DMSwarmSetUpPIC(DM dm)1406 PetscErrorCode DMSwarmSetUpPIC(DM dm)
1407 {
1408 PetscInt dim;
1409 PetscErrorCode ierr;
1410
1411 PetscFunctionBegin;
1412 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1413 if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1414 if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1415 ierr = DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);CHKERRQ(ierr);
1416 ierr = DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);CHKERRQ(ierr);
1417 PetscFunctionReturn(0);
1418 }
1419
1420 /*@
1421 DMSwarmSetType - Set particular flavor of DMSwarm
1422
1423 Collective on dm
1424
1425 Input parameters:
1426 + dm - the DMSwarm
1427 - stype - the DMSwarm type (e.g. DMSWARM_PIC)
1428
1429 Level: advanced
1430
1431 .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
1432 @*/
DMSwarmSetType(DM dm,DMSwarmType stype)1433 PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
1434 {
1435 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1436 PetscErrorCode ierr;
1437
1438 PetscFunctionBegin;
1439 swarm->swarm_type = stype;
1440 if (swarm->swarm_type == DMSWARM_PIC) {
1441 ierr = DMSwarmSetUpPIC(dm);CHKERRQ(ierr);
1442 }
1443 PetscFunctionReturn(0);
1444 }
1445
DMSetup_Swarm(DM dm)1446 PetscErrorCode DMSetup_Swarm(DM dm)
1447 {
1448 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1449 PetscErrorCode ierr;
1450 PetscMPIInt rank;
1451 PetscInt p,npoints,*rankval;
1452
1453 PetscFunctionBegin;
1454 if (swarm->issetup) PetscFunctionReturn(0);
1455
1456 swarm->issetup = PETSC_TRUE;
1457
1458 if (swarm->swarm_type == DMSWARM_PIC) {
1459 /* check dmcell exists */
1460 if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");
1461
1462 if (swarm->dmcell->ops->locatepointssubdomain) {
1463 /* check methods exists for exact ownership identificiation */
1464 ierr = PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");CHKERRQ(ierr);
1465 swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1466 } else {
1467 /* check methods exist for point location AND rank neighbor identification */
1468 if (swarm->dmcell->ops->locatepoints) {
1469 ierr = PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n");CHKERRQ(ierr);
1470 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1471
1472 if (swarm->dmcell->ops->getneighbors) {
1473 ierr = PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n");CHKERRQ(ierr);
1474 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1475
1476 swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1477 }
1478 }
1479
1480 ierr = DMSwarmFinalizeFieldRegister(dm);CHKERRQ(ierr);
1481
1482 /* check some fields were registered */
1483 if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
1484
1485 /* check local sizes were set */
1486 if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
1487
1488 /* initialize values in pid and rank placeholders */
1489 /* TODO: [pid - use MPI_Scan] */
1490 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
1491 ierr = DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
1492 ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
1493 for (p=0; p<npoints; p++) {
1494 rankval[p] = (PetscInt)rank;
1495 }
1496 ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
1497 PetscFunctionReturn(0);
1498 }
1499
1500 extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1501
DMDestroy_Swarm(DM dm)1502 PetscErrorCode DMDestroy_Swarm(DM dm)
1503 {
1504 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1505 PetscErrorCode ierr;
1506
1507 PetscFunctionBegin;
1508 ierr = DMSwarmDataBucketDestroy(&swarm->db);CHKERRQ(ierr);
1509 if (swarm->sort_context) {
1510 ierr = DMSwarmSortDestroy(&swarm->sort_context);CHKERRQ(ierr);
1511 }
1512 ierr = PetscFree(swarm);CHKERRQ(ierr);
1513 PetscFunctionReturn(0);
1514 }
1515
DMSwarmView_Draw(DM dm,PetscViewer viewer)1516 PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1517 {
1518 DM cdm;
1519 PetscDraw draw;
1520 PetscReal *coords, oldPause;
1521 PetscInt Np, p, bs;
1522 PetscErrorCode ierr;
1523
1524 PetscFunctionBegin;
1525 ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);CHKERRQ(ierr);
1526 ierr = DMSwarmGetCellDM(dm, &cdm);CHKERRQ(ierr);
1527 ierr = PetscDrawGetPause(draw, &oldPause);CHKERRQ(ierr);
1528 ierr = PetscDrawSetPause(draw, 0.0);CHKERRQ(ierr);
1529 ierr = DMView(cdm, viewer);CHKERRQ(ierr);
1530 ierr = PetscDrawSetPause(draw, oldPause);CHKERRQ(ierr);
1531
1532 ierr = DMSwarmGetLocalSize(dm, &Np);CHKERRQ(ierr);
1533 ierr = DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);CHKERRQ(ierr);
1534 for (p = 0; p < Np; ++p) {
1535 const PetscInt i = p*bs;
1536
1537 ierr = PetscDrawEllipse(draw, coords[i], coords[i+1], 0.01, 0.01, PETSC_DRAW_BLUE);CHKERRQ(ierr);
1538 }
1539 ierr = DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);CHKERRQ(ierr);
1540 ierr = PetscDrawFlush(draw);CHKERRQ(ierr);
1541 ierr = PetscDrawPause(draw);CHKERRQ(ierr);
1542 ierr = PetscDrawSave(draw);CHKERRQ(ierr);
1543 PetscFunctionReturn(0);
1544 }
1545
DMView_Swarm(DM dm,PetscViewer viewer)1546 PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1547 {
1548 DM_Swarm *swarm = (DM_Swarm*)dm->data;
1549 PetscBool iascii,ibinary,ishdf5,isvtk,isdraw;
1550 PetscErrorCode ierr;
1551
1552 PetscFunctionBegin;
1553 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1554 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
1555 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
1556 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);CHKERRQ(ierr);
1557 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
1558 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr);
1559 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);CHKERRQ(ierr);
1560 if (iascii) {
1561 ierr = DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);CHKERRQ(ierr);
1562 } else if (ibinary) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1563 #if defined(PETSC_HAVE_HDF5)
1564 else if (ishdf5) {
1565 ierr = DMSwarmView_HDF5(dm, viewer);CHKERRQ(ierr);
1566 }
1567 #else
1568 else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1569 #endif
1570 else if (isvtk) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1571 else if (isdraw) {
1572 ierr = DMSwarmView_Draw(dm, viewer);CHKERRQ(ierr);
1573 }
1574 PetscFunctionReturn(0);
1575 }
1576
1577 /*MC
1578
1579 DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1580 This implementation was designed for particle methods in which the underlying
1581 data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1582
1583 User data can be represented by DMSwarm through a registering "fields".
1584 To register a field, the user must provide;
1585 (a) a unique name;
1586 (b) the data type (or size in bytes);
1587 (c) the block size of the data.
1588
1589 For example, suppose the application requires a unique id, energy, momentum and density to be stored
1590 on a set of particles. Then the following code could be used
1591
1592 $ DMSwarmInitializeFieldRegister(dm)
1593 $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1594 $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1595 $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1596 $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1597 $ DMSwarmFinalizeFieldRegister(dm)
1598
1599 The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1600 The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1601
1602 To support particle methods, "migration" techniques are provided. These methods migrate data
1603 between MPI-ranks.
1604
1605 DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1606 As a DMSwarm may internally define and store values of different data types,
1607 before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1608 fields should be used to define a Vec object via
1609 DMSwarmVectorDefineField()
1610 The specified field can be changed at any time - thereby permitting vectors
1611 compatible with different fields to be created.
1612
1613 A dual representation of fields in the DMSwarm and a Vec object is permitted via
1614 DMSwarmCreateGlobalVectorFromField()
1615 Here the data defining the field in the DMSwarm is shared with a Vec.
1616 This is inherently unsafe if you alter the size of the field at any time between
1617 calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1618 If the local size of the DMSwarm does not match the local size of the global vector
1619 when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1620
1621 Additional high-level support is provided for Particle-In-Cell methods.
1622 Please refer to the man page for DMSwarmSetType().
1623
1624 Level: beginner
1625
1626 .seealso: DMType, DMCreate(), DMSetType()
1627 M*/
DMCreate_Swarm(DM dm)1628 PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1629 {
1630 DM_Swarm *swarm;
1631 PetscErrorCode ierr;
1632
1633 PetscFunctionBegin;
1634 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1635 ierr = PetscNewLog(dm,&swarm);CHKERRQ(ierr);
1636 dm->data = swarm;
1637
1638 ierr = DMSwarmDataBucketCreate(&swarm->db);CHKERRQ(ierr);
1639 ierr = DMSwarmInitializeFieldRegister(dm);CHKERRQ(ierr);
1640
1641 swarm->vec_field_set = PETSC_FALSE;
1642 swarm->issetup = PETSC_FALSE;
1643 swarm->swarm_type = DMSWARM_BASIC;
1644 swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1645 swarm->collect_type = DMSWARM_COLLECT_BASIC;
1646 swarm->migrate_error_on_missing_point = PETSC_FALSE;
1647
1648 swarm->dmcell = NULL;
1649 swarm->collect_view_active = PETSC_FALSE;
1650 swarm->collect_view_reset_nlocal = -1;
1651
1652 dm->dim = 0;
1653 dm->ops->view = DMView_Swarm;
1654 dm->ops->load = NULL;
1655 dm->ops->setfromoptions = NULL;
1656 dm->ops->clone = NULL;
1657 dm->ops->setup = DMSetup_Swarm;
1658 dm->ops->createlocalsection = NULL;
1659 dm->ops->createdefaultconstraints = NULL;
1660 dm->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1661 dm->ops->createlocalvector = DMCreateLocalVector_Swarm;
1662 dm->ops->getlocaltoglobalmapping = NULL;
1663 dm->ops->createfieldis = NULL;
1664 dm->ops->createcoordinatedm = NULL;
1665 dm->ops->getcoloring = NULL;
1666 dm->ops->creatematrix = NULL;
1667 dm->ops->createinterpolation = NULL;
1668 dm->ops->createinjection = NULL;
1669 dm->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1670 dm->ops->refine = NULL;
1671 dm->ops->coarsen = NULL;
1672 dm->ops->refinehierarchy = NULL;
1673 dm->ops->coarsenhierarchy = NULL;
1674 dm->ops->globaltolocalbegin = NULL;
1675 dm->ops->globaltolocalend = NULL;
1676 dm->ops->localtoglobalbegin = NULL;
1677 dm->ops->localtoglobalend = NULL;
1678 dm->ops->destroy = DMDestroy_Swarm;
1679 dm->ops->createsubdm = NULL;
1680 dm->ops->getdimpoints = NULL;
1681 dm->ops->locatepoints = NULL;
1682 PetscFunctionReturn(0);
1683 }
1684