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>


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