Actual source code: dadd.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <petsc/private/dmdaimpl.h>  /*I   "petscdmda.h"   I*/

  5: /*@
  6:   DMDACreatePatchIS - Creates an index set corresponding to a patch of the DA.

  8:   Not Collective

 10:   Input Parameters:
 11: +  da - the DMDA
 12: .  lower - a matstencil with i, j and k corresponding to the lower corner of the patch
 13: -  upper - a matstencil with i, j and k corresponding to the upper corner of the patch

 15:   Output Parameters:
 16: .  is - the IS corresponding to the patch

 18:   Level: developer

 20: .seealso: DMDACreateDomainDecomposition(), DMDACreateDomainDecompositionScatters()
 21: @*/
 22: PetscErrorCode DMDACreatePatchIS(DM da,MatStencil *lower,MatStencil *upper,IS *is)
 23: {
 24:   PetscInt       ms=0,ns=0,ps=0;
 25:   PetscInt       me=1,ne=1,pe=1;
 26:   PetscInt       mr=0,nr=0,pr=0;
 27:   PetscInt       ii,jj,kk;
 28:   PetscInt       si,sj,sk;
 29:   PetscInt       i,j,k,l,idx;
 30:   PetscInt       base;
 31:   PetscInt       xm=1,ym=1,zm=1;
 32:   const PetscInt *lx,*ly,*lz;
 33:   PetscInt       ox,oy,oz;
 34:   PetscInt       m,n,p,M,N,P,dof;
 35:   PetscInt       nindices;
 36:   PetscInt       *indices;
 37:   DM_DA          *dd = (DM_DA*)da->data;

 41:   /* need to get the sizes of the actual DM rather than the "global" space of a subdomain DM */
 42:   M = dd->M;N = dd->N;P=dd->P;
 43:   m = dd->m;n = dd->n;p=dd->p;
 44:   dof = dd->w;
 45:   DMDAGetOffset(da,&ox,&oy,&oz,NULL,NULL,NULL);
 46:   DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
 47:   nindices = (upper->i - lower->i)*(upper->j - lower->j)*(upper->k - lower->k)*dof;
 48:   PetscMalloc1(nindices,&indices);
 49:   /* start at index 0 on processor 0 */
 50:   mr = 0;
 51:   nr = 0;
 52:   pr = 0;
 53:   ms = 0;
 54:   ns = 0;
 55:   ps = 0;
 56:   if (lx) me = lx[0];
 57:   if (ly) ne = ly[0];
 58:   if (lz) pe = lz[0];
 59:   idx = 0;
 60:   for (k=lower->k-oz;k<upper->k-oz;k++) {
 61:     for (j=lower->j-oy;j < upper->j-oy;j++) {
 62:       for (i=lower->i-ox;i < upper->i-ox;i++) {
 63:         /* "actual" indices rather than ones outside of the domain */
 64:         ii = i;
 65:         jj = j;
 66:         kk = k;
 67:         if (ii < 0) ii = M + ii;
 68:         if (jj < 0) jj = N + jj;
 69:         if (kk < 0) kk = P + kk;
 70:         if (ii > M-1) ii = ii - M;
 71:         if (jj > N-1) jj = jj - N;
 72:         if (kk > P-1) kk = kk - P;
 73:         /* gone out of processor range on x axis */
 74:         while(ii > me-1 || ii < ms) {
 75:           if (mr == m-1) {
 76:             ms = 0;
 77:             me = lx[0];
 78:             mr = 0;
 79:           } else {
 80:             mr++;
 81:             ms = me;
 82:             me += lx[mr];
 83:           }
 84:         }
 85:         /* gone out of processor range on y axis */
 86:         while(jj > ne-1 || jj < ns) {
 87:           if (nr == n-1) {
 88:             ns = 0;
 89:             ne = ly[0];
 90:             nr = 0;
 91:           } else {
 92:             nr++;
 93:             ns = ne;
 94:             ne += ly[nr];
 95:           }
 96:         }
 97:         /* gone out of processor range on z axis */
 98:         while(kk > pe-1 || kk < ps) {
 99:           if (pr == p-1) {
100:             ps = 0;
101:             pe = lz[0];
102:             pr = 0;
103:           } else {
104:             pr++;
105:             ps = pe;
106:             pe += lz[pr];
107:           }
108:         }
109:         /* compute the vector base on owning processor */
110:         xm = me - ms;
111:         ym = ne - ns;
112:         zm = pe - ps;
113:         base = ms*ym*zm + ns*M*zm + ps*M*N;
114:         /* compute the local coordinates on owning processor */
115:         si = ii - ms;
116:         sj = jj - ns;
117:         sk = kk - ps;
118:         for (l=0;l<dof;l++) {
119:           indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk);
120:           idx++;
121:         }
122:       }
123:     }
124:   }
125:   ISCreateGeneral(PETSC_COMM_SELF,idx,indices,PETSC_OWN_POINTER,is);
126:   return(0);
127: }

