Actual source code: dadd.c

petsc-main 2021-04-20
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:   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
 12: -  offproc - indicate whether the returned IS will contain off process indices

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

 17:   Level: developer

 19: Notes:
 20: This routine always returns an IS on the DMDA's comm, if offproc is set to PETSC_TRUE,
 21: the routine returns an IS with all the indices requested regardless of whether these indices
 22: are present on the requesting rank or not. Thus, it is upon the caller to ensure that
 23: the indices returned in this mode are appropriate. If offproc is set to PETSC_FALSE,
 24: the IS only returns the subset of indices that are present on the requesting rank and there
 25: is no duplication of indices.

 27: .seealso: DMDACreateDomainDecomposition(), DMDACreateDomainDecompositionScatters()
 28: @*/
 29: PetscErrorCode DMDACreatePatchIS(DM da,MatStencil *lower,MatStencil *upper,IS *is, PetscBool offproc)
 30: {
 31:   PetscInt       ms=0,ns=0,ps=0;
 32:   PetscInt       mw=0,nw=0,pw=0;
 33:   PetscInt       me=1,ne=1,pe=1;
 34:   PetscInt       mr=0,nr=0,pr=0;
 35:   PetscInt       ii,jj,kk;
 36:   PetscInt       si,sj,sk;
 37:   PetscInt       i,j,k,l,idx=0;
 38:   PetscInt       base;
 39:   PetscInt       xm=1,ym=1,zm=1;
 40:   PetscInt       ox,oy,oz;
 41:   PetscInt       m,n,p,M,N,P,dof;
 42:   const PetscInt *lx,*ly,*lz;
 43:   PetscInt       nindices;
 44:   PetscInt       *indices;
 45:   DM_DA          *dd = (DM_DA*)da->data;
 46:   PetscBool      skip_i=PETSC_TRUE, skip_j=PETSC_TRUE, skip_k=PETSC_TRUE;
 47:   PetscBool      valid_j=PETSC_FALSE, valid_k=PETSC_FALSE; /* DMDA has at least 1 dimension */

 51:   M = dd->M; N = dd->N; P = dd->P;
 52:   m = dd->m; n = dd->n; p = dd->p;
 53:   dof = dd->w;

 55:   nindices = -1;
 56:   if (PetscLikely(upper->i - lower->i)) {
 57:     nindices = nindices*(upper->i - lower->i);
 58:     skip_i=PETSC_FALSE;
 59:   }
 60:   if (N>1) {
 61:     valid_j = PETSC_TRUE;
 62:     if (PetscLikely(upper->j - lower->j)) {
 63:       nindices = nindices*(upper->j - lower->j);
 64:       skip_j=PETSC_FALSE;
 65:     }
 66:   }
 67:   if (P>1) {
 68:     valid_k = PETSC_TRUE;
 69:     if (PetscLikely(upper->k - lower->k)) {
 70:       nindices = nindices*(upper->k - lower->k);
 71:       skip_k=PETSC_FALSE;
 72:     }
 73:   }
 74:   if (PetscLikely(nindices<0)) {
 75:     if (PetscUnlikely(skip_i && skip_j && skip_k)) {
 76:       nindices = 0;
 77:     } else nindices = nindices*(-1);
 78:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Lower and Upper stencils are identical! Please check inputs.");

 80:   PetscMalloc1(nindices*dof,&indices);
 81:   DMDAGetOffset(da,&ox,&oy,&oz,NULL,NULL,NULL);

 83:   if (!valid_k) {
 84:     k = 0;
 85:     upper->k=0;
 86:     lower->k=0;
 87:   }
 88:   if (!valid_j) {
 89:     j = 0;
 90:     upper->j=0;
 91:     lower->j=0;
 92:   }

 94:   if (offproc) {
 95:     DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
 96:     /* start at index 0 on processor 0 */
 97:     mr = 0;
 98:     nr = 0;
 99:     pr = 0;
100:     ms = 0;
101:     ns = 0;
102:     ps = 0;
103:     if (lx) me = lx[0];
104:     if (ly) ne = ly[0];
105:     if (lz) pe = lz[0];
106:     /*
107:        If no indices are to be returned, create an empty is,
108:        this prevents hanging in while loops
109:     */
110:     if (skip_i && skip_j && skip_k) goto createis;
111:     /*
112:        do..while loops to ensure the block gets entered once,
113:        regardless of control condition being met, necessary for
114:        cases when a subset of skip_i/j/k is true
115:     */
116:     if (skip_k) k = upper->k-oz; else k = lower->k-oz;
117:     do {
118:       if (skip_j) j = upper->j-oy; else j = lower->j-oy;
119:       do {
120:         if (skip_i) i = upper->i-ox; else i = lower->i-ox;
121:         do {
122:           /* "actual" indices rather than ones outside of the domain */
123:           ii = i;
124:           jj = j;
125:           kk = k;
126:           if (ii < 0) ii = M + ii;
127:           if (jj < 0) jj = N + jj;
128:           if (kk < 0) kk = P + kk;
129:           if (ii > M-1) ii = ii - M;
130:           if (jj > N-1) jj = jj - N;
131:           if (kk > P-1) kk = kk - P;
132:           /* gone out of processor range on x axis */
133:           while (ii > me-1 || ii < ms) {
134:             if (mr == m-1) {
135:               ms = 0;
136:               me = lx[0];
137:               mr = 0;
138:             } else {
139:               mr++;
140:               ms = me;
141:               me += lx[mr];
142:             }
143:           }
144:           /* gone out of processor range on y axis */
145:           while (jj > ne-1 || jj < ns) {
146:             if (nr == n-1) {
147:               ns = 0;
148:               ne = ly[0];
149:               nr = 0;
150:             } else {
151:               nr++;
152:               ns = ne;
153:               ne += ly[nr];
154:             }
155:           }
156:           /* gone out of processor range on z axis */
157:           while (kk > pe-1 || kk < ps) {
158:             if (pr == p-1) {
159:               ps = 0;
160:               pe = lz[0];
161:               pr = 0;
162:             } else {
163:               pr++;
164:               ps = pe;
165:               pe += lz[pr];
166:             }
167:           }
168:           /* compute the vector base on owning processor */
169:           xm = me - ms;
170:           ym = ne - ns;
171:           zm = pe - ps;
172:           base = ms*ym*zm + ns*M*zm + ps*M*N;
173:           /* compute the local coordinates on owning processor */
174:           si = ii - ms;
175:           sj = jj - ns;
176:           sk = kk - ps;
177:           for (l=0;l<dof;l++) {
178:             indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk);
179:             idx++;
180:           }
181:           i++;
182:         } while (i<upper->i-ox);
183:         j++;
184:       } while (j<upper->j-oy);
185:       k++;
186:     } while (k<upper->k-oz);
187:   }

189:   if (!offproc){
190:     DMDAGetCorners(da, &ms, &ns, &ps, &mw, &nw, &pw);
191:     me = ms + mw;
192:     if (N>1) ne = ns + nw;
193:     if (P>1) pe = ps + pw;
194:     /* Account for DM offsets */
195:     ms = ms - ox; me = me - ox;
196:     ns = ns - oy; ne = ne - oy;
197:     ps = ps - oz; pe = pe - oz;

199:     /* compute the vector base on owning processor */
200:     xm = me - ms;
201:     ym = ne - ns;
202:     zm = pe - ps;
203:     base = ms*ym*zm + ns*M*zm + ps*M*N;
204:     /*
205:        if no indices are to be returned, create an empty is,
206:        this prevents hanging in while loops
207:     */
208:     if (skip_i && skip_j && skip_k) goto createis;
209:     /*
210:        do..while loops to ensure the block gets entered once,
211:        regardless of control condition being met, necessary for
212:        cases when a subset of skip_i/j/k is true
213:     */
214:     if (skip_k) k = upper->k-oz; else k = lower->k-oz;
215:     do {
216:       if (skip_j) j = upper->j-oy; else j = lower->j-oy;
217:       do {
218:         if (skip_i) i = upper->i-ox; else i = lower->i-ox;
219:         do {
220:           if (k>=ps && k<=pe-1) {
221:             if (j>=ns && j<=ne-1) {
222:               if (i>=ms && i<=me-1) {
223:                 /* compute the local coordinates on owning processor */
224:                 si = i - ms;
225:                 sj = j - ns;
226:                 sk = k - ps;
227:                 for (l=0; l<dof; l++) {
228:                   indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk);
229:                   idx++;
230:                 }
231:               }
232:             }
233:           }
234:           i++;
235:         } while (i<upper->i-ox);
236:         j++;
237:       } while (j<upper->j-oy);
238:       k++;
239:     } while (k<upper->k-oz);

241:     PetscRealloc((size_t)(idx*sizeof(PetscInt)), (void*)&indices);
242:   }

