Actual source code: swarm_ex3.c
petsc-3.9.4 2018-09-11
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>
12: PetscErrorCode _DMLocatePoints_DMDARegular_IS(DM dm,Vec pos,IS *iscell)
13: {
14: PetscInt p,n,bs,npoints,si,sj,milocal,mjlocal,mx,my;
15: DM dmregular;
16: PetscInt *cellidx;
17: const PetscScalar *coor;
18: PetscReal dx,dy;
20: 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,(void**)&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;
71: _DMLocatePoints_DMDARegular_IS(dm,pos,&iscell);
72: VecGetLocalSize(pos,&npoints);
73: VecGetBlockSize(pos,&bs);
74: npoints = npoints / bs;
76: PetscMalloc1(npoints, &cells);
77: ISGetIndices(iscell, &boxCells);
79: for (p=0; p<npoints; p++) {
80: cells[p].rank = 0;
81: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
82: cells[p].index = boxCells[p];
83: }
84: ISRestoreIndices(iscell, &boxCells);
85: ISDestroy(&iscell);
86: nfound = npoints;
87: PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
88: ISDestroy(&iscell);
89: return(0);
90: }
92: PetscErrorCode DMGetNeighbors_DMDARegular(DM dm,PetscInt *nneighbors,const PetscMPIInt **neighbors)
93: {
94: DM dmregular;
97: DMGetApplicationContext(dm,(void**)&dmregular);
98: DMGetNeighbors(dmregular,nneighbors,neighbors);
99: return(0);
100: }
102: PetscErrorCode SwarmViewGP(DM dms,const char prefix[])
103: {
104: PetscReal *array;
105: PetscInt *iarray;
106: PetscInt npoints,p,bs;
107: FILE *fp;
108: char name[PETSC_MAX_PATH_LEN];
109: PetscMPIInt rank;
112: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
113: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"%s-rank%d.gp",prefix,rank);
114: fp = fopen(name,"w");
115: if (!fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file %s",name);
116: DMSwarmGetLocalSize(dms,&npoints);
117: DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
118: DMSwarmGetField(dms,"itag",NULL,NULL,(void**)&iarray);
119: for (p=0; p<npoints; p++) {
120: fprintf(fp,"%+1.4e %+1.4e %1.4e\n",array[2*p],array[2*p+1],(double)iarray[p]);
121: }
122: DMSwarmRestoreField(dms,"itag",NULL,NULL,(void**)&iarray);
123: DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
124: fclose(fp);
125: return(0);
126: }
128: /*
129: Create a DMShell and attach a regularly spaced DMDA for point location
130: Override methods for point location
131: */
132: PetscErrorCode ex3_1(void)
133: {
134: DM dms,dmcell,dmregular;
135: PetscMPIInt rank;
136: PetscInt p,bs,nlocal,overlap,mx,tk;
137: PetscReal dx;
138: PetscReal *array,dt;
139: PetscInt *iarray;
140: PetscRandom rand;
144: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
146: /* Create a regularly spaced DMDA */
147: mx = 40;
148: overlap = 0;
149: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,mx,PETSC_DECIDE,PETSC_DECIDE,1,overlap,NULL,NULL,&dmregular);
150: DMSetFromOptions(dmregular);
151: DMSetUp(dmregular);
153: dx = 2.0/((PetscReal)mx);
154: 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);
156: /* Create a DMShell for point location purposes */
157: DMShellCreate(PETSC_COMM_WORLD,&dmcell);
158: DMSetApplicationContext(dmcell,(void*)dmregular);
159: dmcell->ops->locatepoints = DMLocatePoints_DMDARegular;
160: dmcell->ops->getneighbors = DMGetNeighbors_DMDARegular;
162: /* Create the swarm */
163: DMCreate(PETSC_COMM_WORLD,&dms);
164: DMSetType(dms,DMSWARM);
165: DMSetDimension(dms,2);
167: DMSwarmSetType(dms,DMSWARM_PIC);
168: DMSwarmSetCellDM(dms,dmcell);
170: DMSwarmRegisterPetscDatatypeField(dms,"itag",1,PETSC_INT);
171: DMSwarmFinalizeFieldRegister(dms);
172: {
173: PetscInt si,sj,milocal,mjlocal;
174: const PetscScalar *LA_coors;
175: Vec coors;
176: PetscInt cnt;
178: DMDAGetCorners(dmregular,&si,&sj,NULL,&milocal,&mjlocal,NULL);
179: DMGetCoordinates(dmregular,&coors);
180: VecGetArrayRead(coors,&LA_coors);
181: DMSwarmSetLocalSizes(dms,milocal*mjlocal,4);
182: DMSwarmGetLocalSize(dms,&nlocal);
183: DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
184: cnt = 0;
185: PetscRandomCreate(PETSC_COMM_SELF,&rand);
186: PetscRandomSetInterval(rand,-1.0,1.0);
187: for (p=0; p<nlocal; p++) {
188: PetscReal px,py,rx,ry,r2;
190: PetscRandomGetValueReal(rand,&rx);
191: PetscRandomGetValueReal(rand,&ry);
193: px = PetscRealPart(LA_coors[2*p+0]) + 0.1*rx*dx;
194: py = PetscRealPart(LA_coors[2*p+1]) + 0.1*ry*dx;
196: r2 = px*px + py*py;
197: if (r2 < 0.75*0.75) {
198: array[bs*cnt+0] = px;
199: array[bs*cnt+1] = py;
200: cnt++;
201: }
202: }
203: PetscRandomDestroy(&rand);
204: DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
205: VecRestoreArrayRead(coors,&LA_coors);
206: DMSwarmSetLocalSizes(dms,cnt,4);
208: DMSwarmGetLocalSize(dms,&nlocal);
209: DMSwarmGetField(dms,"itag",&bs,NULL,(void**)&iarray);
210: for (p=0; p<nlocal; p++) {
211: iarray[p] = (PetscInt)rank;
212: }
213: DMSwarmRestoreField(dms,"itag",&bs,NULL,(void**)&iarray);
214: }
216: DMView(dms,PETSC_VIEWER_STDOUT_WORLD);
217: SwarmViewGP(dms,"step0");
219: dt = 0.1;
220: for (tk=1; tk<20; tk++) {
221: char prefix[PETSC_MAX_PATH_LEN];
222: PetscPrintf(PETSC_COMM_WORLD,"Step %D \n",tk);
223: /* push points */
224: DMSwarmGetLocalSize(dms,&nlocal);
225: DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
226: for (p=0; p<nlocal; p++) {
227: PetscReal cx,cy,vx,vy;
229: cx = array[2*p];
230: cy = array[2*p+1];
231: vx = cy;
232: vy = -cx;
234: array[2*p ] += dt * vx;
235: array[2*p+1] += dt * vy;
236: }
237: DMSwarmRestoreField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);
239: /* migrate points */
240: DMSwarmMigrate(dms,PETSC_TRUE);
241: /* view points */
242: PetscSNPrintf(prefix,PETSC_MAX_PATH_LEN-1,"step%d",tk);
243: /* should use the regular SwarmView() api, not one for a particular type */
244: SwarmViewGP(dms,prefix);
245: }
246: DMDestroy(&dmregular);
247: DMDestroy(&dmcell);
248: DMDestroy(&dms);
249: return(0);
250: }
252: int main(int argc,char **argv)
253: {
256: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
257: ex3_1();
258: PetscFinalize();
259: return ierr;
260: }
262: /*TEST
264: build:
265: requires: double !complex
267: test:
268: filter: grep -v atomic
269: filter_output: grep -v atomic
270: TEST*/