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)) {
48: cellidx[p] = (mi-si) + (mj-sj) * milocal;
49: }
50: }
51: if (coorx < -1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
52: if (coorx > 1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
53: if (coory < -1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
54: if (coory > 1.0) cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND;
55: }
56: VecRestoreArrayRead(pos,&coor);
57: ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);
58: return 0;
59: }
61: PetscErrorCode DMLocatePoints_DMDARegular(DM dm,Vec pos,DMPointLocationType ltype, PetscSF cellSF)
62: {
63: IS iscell;
64: PetscSFNode *cells;
65: PetscInt p,bs,npoints,nfound;
66: const PetscInt *boxCells;
68: _DMLocatePoints_DMDARegular_IS(dm,pos,&iscell);
69: VecGetLocalSize(pos,&npoints);
70: VecGetBlockSize(pos,&bs);
71: npoints = npoints / bs;
73: PetscMalloc1(npoints, &cells);
74: ISGetIndices(iscell, &boxCells);
76: for (p=0; p<npoints; p++) {
77: cells[p].rank = 0;
78: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
79: cells[p].index = boxCells[p];
80: }
81: ISRestoreIndices(iscell, &boxCells);
82: ISDestroy(&iscell);
83: nfound = npoints;
84: PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
85: ISDestroy(&iscell);
86: return 0;
87: }
89: PetscErrorCode DMGetNeighbors_DMDARegular(DM dm,PetscInt *nneighbors,const PetscMPIInt **neighbors)
90: {
91: DM dmregular;
93: DMGetApplicationContext(dm,&dmregular);
94: DMGetNeighbors(dmregular,nneighbors,neighbors);
95: return 0;
96: }
98: PetscErrorCode SwarmViewGP(DM dms,const char prefix[])
99: {
100: PetscReal *array;
101: PetscInt *iarray;
102: PetscInt npoints,p,bs;
103: FILE *fp;
104: char name[PETSC_MAX_PATH_LEN];
105: PetscMPIInt rank;
107: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
108: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"%s-rank%d.gp",prefix,rank);
109: fp = fopen(name,"w");
111: DMSwarmGetLocalSize(dms,&npoints);
112: DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
113: DMSwarmGetField(dms,"itag",NULL,NULL,(void**)&iarray);
114: for (p=0; p<npoints; p++) {
115: fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array[2*p],array[2*p+1],(double)iarray[p]);
116: }
117: DMSwarmRestoreField(dms,"itag",NULL,NULL,(void**)&iarray);
118: DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
119: fclose(fp);
120: return 0;
121: }
123: /*
124: Create a DMShell and attach a regularly spaced DMDA for point location
125: Override methods for point location
126: */
127: PetscErrorCode ex3_1(void)
128: {
129: DM dms,dmcell,dmregular;
130: PetscMPIInt rank;
131: PetscInt p,bs,nlocal,overlap,mx,tk;
132: PetscReal dx;
133: PetscReal *array,dt;
134: PetscInt *iarray;
135: PetscRandom rand;
137: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
139: /* Create a regularly spaced DMDA */
140: mx = 40;
141: overlap = 0;
142: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,mx,PETSC_DECIDE,PETSC_DECIDE,1,overlap,NULL,NULL,&dmregular);
143: DMSetFromOptions(dmregular);
144: DMSetUp(dmregular);
146: dx = 2.0/((PetscReal)mx);
147: 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);
149: /* Create a DMShell for point location purposes */
150: DMShellCreate(PETSC_COMM_WORLD,&dmcell);
151: DMSetApplicationContext(dmcell,dmregular);
152: dmcell->ops->locatepoints = DMLocatePoints_DMDARegular;
153: dmcell->ops->getneighbors = DMGetNeighbors_DMDARegular;
155: /* Create the swarm */
156: DMCreate(PETSC_COMM_WORLD,&dms);
157: DMSetType(dms,DMSWARM);
158: DMSetDimension(dms,2);
160: DMSwarmSetType(dms,DMSWARM_PIC);
161: DMSwarmSetCellDM(dms,dmcell);
163: DMSwarmRegisterPetscDatatypeField(dms,"itag",1,PETSC_INT);
164: DMSwarmFinalizeFieldRegister(dms);
165: {
166: PetscInt si,sj,milocal,mjlocal;
167: const PetscScalar *LA_coors;
168: Vec coors;
169: PetscInt cnt;
171: DMDAGetCorners(dmregular,&si,&sj,NULL,&milocal,&mjlocal,NULL);
172: DMGetCoordinates(dmregular,&coors);
173: VecGetArrayRead(coors,&LA_coors);
174: DMSwarmSetLocalSizes(dms,milocal*mjlocal,4);
175: DMSwarmGetLocalSize(dms,&nlocal);
176: DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
177: cnt = 0;
178: PetscRandomCreate(PETSC_COMM_SELF,&rand);
179: PetscRandomSetInterval(rand,-1.0,1.0);
180: for (p=0; p<nlocal; p++) {
181: PetscReal px,py,rx,ry,r2;
183: PetscRandomGetValueReal(rand,&rx);
184: PetscRandomGetValueReal(rand,&ry);
186: px = PetscRealPart(LA_coors[2*p+0]) + 0.1*rx*dx;
187: py = PetscRealPart(LA_coors[2*p+1]) + 0.1*ry*dx;
189: r2 = px*px + py*py;
190: if (r2 < 0.75*0.75) {
191: array[bs*cnt+0] = px;
192: array[bs*cnt+1] = py;
193: cnt++;
194: }
195: }
196: PetscRandomDestroy(&rand);
197: DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
198: VecRestoreArrayRead(coors,&LA_coors);
199: DMSwarmSetLocalSizes(dms,cnt,4);
201: DMSwarmGetLocalSize(dms,&nlocal);
202: DMSwarmGetField(dms,"itag",&bs,NULL,(void**)&iarray);
203: for (p=0; p<nlocal; p++) {
204: iarray[p] = (PetscInt)rank;
205: }
206: DMSwarmRestoreField(dms,"itag",&bs,NULL,(void**)&iarray);
207: }
209: DMView(dms,PETSC_VIEWER_STDOUT_WORLD);
210: SwarmViewGP(dms,"step0");
212: dt = 0.1;
213: for (tk=1; tk<20; tk++) {
214: char prefix[PETSC_MAX_PATH_LEN];
215: PetscPrintf(PETSC_COMM_WORLD,"Step %D \n",tk);
216: /* push points */
217: DMSwarmGetLocalSize(dms,&nlocal);
218: DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
219: for (p=0; p<nlocal; p++) {
220: PetscReal cx,cy,vx,vy;
222: cx = array[2*p];
223: cy = array[2*p+1];
224: vx = cy;
225: vy = -cx;
227: array[2*p ] += dt * vx;
228: array[2*p+1] += dt * vy;
229: }
230: DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
232: /* migrate points */
233: DMSwarmMigrate(dms,PETSC_TRUE);
234: /* view points */
235: PetscSNPrintf(prefix,PETSC_MAX_PATH_LEN-1,"step%d",tk);
236: /* should use the regular SwarmView() api, not one for a particular type */
237: SwarmViewGP(dms,prefix);
238: }
239: DMDestroy(&dmregular);
240: DMDestroy(&dmcell);
241: DMDestroy(&dms);
242: return 0;
243: }
245: int main(int argc,char **argv)
246: {
248: PetscInitialize(&argc,&argv,(char*)0,help);
249: ex3_1();
250: PetscFinalize();
251: return 0;
252: }
254: /*TEST
256: build:
257: requires: double !complex
259: test:
260: filter: grep -v atomic
261: filter_output: grep -v atomic
262: TEST*/