244:   createis:
245:   ISCreateGeneral(PetscObjectComm((PetscObject)da),idx,indices,PETSC_OWN_POINTER,is);
246:   return(0);
247: }

249: PetscErrorCode DMDASubDomainDA_Private(DM dm, PetscInt *nlocal, DM **sdm)
250: {
251:   DM             *da;
252:   PetscInt       dim,size,i,j,k,idx;
254:   DMDALocalInfo  info;
255:   PetscInt       xsize,ysize,zsize;
256:   PetscInt       xo,yo,zo;
257:   PetscInt       xs,ys,zs;
258:   PetscInt       xm=1,ym=1,zm=1;
259:   PetscInt       xol,yol,zol;
260:   PetscInt       m=1,n=1,p=1;
261:   PetscInt       M,N,P;
262:   PetscInt       pm,mtmp;

265:   DMDAGetLocalInfo(dm,&info);
266:   DMDAGetOverlap(dm,&xol,&yol,&zol);
267:   DMDAGetNumLocalSubDomains(dm,&size);
268:   PetscMalloc1(size,&da);

270:   if (nlocal) *nlocal = size;

272:   dim = info.dim;

274:   M = info.xm;
275:   N = info.ym;
276:   P = info.zm;

278:   if (dim == 1) {
279:     m = size;
280:   } else if (dim == 2) {
281:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N)));
282:     while (m > 0) {
283:       n = size/m;
284:       if (m*n*p == size) break;
285:       m--;
286:     }
287:   } else if (dim == 3) {
288:     n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.))); if (!n) n = 1;
289:     while (n > 0) {
290:       pm = size/n;
291:       if (n*pm == size) break;
292:       n--;
293:     }
294:     if (!n) n = 1;
295:     m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
296:     if (!m) m = 1;
297:     while (m > 0) {
298:       p = size/(m*n);
299:       if (m*n*p == size) break;
300:       m--;
301:     }
302:     if (M > P && m < p) {mtmp = m; m = p; p = mtmp;}
303:   }