131: PetscErrorCode DMDASubDomainDA_Private(DM dm, PetscInt *nlocal, DM **sdm)
132: {
133:   DM             *da;
134:   PetscInt       dim,size,i,j,k,idx;
136:   DMDALocalInfo  info;
137:   PetscInt       xsize,ysize,zsize;
138:   PetscInt       xo,yo,zo;
139:   PetscInt       xs,ys,zs;
140:   PetscInt       xm=1,ym=1,zm=1;
141:   PetscInt       xol,yol,zol;
142:   PetscInt       m=1,n=1,p=1;
143:   PetscInt       M,N,P;
144:   PetscInt       pm,mtmp;

147:   DMDAGetLocalInfo(dm,&info);
148:   DMDAGetOverlap(dm,&xol,&yol,&zol);
149:   DMDAGetNumLocalSubDomains(dm,&size);
150:   PetscMalloc1(size,&da);

152:   if (nlocal) *nlocal = size;

154:   dim = info.dim;

156:   M = info.xm;
157:   N = info.ym;
158:   P = info.zm;

160:   if (dim == 1) {
161:     m = size;
162:   } else if (dim == 2) {
163:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N)));
164:     while (m > 0) {
165:       n = size/m;
166:       if (m*n*p == size) break;
167:       m--;
168:     }
169:   } else if (dim == 3) {
170:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.)));    if (!n) n = 1;
171:     while (n > 0) {
172:       pm = size/n;
173:       if (n*pm == size) break;
174:       n--;
175:     }
176:     if (!n) n = 1;
177:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
178:     if (!m) m = 1;
179:     while (m > 0) {
180:       p = size/(m*n);
181:       if (m*n*p == size) break;
182:       m--;
183:     }
184:     if (M > P && m < p) {mtmp = m; m = p; p = mtmp;}
185:   }

