Actual source code: dadd.c
petsc-3.11.4 2019-09-28
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: }