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