Actual source code: dadd.c
petsc-3.6.1 2015-08-06
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: }