Actual source code: dadd.c
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 */
49: M = dd->M; N = dd->N; P = dd->P;
50: m = dd->m; n = dd->n; p = dd->p;
51: dof = dd->w;
53: nindices = -1;
54: if (PetscLikely(upper->i - lower->i)) {
55: nindices = nindices*(upper->i - lower->i);
56: skip_i=PETSC_FALSE;
57: }
58: if (N>1) {
59: valid_j = PETSC_TRUE;
60: if (PetscLikely(upper->j - lower->j)) {
61: nindices = nindices*(upper->j - lower->j);
62: skip_j=PETSC_FALSE;
63: }
64: }
65: if (P>1) {
66: valid_k = PETSC_TRUE;
67: if (PetscLikely(upper->k - lower->k)) {
68: nindices = nindices*(upper->k - lower->k);
69: skip_k=PETSC_FALSE;
70: }
71: }
72: if (PetscLikely(nindices<0)) {
73: if (PetscUnlikely(skip_i && skip_j && skip_k)) {
74: nindices = 0;
75: } else nindices = nindices*(-1);
76: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Lower and Upper stencils are identical! Please check inputs.");
78: PetscMalloc1(nindices*dof,&indices);
79: DMDAGetOffset(da,&ox,&oy,&oz,NULL,NULL,NULL);
81: if (!valid_k) {
82: k = 0;
83: upper->k=0;
84: lower->k=0;
85: }
86: if (!valid_j) {
87: j = 0;
88: upper->j=0;
89: lower->j=0;
90: }
92: if (offproc) {
93: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
94: /* start at index 0 on processor 0 */
95: mr = 0;
96: nr = 0;
97: pr = 0;
98: ms = 0;
99: ns = 0;
100: ps = 0;
101: if (lx) me = lx[0];
102: if (ly) ne = ly[0];
103: if (lz) pe = lz[0];
104: /*
105: If no indices are to be returned, create an empty is,
106: this prevents hanging in while loops
107: */
108: if (skip_i && skip_j && skip_k) goto createis;
109: /*
110: do..while loops to ensure the block gets entered once,
111: regardless of control condition being met, necessary for
112: cases when a subset of skip_i/j/k is true
113: */
114: if (skip_k) k = upper->k-oz; else k = lower->k-oz;
115: do {
116: if (skip_j) j = upper->j-oy; else j = lower->j-oy;
117: do {
118: if (skip_i) i = upper->i-ox; else i = lower->i-ox;
119: do {
120: /* "actual" indices rather than ones outside of the domain */
121: ii = i;
122: jj = j;
123: kk = k;
124: if (ii < 0) ii = M + ii;
125: if (jj < 0) jj = N + jj;
126: if (kk < 0) kk = P + kk;
127: if (ii > M-1) ii = ii - M;
128: if (jj > N-1) jj = jj - N;
129: if (kk > P-1) kk = kk - P;
130: /* gone out of processor range on x axis */
131: while (ii > me-1 || ii < ms) {
132: if (mr == m-1) {
133: ms = 0;
134: me = lx[0];
135: mr = 0;
136: } else {
137: mr++;
138: ms = me;
139: me += lx[mr];
140: }
141: }
142: /* gone out of processor range on y axis */
143: while (jj > ne-1 || jj < ns) {
144: if (nr == n-1) {
145: ns = 0;
146: ne = ly[0];
147: nr = 0;
148: } else {
149: nr++;
150: ns = ne;
151: ne += ly[nr];
152: }
153: }
154: /* gone out of processor range on z axis */
155: while (kk > pe-1 || kk < ps) {
156: if (pr == p-1) {
157: ps = 0;
158: pe = lz[0];
159: pr = 0;
160: } else {
161: pr++;
162: ps = pe;
163: pe += lz[pr];
164: }
165: }
166: /* compute the vector base on owning processor */
167: xm = me - ms;
168: ym = ne - ns;
169: zm = pe - ps;
170: base = ms*ym*zm + ns*M*zm + ps*M*N;
171: /* compute the local coordinates on owning processor */
172: si = ii - ms;
173: sj = jj - ns;
174: sk = kk - ps;
175: for (l=0;l<dof;l++) {
176: indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk);
177: idx++;
178: }
179: i++;
180: } while (i<upper->i-ox);
181: j++;
182: } while (j<upper->j-oy);
183: k++;
184: } while (k<upper->k-oz);
185: }
187: if (!offproc) {
188: DMDAGetCorners(da, &ms, &ns, &ps, &mw, &nw, &pw);
189: me = ms + mw;
190: if (N>1) ne = ns + nw;
191: if (P>1) pe = ps + pw;
192: /* Account for DM offsets */
193: ms = ms - ox; me = me - ox;
194: ns = ns - oy; ne = ne - oy;
195: ps = ps - oz; pe = pe - oz;
197: /* compute the vector base on owning processor */
198: xm = me - ms;
199: ym = ne - ns;
200: zm = pe - ps;
201: base = ms*ym*zm + ns*M*zm + ps*M*N;
202: /*
203: if no indices are to be returned, create an empty is,
204: this prevents hanging in while loops
205: */
206: if (skip_i && skip_j && skip_k) goto createis;
207: /*
208: do..while loops to ensure the block gets entered once,
209: regardless of control condition being met, necessary for
210: cases when a subset of skip_i/j/k is true
211: */
212: if (skip_k) k = upper->k-oz; else k = lower->k-oz;
213: do {
214: if (skip_j) j = upper->j-oy; else j = lower->j-oy;
215: do {
216: if (skip_i) i = upper->i-ox; else i = lower->i-ox;
217: do {
218: if (k>=ps && k<=pe-1) {
219: if (j>=ns && j<=ne-1) {
220: if (i>=ms && i<=me-1) {
221: /* compute the local coordinates on owning processor */
222: si = i - ms;
223: sj = j - ns;
224: sk = k - ps;
225: for (l=0; l<dof; l++) {
226: indices[idx] = l + dof*(base + si + xm*sj + xm*ym*sk);
227: idx++;
228: }
229: }
230: }
231: }
232: i++;
233: } while (i<upper->i-ox);
234: j++;
235: } while (j<upper->j-oy);
236: k++;
237: } while (k<upper->k-oz);
239: PetscRealloc((size_t)(idx*sizeof(PetscInt)), (void*)&indices);
240: }
242: createis:
243: ISCreateGeneral(PetscObjectComm((PetscObject)da),idx,indices,PETSC_OWN_POINTER,is);
244: return 0;
245: }
247: PetscErrorCode DMDASubDomainDA_Private(DM dm, PetscInt *nlocal, DM **sdm)
248: {
249: DM *da;
250: PetscInt dim,size,i,j,k,idx;
251: DMDALocalInfo info;
252: PetscInt xsize,ysize,zsize;
253: PetscInt xo,yo,zo;
254: PetscInt xs,ys,zs;
255: PetscInt xm=1,ym=1,zm=1;
256: PetscInt xol,yol,zol;
257: PetscInt m=1,n=1,p=1;
258: PetscInt M,N,P;
259: PetscInt pm,mtmp;
261: DMDAGetLocalInfo(dm,&info);
262: DMDAGetOverlap(dm,&xol,&yol,&zol);
263: DMDAGetNumLocalSubDomains(dm,&size);
264: PetscMalloc1(size,&da);
266: if (nlocal) *nlocal = size;
268: dim = info.dim;
270: M = info.xm;
271: N = info.ym;
272: P = info.zm;
274: if (dim == 1) {
275: m = size;
276: } else if (dim == 2) {
277: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N)));
278: while (m > 0) {
279: n = size/m;
280: if (m*n*p == size) break;
281: m--;
282: }
283: } else if (dim == 3) {
284: n = (PetscInt)(0.5 + PetscPowReal(((PetscReal)N*N)*((PetscReal)size)/((PetscReal)P*M),(PetscReal)(1./3.))); if (!n) n = 1;
285: while (n > 0) {
286: pm = size/n;
287: if (n*pm == size) break;
288: n--;
289: }
290: if (!n) n = 1;
291: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)P*n)));
292: if (!m) m = 1;
293: while (m > 0) {
294: p = size/(m*n);
295: if (m*n*p == size) break;
296: m--;
297: }
298: if (M > P && m < p) {mtmp = m; m = p; p = mtmp;}
299: }
301: zs = info.zs;
302: idx = 0;
303: for (k = 0; k < p; k++) {
304: ys = info.ys;
305: for (j = 0; j < n; j++) {
306: xs = info.xs;
307: for (i = 0; i < m; i++) {
308: if (dim == 1) {
309: xm = M/m + ((M % m) > i);
310: } else if (dim == 2) {
311: xm = M/m + ((M % m) > i);
312: ym = N/n + ((N % n) > j);
313: } else if (dim == 3) {
314: xm = M/m + ((M % m) > i);
315: ym = N/n + ((N % n) > j);
316: zm = P/p + ((P % p) > k);
317: }
319: xsize = xm;
320: ysize = ym;
321: zsize = zm;
322: xo = xs;
323: yo = ys;
324: zo = zs;
326: DMDACreate(PETSC_COMM_SELF,&(da[idx]));
327: DMSetOptionsPrefix(da[idx],"sub_");
328: DMSetDimension(da[idx], info.dim);
329: DMDASetDof(da[idx], info.dof);
331: DMDASetStencilType(da[idx],info.st);
332: DMDASetStencilWidth(da[idx],info.sw);
334: if (info.bx == DM_BOUNDARY_PERIODIC || (xs != 0)) {
335: xsize += xol;
336: xo -= xol;
337: }
338: if (info.by == DM_BOUNDARY_PERIODIC || (ys != 0)) {
339: ysize += yol;
340: yo -= yol;
341: }
342: if (info.bz == DM_BOUNDARY_PERIODIC || (zs != 0)) {
343: zsize += zol;
344: zo -= zol;
345: }
347: if (info.bx == DM_BOUNDARY_PERIODIC || (xs+xm != info.mx)) xsize += xol;
348: if (info.by == DM_BOUNDARY_PERIODIC || (ys+ym != info.my)) ysize += yol;
349: if (info.bz == DM_BOUNDARY_PERIODIC || (zs+zm != info.mz)) zsize += zol;
351: if (info.bx != DM_BOUNDARY_PERIODIC) {
352: if (xo < 0) {
353: xsize += xo;
354: xo = 0;
355: }
356: if (xo+xsize > info.mx-1) {
357: xsize -= xo+xsize - info.mx;
358: }
359: }
360: if (info.by != DM_BOUNDARY_PERIODIC) {
361: if (yo < 0) {
362: ysize += yo;
363: yo = 0;
364: }
365: if (yo+ysize > info.my-1) {
366: ysize -= yo+ysize - info.my;
367: }
368: }
369: if (info.bz != DM_BOUNDARY_PERIODIC) {
370: if (zo < 0) {
371: zsize += zo;
372: zo = 0;
373: }
374: if (zo+zsize > info.mz-1) {
375: zsize -= zo+zsize - info.mz;
376: }
377: }
379: DMDASetSizes(da[idx], xsize, ysize, zsize);
380: DMDASetNumProcs(da[idx], 1, 1, 1);
381: DMDASetBoundaryType(da[idx], DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_GHOSTED);
383: /* set up as a block instead */
384: DMSetUp(da[idx]);
386: /* nonoverlapping region */
387: DMDASetNonOverlappingRegion(da[idx],xs,ys,zs,xm,ym,zm);
389: /* this alters the behavior of DMDAGetInfo, DMDAGetLocalInfo, DMDAGetCorners, and DMDAGetGhostedCorners and should be used with care */
390: DMDASetOffset(da[idx],xo,yo,zo,info.mx,info.my,info.mz);
391: xs += xm;
392: idx++;
393: }
394: ys += ym;
395: }
396: zs += zm;
397: }
398: *sdm = da;
399: return 0;
400: }
402: /*
403: Fills the local vector problem on the subdomain from the global problem.
405: Right now this assumes one subdomain per processor.
407: */
408: PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM dm,PetscInt nsubdms,DM *subdms,VecScatter **iscat,VecScatter **oscat, VecScatter **lscat)
409: {
410: DMDALocalInfo info,subinfo;
411: DM subdm;
412: MatStencil upper,lower;
413: IS idis,isis,odis,osis,gdis;
414: Vec svec,dvec,slvec;
415: PetscInt xm,ym,zm,xs,ys,zs;
416: PetscInt i;
417: PetscBool patchis_offproc = PETSC_TRUE;
419: /* allocate the arrays of scatters */
420: if (iscat) PetscMalloc1(nsubdms,iscat);
421: if (oscat) PetscMalloc1(nsubdms,oscat);
422: if (lscat) PetscMalloc1(nsubdms,lscat);
424: DMDAGetLocalInfo(dm,&info);
425: for (i = 0; i < nsubdms; i++) {
426: subdm = subdms[i];
427: DMDAGetLocalInfo(subdm,&subinfo);
428: DMDAGetNonOverlappingRegion(subdm,&xs,&ys,&zs,&xm,&ym,&zm);
430: /* create the global and subdomain index sets for the inner domain */
431: lower.i = xs;
432: lower.j = ys;
433: lower.k = zs;
434: upper.i = xs+xm;
435: upper.j = ys+ym;
436: upper.k = zs+zm;
437: DMDACreatePatchIS(dm,&lower,&upper,&idis,patchis_offproc);
438: DMDACreatePatchIS(subdm,&lower,&upper,&isis,patchis_offproc);
440: /* create the global and subdomain index sets for the outer subdomain */
441: lower.i = subinfo.xs;
442: lower.j = subinfo.ys;
443: lower.k = subinfo.zs;
444: upper.i = subinfo.xs+subinfo.xm;
445: upper.j = subinfo.ys+subinfo.ym;
446: upper.k = subinfo.zs+subinfo.zm;
447: DMDACreatePatchIS(dm,&lower,&upper,&odis,patchis_offproc);
448: DMDACreatePatchIS(subdm,&lower,&upper,&osis,patchis_offproc);
450: /* global and subdomain ISes for the local indices of the subdomain */
451: /* todo - make this not loop over at nonperiodic boundaries, which will be more involved */
452: lower.i = subinfo.gxs;
453: lower.j = subinfo.gys;
454: lower.k = subinfo.gzs;
455: upper.i = subinfo.gxs+subinfo.gxm;
456: upper.j = subinfo.gys+subinfo.gym;
457: upper.k = subinfo.gzs+subinfo.gzm;
458: DMDACreatePatchIS(dm,&lower,&upper,&gdis,patchis_offproc);
460: /* form the scatter */
461: DMGetGlobalVector(dm,&dvec);
462: DMGetGlobalVector(subdm,&svec);
463: DMGetLocalVector(subdm,&slvec);
465: if (iscat) VecScatterCreate(dvec,idis,svec,isis,&(*iscat)[i]);
466: if (oscat) VecScatterCreate(dvec,odis,svec,osis,&(*oscat)[i]);
467: if (lscat) VecScatterCreate(dvec,gdis,slvec,NULL,&(*lscat)[i]);
469: DMRestoreGlobalVector(dm,&dvec);
470: DMRestoreGlobalVector(subdm,&svec);
471: DMRestoreLocalVector(subdm,&slvec);
473: ISDestroy(&idis);
474: ISDestroy(&isis);
476: ISDestroy(&odis);
477: ISDestroy(&osis);
479: ISDestroy(&gdis);
480: }
481: return 0;
482: }
484: PetscErrorCode DMDASubDomainIS_Private(DM dm,PetscInt n,DM *subdm,IS **iis,IS **ois)
485: {
486: PetscInt i;
487: DMDALocalInfo info,subinfo;
488: MatStencil lower,upper;
489: PetscBool patchis_offproc = PETSC_TRUE;
491: DMDAGetLocalInfo(dm,&info);
492: if (iis) PetscMalloc1(n,iis);
493: if (ois) PetscMalloc1(n,ois);
495: for (i = 0;i < n; i++) {
496: DMDAGetLocalInfo(subdm[i],&subinfo);
497: if (iis) {
498: /* create the inner IS */
499: lower.i = info.xs;
500: lower.j = info.ys;
501: lower.k = info.zs;
502: upper.i = info.xs+info.xm;
503: upper.j = info.ys+info.ym;
504: upper.k = info.zs+info.zm;
505: DMDACreatePatchIS(dm,&lower,&upper,&(*iis)[i],patchis_offproc);
506: }
508: if (ois) {
509: /* create the outer IS */
510: lower.i = subinfo.xs;
511: lower.j = subinfo.ys;
512: lower.k = subinfo.zs;
513: upper.i = subinfo.xs+subinfo.xm;
514: upper.j = subinfo.ys+subinfo.ym;
515: upper.k = subinfo.zs+subinfo.zm;
516: DMDACreatePatchIS(dm,&lower,&upper,&(*ois)[i],patchis_offproc);
517: }
518: }
519: return 0;
520: }
522: PetscErrorCode DMCreateDomainDecomposition_DA(DM dm,PetscInt *len,char ***names,IS **iis,IS **ois,DM **subdm)
523: {
524: DM *sdm;
525: PetscInt n,i;
527: DMDASubDomainDA_Private(dm,&n,&sdm);
528: if (names) {
529: PetscMalloc1(n,names);
530: for (i=0;i<n;i++) (*names)[i] = NULL;
531: }
532: DMDASubDomainIS_Private(dm,n,sdm,iis,ois);
533: if (subdm) *subdm = sdm;
534: else {
535: for (i=0;i<n;i++) {
536: DMDestroy(&sdm[i]);
537: }
538: }
539: if (len) *len = n;
540: return 0;
541: }