Actual source code: dadd.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
  1:  #include <petsc/private/dmdaimpl.h>

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

  6:   Not Collective

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

 13:   Output Parameters:
 14: .  is - the IS corresponding to the patch

 16:   Level: developer

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

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

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

143:   DMDAGetLocalInfo(dm,&info);
144:   DMDAGetOverlap(dm,&xol,&yol,&zol);
145:   DMDAGetNumLocalSubDomains(dm,&size);
146:   PetscMalloc1(size,&da);

148:   if (nlocal) *nlocal = size;

150:   dim = info.dim;

152:   M = info.xm;
153:   N = info.ym;
154:   P = info.zm;

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

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

201:         xsize = xm;
202:         ysize = ym;
203:         zsize = zm;
204:         xo = xs;
205:         yo = ys;
206:         zo = zs;

208:         DMDACreate(PETSC_COMM_SELF,&(da[idx]));
209:         DMSetOptionsPrefix(da[idx],"sub_");
210:         DMSetDimension(da[idx], info.dim);
211:         DMDASetDof(da[idx], info.dof);

213:         DMDASetStencilType(da[idx],info.st);
214:         DMDASetStencilWidth(da[idx],info.sw);

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

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

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

261:         DMDASetSizes(da[idx], xsize, ysize, zsize);
262:         DMDASetNumProcs(da[idx], 1, 1, 1);
263:         DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED);

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

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

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

284: /*
285:  Fills the local vector problem on the subdomain from the global problem.

287:  Right now this assumes one subdomain per processor.

289:  */
290: PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat)
291: {
293:   DMDALocalInfo  info,subinfo;
294:   DM             subdm;
295:   MatStencil     upper,lower;
296:   IS             idis,isis,odis,osis,gdis;
297:   Vec            svec,dvec,slvec;
298:   PetscInt       xm,ym,zm,xs,ys,zs;
299:   PetscInt       i;


303:   /* allocate the arrays of scatters */
304:   if (iscat) {PetscMalloc1(nsubdms,iscat);}
305:   if (oscat) {PetscMalloc1(nsubdms,oscat);}
306:   if (lscat) {PetscMalloc1(nsubdms,lscat);}

308:   DMDAGetLocalInfo(dm,&info);
309:   for (i = 0; i < nsubdms; i++) {
310:     subdm = subdms[i];
311:     DMDAGetLocalInfo(subdm,&subinfo);
312:     DMDAGetNonOverlappingRegion(subdm,&xs,&ys,&zs,&xm,&ym,&zm);

314:     /* create the global and subdomain index sets for the inner domain */
315:     lower.i = xs;
316:     lower.j = ys;
317:     lower.k = zs;
318:     upper.i = xs+xm;
319:     upper.j = ys+ym;
320:     upper.k = zs+zm;
321:     DMDACreatePatchIS(dm,&lower,&upper,&idis);
322:     DMDACreatePatchIS(subdm,&lower,&upper,&isis);

324:     /* create the global and subdomain index sets for the outer subdomain */
325:     lower.i = subinfo.xs;
326:     lower.j = subinfo.ys;
327:     lower.k = subinfo.zs;
328:     upper.i = subinfo.xs+subinfo.xm;
329:     upper.j = subinfo.ys+subinfo.ym;
330:     upper.k = subinfo.zs+subinfo.zm;
331:     DMDACreatePatchIS(dm,&lower,&upper,&odis);
332:     DMDACreatePatchIS(subdm,&lower,&upper,&osis);

334:     /* global and subdomain ISes for the local indices of the subdomain */
335:     /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */
336:     lower.i = subinfo.gxs;
337:     lower.j = subinfo.gys;
338:     lower.k = subinfo.gzs;
339:     upper.i = subinfo.gxs+subinfo.gxm;
340:     upper.j = subinfo.gys+subinfo.gym;
341:     upper.k = subinfo.gzs+subinfo.gzm;

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

345:     /* form the scatter */
346:     DMGetGlobalVector(dm,&dvec);
347:     DMGetGlobalVector(subdm,&svec);
348:     DMGetLocalVector(subdm,&slvec);

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

354:     DMRestoreGlobalVector(dm,&dvec);
355:     DMRestoreGlobalVector(subdm,&svec);
356:     DMRestoreLocalVector(subdm,&slvec);

358:     ISDestroy(&idis);
359:     ISDestroy(&isis);

361:     ISDestroy(&odis);
362:     ISDestroy(&osis);

364:     ISDestroy(&gdis);
365:   }
366:   return(0);
367: }

369: PetscErrorCode DMDASubDomainIS_Private(DM dm,PetscInt n,DM *subdm,IS **iis,IS **ois)
370: {
372:   PetscInt       i;
373:   DMDALocalInfo  info,subinfo;
374:   MatStencil     lower,upper;

377:   DMDAGetLocalInfo(dm,&info);
378:   if (iis) {PetscMalloc1(n,iis);}
379:   if (ois) {PetscMalloc1(n,ois);}

381:   for (i = 0;i < n; i++) {
382:     DMDAGetLocalInfo(subdm[i],&subinfo);
383:     if (iis) {
384:       /* create the inner IS */
385:       lower.i = info.xs;
386:       lower.j = info.ys;
387:       lower.k = info.zs;
388:       upper.i = info.xs+info.xm;
389:       upper.j = info.ys+info.ym;
390:       upper.k = info.zs+info.zm;
391:       DMDACreatePatchIS(dm,&lower,&upper,&(*iis)[i]);
392:     }

394:     if (ois) {
395:       /* create the outer IS */
396:       lower.i = subinfo.xs;
397:       lower.j = subinfo.ys;
398:       lower.k = subinfo.zs;
399:       upper.i = subinfo.xs+subinfo.xm;
400:       upper.j = subinfo.ys+subinfo.ym;
401:       upper.k = subinfo.zs+subinfo.zm;
402:       DMDACreatePatchIS(dm,&lower,&upper,&(*ois)[i]);
403:     }
404:   }
405:   return(0);
406: }

408: PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm)
409: {
411:   DM             *sdm;
412:   PetscInt       n,i;

415:   DMDASubDomainDA_Private(dm,&n,&sdm);
416:   if (names) {
417:     PetscMalloc1(n,names);
418:     for (i=0;i<n;i++) (*names)[i] = 0;
419:   }
420:   DMDASubDomainIS_Private(dm,n,sdm,iis,ois);
421:   if (subdm) *subdm = sdm;
422:   else {
423:     for (i=0;i<n;i++) {
424:       DMDestroy(&sdm[i]);
425:     }
426:   }
427:   if (len) *len = n;
428:   return(0);
429: }