305:   zs = info.zs;
306:   idx = 0;
307:   for (k = 0; k < p; k++) {
308:     ys = info.ys;
309:     for (j = 0; j < n; j++) {
310:       xs = info.xs;
311:       for (i = 0; i < m; i++) {
312:         if (dim == 1) {
313:           xm = M/m + ((M % m) > i);
314:         } else if (dim == 2) {
315:           xm = M/m + ((M % m) > i);
316:           ym = N/n + ((N % n) > j);
317:         } else if (dim == 3) {
318:           xm = M/m + ((M % m) > i);
319:           ym = N/n + ((N % n) > j);
320:           zm = P/p + ((P % p) > k);
321:         }

323:         xsize = xm;
324:         ysize = ym;
325:         zsize = zm;
326:         xo = xs;
327:         yo = ys;
328:         zo = zs;

330:         DMDACreate(PETSC_COMM_SELF,&(da[idx]));
331:         DMSetOptionsPrefix(da[idx],"sub_");
332:         DMSetDimension(da[idx], info.dim);
333:         DMDASetDof(da[idx], info.dof);

335:         DMDASetStencilType(da[idx],info.st);
336:         DMDASetStencilWidth(da[idx],info.sw);

338:         if (info.bx == DM_BOUNDARY_PERIODIC || (xs != 0)) {
339:           xsize += xol;
340:           xo    -= xol;
341:         }
342:         if (info.by == DM_BOUNDARY_PERIODIC || (ys != 0)) {
343:           ysize += yol;
344:           yo    -= yol;
345:         }
346:         if (info.bz == DM_BOUNDARY_PERIODIC || (zs != 0)) {
347:           zsize += zol;
348:           zo    -= zol;
349:         }

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

355:         if (info.bx != DM_BOUNDARY_PERIODIC) {
356:           if (xo < 0) {
357:             xsize += xo;
358:             xo = 0;
359:           }
360:           if (xo+xsize > info.mx-1) {
361:             xsize -= xo+xsize - info.mx;
362:           }
363:         }
364:         if (info.by != DM_BOUNDARY_PERIODIC) {
365:           if (yo < 0) {
366:             ysize += yo;
367:             yo = 0;
368:           }
369:           if (yo+ysize > info.my-1) {
370:             ysize -= yo+ysize - info.my;
371:           }
372:         }
373:         if (info.bz != DM_BOUNDARY_PERIODIC) {
374:           if (zo < 0) {
375:             zsize += zo;
376:             zo = 0;
377:           }
378:           if (zo+zsize > info.mz-1) {
379:             zsize -= zo+zsize - info.mz;
380:           }
381:         }

383:         DMDASetSizes(da[idx], xsize, ysize, zsize);
384:         DMDASetNumProcs(da[idx], 1, 1, 1);
385:         DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED);

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

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

393:         /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */
394:         DMDASetOffset(da[idx],xo,yo,zo,info.mx,info.my,info.mz);
395:         xs += xm;
396:         idx++;
397:       }
398:       ys += ym;
399:     }
400:     zs += zm;
401:   }
402:   *sdm = da;
403:   return(0);
404: }