187:   zs = info.zs;
188:   idx = 0;
189:   for (k = 0; k < p; k++) {
190:     ys = info.ys;
191:     for (j = 0; j < n; j++) {
192:       xs = info.xs;
193:       for (i = 0; i < m; i++) {
194:         if (dim == 1) {
195:           xm = M/m + ((M % m) > i);
196:         } else if (dim == 2) {
197:           xm = M/m + ((M % m) > i);
198:           ym = N/n + ((N % n) > j);
199:         } else if (dim == 3) {
200:           xm = M/m + ((M % m) > i);
201:           ym = N/n + ((N % n) > j);
202:           zm = P/p + ((P % p) > k);
203:         }

205:         xsize = xm;
206:         ysize = ym;
207:         zsize = zm;
208:         xo = xs;
209:         yo = ys;
210:         zo = zs;

212:         DMDACreate(PETSC_COMM_SELF,&(da[idx]));
213:         DMSetOptionsPrefix(da[idx],"sub_");
214:         DMSetDimension(da[idx], info.dim);
215:         DMDASetDof(da[idx], info.dof);

217:         DMDASetStencilType(da[idx],info.st);
218:         DMDASetStencilWidth(da[idx],info.sw);

220:         if (info.bx == DM_BOUNDARY_PERIODIC || (xs != 0)) {
221:           xsize += xol;
222:           xo    -= xol;
223:         }
224:         if (info.by == DM_BOUNDARY_PERIODIC || (ys != 0)) {
225:           ysize += yol;
226:           yo    -= yol;
227:         }
228:         if (info.bz == DM_BOUNDARY_PERIODIC || (zs != 0)) {
229:           zsize += zol;
230:           zo    -= zol;
231:         }

233:         if (info.bx == DM_BOUNDARY_PERIODIC || (xs+xm != info.mx)) xsize += xol;
234:         if (info.by == DM_BOUNDARY_PERIODIC || (ys+ym != info.my)) ysize += yol;
235:         if (info.bz == DM_BOUNDARY_PERIODIC || (zs+zm != info.mz)) zsize += zol;

237:         if (info.bx != DM_BOUNDARY_PERIODIC) {
238:           if (xo < 0) {
239:             xsize += xo;
240:             xo = 0;
241:           }
242:           if (xo+xsize > info.mx-1) {
243:             xsize -= xo+xsize - info.mx;
244:           }
245:         }
246:         if (info.by != DM_BOUNDARY_PERIODIC) {
247:           if (yo < 0) {
248:             ysize += yo;
249:             yo = 0;
250:           }
251:           if (yo+ysize > info.my-1) {
252:             ysize -= yo+ysize - info.my;
253:           }
254:         }
255:         if (info.bz != DM_BOUNDARY_PERIODIC) {
256:           if (zo < 0) {
257:             zsize += zo;
258:             zo = 0;
259:           }
260:           if (zo+zsize > info.mz-1) {
261:             zsize -= zo+zsize - info.mz;
262:           }
263:         }

265:         DMDASetSizes(da[idx], xsize, ysize, zsize);
266:         DMDASetNumProcs(da[idx], 1, 1, 1);
267:         DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED);

269:         /* set up as a block instead */
270:         DMSetUp(da[idx]);

272:         /* nonoverlapping region */
273:         DMDASetNonOverlappingRegion(da[idx],xs,ys,zs,xm,ym,zm);

275:         /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */
276:         DMDASetOffset(da[idx],xo,yo,zo,info.mx,info.my,info.mz);
277:         xs += xm;
278:         idx++;
279:       }
280:       ys += ym;
281:     }
282:     zs += zm;
283:   }
284:   *sdm = da;
285:   return(0);
286: }

290: /*
291:  Fills the local vector problem on the subdomain from the global problem.

293:  Right now this assumes one subdomain per processor.

295:  */
296: PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat)
297: {
299:   DMDALocalInfo  info,subinfo;
300:   DM             subdm;
301:   MatStencil     upper,lower;
302:   IS             idis,isis,odis,osis,gdis;
303:   Vec            svec,dvec,slvec;
304:   PetscInt       xm,ym,zm,xs,ys,zs;
305:   PetscInt       i;


309:   /* allocate the arrays of scatters */
310:   if (iscat) {PetscMalloc1(nsubdms,iscat);}
311:   if (oscat) {PetscMalloc1(nsubdms,oscat);}
312:   if (lscat) {PetscMalloc1(nsubdms,lscat);}

314:   DMDAGetLocalInfo(dm,&info);
315:   for (i = 0; i < nsubdms; i++) {
316:     subdm = subdms[i];
317:     DMDAGetLocalInfo(subdm,&subinfo);
318:     DMDAGetNonOverlappingRegion(subdm,&xs,&ys,&zs,&xm,&ym,&zm);

320:     /* create the global and subdomain index sets for the inner domain */
321:     lower.i = xs;
322:     lower.j = ys;
323:     lower.k = zs;
324:     upper.i = xs+xm;
325:     upper.j = ys+ym;
326:     upper.k = zs+zm;
327:     DMDACreatePatchIS(dm,&lower,&upper,&idis);
328:     DMDACreatePatchIS(subdm,&lower,&upper,&isis);

330:     /* create the global and subdomain index sets for the outer subdomain */
331:     lower.i = subinfo.xs;
332:     lower.j = subinfo.ys;
333:     lower.k = subinfo.zs;
334:     upper.i = subinfo.xs+subinfo.xm;
335:     upper.j = subinfo.ys+subinfo.ym;
336:     upper.k = subinfo.zs+subinfo.zm;
337:     DMDACreatePatchIS(dm,&lower,&upper,&odis);
338:     DMDACreatePatchIS(subdm,&lower,&upper,&osis);

340:     /* global and subdomain ISes for the local indices of the subdomain */
341:     /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */
342:     lower.i = subinfo.gxs;
343:     lower.j = subinfo.gys;
344:     lower.k = subinfo.gzs;
345:     upper.i = subinfo.gxs+subinfo.gxm;
346:     upper.j = subinfo.gys+subinfo.gym;
347:     upper.k = subinfo.gzs+subinfo.gzm;

349:     DMDACreatePatchIS(dm,&lower,&upper,&gdis);

351:     /* form the scatter */
352:     DMGetGlobalVector(dm,&dvec);
353:     DMGetGlobalVector(subdm,&svec);
354:     DMGetLocalVector(subdm,&slvec);

356:     if (iscat) {VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[i]);}
357:     if (oscat) {VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[i]);}
358:     if (lscat) {VecScatterCreate(dvec,gdis,slvec,NULL,&(*lscat)[i]);}

