Actual source code: swarm_ex3.c
2: static char help[] = "Tests DMSwarm with DMShell\n\n";
4: #include <petscsf.h>
5: #include <petscdm.h>
6: #include <petscdmshell.h>
7: #include <petscdmda.h>
8: #include <petscdmswarm.h>
9: #include <petsc/private/dmimpl.h>
11: PetscErrorCode _DMLocatePoints_DMDARegular_IS(DM dm, Vec pos, IS *iscell)
12: {
13: PetscInt p, n, bs, npoints, si, sj, milocal, mjlocal, mx, my;
14: DM dmregular;
15: PetscInt *cellidx;
16: const PetscScalar *coor;
17: PetscReal dx, dy;
18: PetscMPIInt rank;
20: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
21: VecGetLocalSize(pos, &n);
22: VecGetBlockSize(pos, &bs);
23: npoints = n / bs;
25: PetscMalloc1(npoints, &cellidx);
26: DMGetApplicationContext(dm, &dmregular);
27: DMDAGetCorners(dmregular, &si, &sj, NULL, &milocal, &mjlocal, NULL);
28: DMDAGetInfo(dmregular, NULL, &mx, &my, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
30: dx = 2.0 / ((PetscReal)mx);
31: dy = 2.0 / ((PetscReal)my);
33: VecGetArrayRead(pos, &coor);
34: for (p = 0; p < npoints; p++) {
35: PetscReal coorx, coory;
36: PetscInt mi, mj;
38: coorx = PetscRealPart(coor[2 * p]);
39: coory = PetscRealPart(coor[2 * p + 1]);
41: mi = (PetscInt)((coorx - (-1.0)) / dx);
42: mj = (PetscInt)((coory - (-1.0)) / dy);
44: cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
46: if ((mj >= sj) && (mj < sj + mjlocal)) {
47: if ((mi >= si) && (mi < si + milocal)) cellidx[p] = (mi - si) + (mj - sj) * milocal;
48: }
49: if (coorx < -1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
50: if (coorx > 1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
51: if (coory < -1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
52: if (coory > 1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
53: }
54: VecRestoreArrayRead(pos, &coor);
55: ISCreateGeneral(PETSC_COMM_SELF, npoints, cellidx, PETSC_OWN_POINTER, iscell);
56: return 0;
57: }
59: PetscErrorCode DMLocatePoints_DMDARegular(DM dm, Vec pos, DMPointLocationType ltype, PetscSF cellSF)
60: {
61: IS iscell;
62: PetscSFNode *cells;
63: PetscInt p, bs, npoints, nfound;
64: const PetscInt *boxCells;
66: _DMLocatePoints_DMDARegular_IS(dm, pos, &iscell);
67: VecGetLocalSize(pos, &npoints);
68: VecGetBlockSize(pos, &bs);
69: npoints = npoints / bs;
71: PetscMalloc1(npoints, &cells);
72: ISGetIndices(iscell, &boxCells);
74: for (p = 0; p < npoints; p++) {
75: cells[p].rank = 0;
76: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
77: cells[p].index = boxCells[p];
78: }
79: ISRestoreIndices(iscell, &boxCells);
80: ISDestroy(&iscell);
81: nfound = npoints;
82: PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
83: ISDestroy(&iscell);
84: return 0;
85: }
87: PetscErrorCode DMGetNeighbors_DMDARegular(DM dm, PetscInt *nneighbors, const PetscMPIInt **neighbors)
88: {
89: DM dmregular;
91: DMGetApplicationContext(dm, &dmregular);
92: DMGetNeighbors(dmregular, nneighbors, neighbors);
93: return 0;
94: }
96: PetscErrorCode SwarmViewGP(DM dms, const char prefix[])
97: {
98: PetscReal *array;
99: PetscInt *iarray;
100: PetscInt npoints, p, bs;
101: FILE *fp;
102: char name[PETSC_MAX_PATH_LEN];
103: PetscMPIInt rank;
105: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
106: PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "%s-rank%d.gp", prefix, rank);
107: fp = fopen(name, "w");
109: DMSwarmGetLocalSize(dms, &npoints);
110: DMSwarmGetField(dms, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
111: DMSwarmGetField(dms, "itag", NULL, NULL, (void **)&iarray);
112: for (p = 0; p < npoints; p++) fprintf(fp, "%+1.4e %+1.4e %1.4e\n", array[2 * p], array[2 * p + 1], (double)iarray[p]);
113: DMSwarmRestoreField(dms, "itag", NULL, NULL, (void **)&iarray);
114: DMSwarmRestoreField(dms, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
115: fclose(fp);
116: return 0;
117: }
119: /*
120: Create a DMShell and attach a regularly spaced DMDA for point location
121: Override methods for point location
122: */
123: PetscErrorCode ex3_1(void)
124: {
125: DM dms, dmcell, dmregular;
126: PetscMPIInt rank;
127: PetscInt p, bs, nlocal, overlap, mx, tk;
128: PetscReal dx;
129: PetscReal *array, dt;
130: PetscInt *iarray;
131: PetscRandom rand;
133: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
135: /* Create a regularly spaced DMDA */
136: mx = 40;
137: overlap = 0;
138: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, mx, PETSC_DECIDE, PETSC_DECIDE, 1, overlap, NULL, NULL, &dmregular);
139: DMSetFromOptions(dmregular);
140: DMSetUp(dmregular);
142: dx = 2.0 / ((PetscReal)mx);
143: DMDASetUniformCoordinates(dmregular, -1.0 + 0.5 * dx, 1.0 - 0.5 * dx, -1.0 + 0.5 * dx, 1.0 - 0.5 * dx, -1.0, 1.0);
145: /* Create a DMShell for point location purposes */
146: DMShellCreate(PETSC_COMM_WORLD, &dmcell);
147: DMSetApplicationContext(dmcell, dmregular);
148: dmcell->ops->locatepoints = DMLocatePoints_DMDARegular;
149: dmcell->ops->getneighbors = DMGetNeighbors_DMDARegular;
151: /* Create the swarm */
152: DMCreate(PETSC_COMM_WORLD, &dms);
153: DMSetType(dms, DMSWARM);
154: DMSetDimension(dms, 2);
155: PetscObjectSetName((PetscObject)dms, "Particles");
157: DMSwarmSetType(dms, DMSWARM_PIC);
158: DMSwarmSetCellDM(dms, dmcell);
160: DMSwarmRegisterPetscDatatypeField(dms, "itag", 1, PETSC_INT);
161: DMSwarmFinalizeFieldRegister(dms);
162: {
163: PetscInt si, sj, milocal, mjlocal;
164: const PetscScalar *LA_coors;
165: Vec coors;
166: PetscInt cnt;
168: DMDAGetCorners(dmregular, &si, &sj, NULL, &milocal, &mjlocal, NULL);
169: DMGetCoordinates(dmregular, &coors);
170: VecGetArrayRead(coors, &LA_coors);
171: DMSwarmSetLocalSizes(dms, milocal * mjlocal, 4);
172: DMSwarmGetLocalSize(dms, &nlocal);
173: DMSwarmGetField(dms, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
174: cnt = 0;
175: PetscRandomCreate(PETSC_COMM_SELF, &rand);
176: PetscRandomSetInterval(rand, -1.0, 1.0);
177: for (p = 0; p < nlocal; p++) {
178: PetscReal px, py, rx, ry, r2;
180: PetscRandomGetValueReal(rand, &rx);
181: PetscRandomGetValueReal(rand, &ry);
183: px = PetscRealPart(LA_coors[2 * p + 0]) + 0.1 * rx * dx;
184: py = PetscRealPart(LA_coors[2 * p + 1]) + 0.1 * ry * dx;
186: r2 = px * px + py * py;
187: if (r2 < 0.75 * 0.75) {
188: array[bs * cnt + 0] = px;
189: array[bs * cnt + 1] = py;
190: cnt++;
191: }
192: }
193: PetscRandomDestroy(&rand);
194: DMSwarmRestoreField(dms, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
195: VecRestoreArrayRead(coors, &LA_coors);
196: DMSwarmSetLocalSizes(dms, cnt, 4);
198: DMSwarmGetLocalSize(dms, &nlocal);
199: DMSwarmGetField(dms, "itag", &bs, NULL, (void **)&iarray);
200: for (p = 0; p < nlocal; p++) iarray[p] = (PetscInt)rank;
201: DMSwarmRestoreField(dms, "itag", &bs, NULL, (void **)&iarray);
202: }
204: DMView(dms, PETSC_VIEWER_STDOUT_WORLD);
205: SwarmViewGP(dms, "step0");
207: dt = 0.1;
208: for (tk = 1; tk < 20; tk++) {
209: char prefix[PETSC_MAX_PATH_LEN];
210: PetscPrintf(PETSC_COMM_WORLD, "Step %" PetscInt_FMT " \n", tk);
211: /* push points */
212: DMSwarmGetLocalSize(dms, &nlocal);
213: DMSwarmGetField(dms, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
214: for (p = 0; p < nlocal; p++) {
215: PetscReal cx, cy, vx, vy;
217: cx = array[2 * p];
218: cy = array[2 * p + 1];
219: vx = cy;
220: vy = -cx;
222: array[2 * p] += dt * vx;
223: array[2 * p + 1] += dt * vy;
224: }
225: DMSwarmRestoreField(dms, DMSwarmPICField_coor, &bs, NULL, (void **)&array);
227: /* migrate points */
228: DMSwarmMigrate(dms, PETSC_TRUE);
229: /* view points */
230: PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN - 1, "step%" PetscInt_FMT, tk);
231: /* should use the regular SwarmView() api, not one for a particular type */
232: SwarmViewGP(dms, prefix);
233: }
234: DMDestroy(&dmregular);
235: DMDestroy(&dmcell);
236: DMDestroy(&dms);
237: return 0;
238: }
240: int main(int argc, char **argv)
241: {
243: PetscInitialize(&argc, &argv, (char *)0, help);
244: ex3_1();
245: PetscFinalize();
246: return 0;
247: }
249: /*TEST
251: build:
252: requires: double !complex
254: test:
255: filter: grep -v atomic
256: filter_output: grep -v atomic
257: TEST*/