406: /*
407:    Fills the local vector problem on the subdomain from the global problem.

409:    Right now this assumes one subdomain per processor.

411: */
412: PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat)
413: {
415:   DMDALocalInfo  info,subinfo;
416:   DM             subdm;
417:   MatStencil     upper,lower;
418:   IS             idis,isis,odis,osis,gdis;
419:   Vec            svec,dvec,slvec;
420:   PetscInt       xm,ym,zm,xs,ys,zs;
421:   PetscInt       i;
422:   PetscBool      patchis_offproc = PETSC_TRUE;

425:   /* allocate the arrays of scatters */
426:   if (iscat) {PetscMalloc1(nsubdms,iscat);}
427:   if (oscat) {PetscMalloc1(nsubdms,oscat);}
428:   if (lscat) {PetscMalloc1(nsubdms,lscat);}

430:   DMDAGetLocalInfo(dm,&info);
431:   for (i = 0; i < nsubdms; i++) {
432:     subdm = subdms[i];
433:     DMDAGetLocalInfo(subdm,&subinfo);
434:     DMDAGetNonOverlappingRegion(subdm,&xs,&ys,&zs,&xm,&ym,&zm);

436:     /* create the global and subdomain index sets for the inner domain */
437:     lower.i = xs;
438:     lower.j = ys;
439:     lower.k = zs;
440:     upper.i = xs+xm;
441:     upper.j = ys+ym;
442:     upper.k = zs+zm;
443:     DMDACreatePatchIS(dm,&lower,&upper,&idis,patchis_offproc);
444:     DMDACreatePatchIS(subdm,&lower,&upper,&isis,patchis_offproc);

446:     /* create the global and subdomain index sets for the outer subdomain */
447:     lower.i = subinfo.xs;
448:     lower.j = subinfo.ys;
449:     lower.k = subinfo.zs;
450:     upper.i = subinfo.xs+subinfo.xm;
451:     upper.j = subinfo.ys+subinfo.ym;
452:     upper.k = subinfo.zs+subinfo.zm;
453:     DMDACreatePatchIS(dm,&lower,&upper,&odis,patchis_offproc);
454:     DMDACreatePatchIS(subdm,&lower,&upper,&osis,patchis_offproc);

456:     /* global and subdomain ISes for the local indices of the subdomain */
457:     /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */
458:     lower.i = subinfo.gxs;
459:     lower.j = subinfo.gys;
460:     lower.k = subinfo.gzs;
461:     upper.i = subinfo.gxs+subinfo.gxm;
462:     upper.j = subinfo.gys+subinfo.gym;
463:     upper.k = subinfo.gzs+subinfo.gzm;
464:     DMDACreatePatchIS(dm,&lower,&upper,&gdis,patchis_offproc);

466:     /* form the scatter */
467:     DMGetGlobalVector(dm,&dvec);
468:     DMGetGlobalVector(subdm,&svec);
469:     DMGetLocalVector(subdm,&slvec);

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

475:     DMRestoreGlobalVector(dm,&dvec);
476:     DMRestoreGlobalVector(subdm,&svec);
477:     DMRestoreLocalVector(subdm,&slvec);

479:     ISDestroy(&idis);
480:     ISDestroy(&isis);

482:     ISDestroy(&odis);
483:     ISDestroy(&osis);

485:     ISDestroy(&gdis);
486:   }
487:   return(0);
488: }