360:     DMRestoreGlobalVector(dm,&dvec);
361:     DMRestoreGlobalVector(subdm,&svec);
362:     DMRestoreLocalVector(subdm,&slvec);

364:     ISDestroy(&idis);
365:     ISDestroy(&isis);

367:     ISDestroy(&odis);
368:     ISDestroy(&osis);

370:     ISDestroy(&gdis);
371:   }
372:   return(0);
373: }

377: PetscErrorCode DMDASubDomainIS_Private(DM dm,PetscInt n,DM *subdm,IS **iis,IS **ois)
378: {
380:   PetscInt       i;
381:   DMDALocalInfo  info,subinfo;
382:   MatStencil     lower,upper;

385:   DMDAGetLocalInfo(dm,&info);
386:   if (iis) {PetscMalloc1(n,iis);}
387:   if (ois) {PetscMalloc1(n,ois);}

389:   for (i = 0;i < n; i++) {
390:     DMDAGetLocalInfo(subdm[i],&subinfo);
391:     if (iis) {
392:       /* create the inner IS */
393:       lower.i = info.xs;
394:       lower.j = info.ys;
395:       lower.k = info.zs;
396:       upper.i = info.xs+info.xm;
397:       upper.j = info.ys+info.ym;
398:       upper.k = info.zs+info.zm;
399:       DMDACreatePatchIS(dm,&lower,&upper,&(*iis)[i]);
400:     }

402:     if (ois) {
403:       /* create the outer IS */
404:       lower.i = subinfo.xs;
405:       lower.j = subinfo.ys;
406:       lower.k = subinfo.zs;
407:       upper.i = subinfo.xs+subinfo.xm;
408:       upper.j = subinfo.ys+subinfo.ym;
409:       upper.k = subinfo.zs+subinfo.zm;
410:       DMDACreatePatchIS(dm,&lower,&upper,&(*ois)[i]);
411:     }
412:   }
413:   return(0);
414: }

418: PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm)
419: {
421:   DM             *sdm;
422:   PetscInt       n,i;

425:   DMDASubDomainDA_Private(dm,&n,&sdm);
426:   if (names) {
427:     PetscMalloc1(n,names);
428:     for (i=0;i<n;i++) (*names)[i] = 0;
429:   }
430:   DMDASubDomainIS_Private(dm,n,sdm,iis,ois);
431:   if (subdm) *subdm = sdm;
432:   else {
433:     for (i=0;i<n;i++) {
434:       DMDestroy(&sdm[i]);
435:     }
436:   }
437:   if (len) *len = n;
438:   return(0);
439: }