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*/