490: PetscErrorCode DMDASubDomainIS_Private(DM dm,PetscInt n,DM *subdm,IS **iis,IS **ois)
491: {
493:   PetscInt       i;
494:   DMDALocalInfo  info,subinfo;
495:   MatStencil     lower,upper;
496:   PetscBool      patchis_offproc = PETSC_TRUE;

499:   DMDAGetLocalInfo(dm,&info);
500:   if (iis) {PetscMalloc1(n,iis);}
501:   if (ois) {PetscMalloc1(n,ois);}

503:   for (i = 0;i < n; i++) {
504:     DMDAGetLocalInfo(subdm[i],&subinfo);
505:     if (iis) {
506:       /* create the inner IS */
507:       lower.i = info.xs;
508:       lower.j = info.ys;
509:       lower.k = info.zs;
510:       upper.i = info.xs+info.xm;
511:       upper.j = info.ys+info.ym;
512:       upper.k = info.zs+info.zm;
513:       DMDACreatePatchIS(dm,&lower,&upper,&(*iis)[i],patchis_offproc);
514:     }

516:     if (ois) {
517:       /* create the outer IS */
518:       lower.i = subinfo.xs;
519:       lower.j = subinfo.ys;
520:       lower.k = subinfo.zs;
521:       upper.i = subinfo.xs+subinfo.xm;
522:       upper.j = subinfo.ys+subinfo.ym;
523:       upper.k = subinfo.zs+subinfo.zm;
524:       DMDACreatePatchIS(dm,&lower,&upper,&(*ois)[i],patchis_offproc);
525:     }
526:   }
527:   return(0);
528: }

530: PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm)
531: {
533:   DM             *sdm;
534:   PetscInt       n,i;

537:   DMDASubDomainDA_Private(dm,&n,&sdm);
538:   if (names) {
539:     PetscMalloc1(n,names);
540:     for (i=0;i<n;i++) (*names)[i] = NULL;
541:   }
542:   DMDASubDomainIS_Private(dm,n,sdm,iis,ois);
543:   if (subdm) *subdm = sdm;
544:   else {
545:     for (i=0;i<n;i++) {
546:       DMDestroy(&sdm[i]);
547:     }
548:   }
549:   if (len) *len = n;
550:   return(0);
551: }