Actual source code: da2.c
petsc-3.8.4 2018-03-24
2: #include <petsc/private/dmdaimpl.h>
3: #include <petscdraw.h>
5: static PetscErrorCode DMView_DA_2d(DM da,PetscViewer viewer)
6: {
8: PetscMPIInt rank;
9: PetscBool iascii,isdraw,isglvis,isbinary;
10: DM_DA *dd = (DM_DA*)da->data;
11: #if defined(PETSC_HAVE_MATLAB_ENGINE)
12: PetscBool ismatlab;
13: #endif
16: MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
18: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
19: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
20: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
21: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
22: #if defined(PETSC_HAVE_MATLAB_ENGINE)
23: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
24: #endif
25: if (iascii) {
26: PetscViewerFormat format;
28: PetscViewerGetFormat(viewer, &format);
29: if (format == PETSC_VIEWER_LOAD_BALANCE) {
30: PetscInt i,nmax = 0,nmin = PETSC_MAX_INT,navg = 0,*nz,nzlocal;
31: DMDALocalInfo info;
32: PetscMPIInt size;
33: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
34: DMDAGetLocalInfo(da,&info);
35: nzlocal = info.xm*info.ym;
36: PetscMalloc1(size,&nz);
37: MPI_Allgather(&nzlocal,1,MPIU_INT,nz,1,MPIU_INT,PetscObjectComm((PetscObject)da));
38: for (i=0; i<(PetscInt)size; i++) {
39: nmax = PetscMax(nmax,nz[i]);
40: nmin = PetscMin(nmin,nz[i]);
41: navg += nz[i];
42: }
43: PetscFree(nz);
44: navg = navg/size;
45: PetscViewerASCIIPrintf(viewer," Load Balance - Grid Points: Min %D avg %D max %D\n",nmin,navg,nmax);
46: return(0);
47: }
48: if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL && format != PETSC_VIEWER_ASCII_GLVIS) {
49: DMDALocalInfo info;
50: DMDAGetLocalInfo(da,&info);
51: PetscViewerASCIIPushSynchronized(viewer);
52: PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D m %D n %D w %D s %D\n",rank,dd->M,dd->N,dd->m,dd->n,dd->w,dd->s);
53: PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D\n",info.xs,info.xs+info.xm,info.ys,info.ys+info.ym);
54: PetscViewerFlush(viewer);
55: PetscViewerASCIIPopSynchronized(viewer);
56: } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
57: DMView_DA_GLVis(da,viewer);
58: } else {
59: DMView_DA_VTK(da,viewer);
60: }
61: } else if (isdraw) {
62: PetscDraw draw;
63: double ymin = -1*dd->s-1,ymax = dd->N+dd->s;
64: double xmin = -1*dd->s-1,xmax = dd->M+dd->s;
65: double x,y;
66: PetscInt base;
67: const PetscInt *idx;
68: char node[10];
69: PetscBool isnull;
71: PetscViewerDrawGetDraw(viewer,0,&draw);
72: PetscDrawIsNull(draw,&isnull);
73: if (isnull) return(0);
75: PetscDrawCheckResizedWindow(draw);
76: PetscDrawClear(draw);
77: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
79: PetscDrawCollectiveBegin(draw);
80: /* first processor draw all node lines */
81: if (!rank) {
82: ymin = 0.0; ymax = dd->N - 1;
83: for (xmin=0; xmin<dd->M; xmin++) {
84: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
85: }
86: xmin = 0.0; xmax = dd->M - 1;
87: for (ymin=0; ymin<dd->N; ymin++) {
88: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
89: }
90: }
91: PetscDrawCollectiveEnd(draw);
92: PetscDrawFlush(draw);
93: PetscDrawPause(draw);
95: PetscDrawCollectiveBegin(draw);
96: /* draw my box */
97: xmin = dd->xs/dd->w; xmax =(dd->xe-1)/dd->w; ymin = dd->ys; ymax = dd->ye - 1;
98: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
99: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
100: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
101: PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);
102: /* put in numbers */
103: base = (dd->base)/dd->w;
104: for (y=ymin; y<=ymax; y++) {
105: for (x=xmin; x<=xmax; x++) {
106: PetscSNPrintf(node,sizeof(node),"%d",(int)base++);
107: PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
108: }
109: }
110: PetscDrawCollectiveEnd(draw);
111: PetscDrawFlush(draw);
112: PetscDrawPause(draw);
114: PetscDrawCollectiveBegin(draw);
115: /* overlay ghost numbers, useful for error checking */
116: ISLocalToGlobalMappingGetBlockIndices(da->ltogmap,&idx);
117: base = 0; xmin = dd->Xs; xmax = dd->Xe; ymin = dd->Ys; ymax = dd->Ye;
118: for (y=ymin; y<ymax; y++) {
119: for (x=xmin; x<xmax; x++) {
120: if ((base % dd->w) == 0) {
121: PetscSNPrintf(node,sizeof(node),"%d",(int)(idx[base/dd->w]));
122: PetscDrawString(draw,x/dd->w,y,PETSC_DRAW_BLUE,node);
123: }
124: base++;
125: }
126: }
127: ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap,&idx);
128: PetscDrawCollectiveEnd(draw);
129: PetscDrawFlush(draw);
130: PetscDrawPause(draw);
131: PetscDrawSave(draw);
132: } else if (isglvis) {
133: DMView_DA_GLVis(da,viewer);
134: } else if (isbinary) {
135: DMView_DA_Binary(da,viewer);
136: #if defined(PETSC_HAVE_MATLAB_ENGINE)
137: } else if (ismatlab) {
138: DMView_DA_Matlab(da,viewer);
139: #endif
140: }
141: return(0);
142: }
144: /*
145: M is number of grid points
146: m is number of processors
148: */
149: PetscErrorCode DMDASplitComm2d(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt sw,MPI_Comm *outcomm)
150: {
152: PetscInt m,n = 0,x = 0,y = 0;
153: PetscMPIInt size,csize,rank;
156: MPI_Comm_size(comm,&size);
157: MPI_Comm_rank(comm,&rank);
159: csize = 4*size;
160: do {
161: if (csize % 4) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot split communicator of size %d tried %d %D %D",size,csize,x,y);
162: csize = csize/4;
164: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)csize)/((PetscReal)N)));
165: if (!m) m = 1;
166: while (m > 0) {
167: n = csize/m;
168: if (m*n == csize) break;
169: m--;
170: }
171: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
173: x = M/m + ((M % m) > ((csize-1) % m));
174: y = (N + (csize-1)/m)/n;
175: } while ((x < 4 || y < 4) && csize > 1);
176: if (size != csize) {
177: MPI_Group entire_group,sub_group;
178: PetscMPIInt i,*groupies;
180: MPI_Comm_group(comm,&entire_group);
181: PetscMalloc1(csize,&groupies);
182: for (i=0; i<csize; i++) {
183: groupies[i] = (rank/csize)*csize + i;
184: }
185: MPI_Group_incl(entire_group,csize,groupies,&sub_group);
186: PetscFree(groupies);
187: MPI_Comm_create(comm,sub_group,outcomm);
188: MPI_Group_free(&entire_group);
189: MPI_Group_free(&sub_group);
190: PetscInfo1(0,"DMDASplitComm2d:Creating redundant coarse problems of size %d\n",csize);
191: } else {
192: *outcomm = comm;
193: }
194: return(0);
195: }
197: #if defined(new)
198: /*
199: DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local
200: function lives on a DMDA
202: y ~= (F(u + ha) - F(u))/h,
203: where F = nonlinear function, as set by SNESSetFunction()
204: u = current iterate
205: h = difference interval
206: */
207: PetscErrorCode DMDAGetDiagonal_MFFD(DM da,Vec U,Vec a)
208: {
209: PetscScalar h,*aa,*ww,v;
210: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
212: PetscInt gI,nI;
213: MatStencil stencil;
214: DMDALocalInfo info;
217: (*ctx->func)(0,U,a,ctx->funcctx);
218: (*ctx->funcisetbase)(U,ctx->funcctx);
220: VecGetArray(U,&ww);
221: VecGetArray(a,&aa);
223: nI = 0;
224: h = ww[gI];
225: if (h == 0.0) h = 1.0;
226: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
227: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
228: h *= epsilon;
230: ww[gI] += h;
231: (*ctx->funci)(i,w,&v,ctx->funcctx);
232: aa[nI] = (v - aa[nI])/h;
233: ww[gI] -= h;
234: nI++;
236: VecRestoreArray(U,&ww);
237: VecRestoreArray(a,&aa);
238: return(0);
239: }
240: #endif
242: PetscErrorCode DMSetUp_DA_2D(DM da)
243: {
244: DM_DA *dd = (DM_DA*)da->data;
245: const PetscInt M = dd->M;
246: const PetscInt N = dd->N;
247: PetscInt m = dd->m;
248: PetscInt n = dd->n;
249: const PetscInt dof = dd->w;
250: const PetscInt s = dd->s;
251: DMBoundaryType bx = dd->bx;
252: DMBoundaryType by = dd->by;
253: DMDAStencilType stencil_type = dd->stencil_type;
254: PetscInt *lx = dd->lx;
255: PetscInt *ly = dd->ly;
256: MPI_Comm comm;
257: PetscMPIInt rank,size;
258: PetscInt xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,IXs,IXe,IYs,IYe;
259: PetscInt up,down,left,right,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn;
260: PetscInt xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count;
261: PetscInt s_x,s_y; /* s proportionalized to w */
262: PetscInt sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0;
263: Vec local,global;
264: VecScatter gtol;
265: IS to,from;
266: PetscErrorCode ierr;
269: if (stencil_type == DMDA_STENCIL_BOX && (bx == DM_BOUNDARY_MIRROR || by == DM_BOUNDARY_MIRROR)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Mirror boundary and box stencil");
270: PetscObjectGetComm((PetscObject)da,&comm);
271: #if !defined(PETSC_USE_64BIT_INDICES)
272: if (((PetscInt64) M)*((PetscInt64) N)*((PetscInt64) dof) > (PetscInt64) PETSC_MPI_INT_MAX) SETERRQ3(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D (dof) is too large for 32 bit indices",M,N,dof);
273: #endif
275: if (dof < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
276: if (s < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);
278: MPI_Comm_size(comm,&size);
279: MPI_Comm_rank(comm,&rank);
281: dd->p = 1;
282: if (m != PETSC_DECIDE) {
283: if (m < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
284: else if (m > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
285: }
286: if (n != PETSC_DECIDE) {
287: if (n < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
288: else if (n > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
289: }
291: if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
292: if (n != PETSC_DECIDE) {
293: m = size/n;
294: } else if (m != PETSC_DECIDE) {
295: n = size/m;
296: } else {
297: /* try for squarish distribution */
298: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N)));
299: if (!m) m = 1;
300: while (m > 0) {
301: n = size/m;
302: if (m*n == size) break;
303: m--;
304: }
305: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
306: }
307: if (m*n != size) SETERRQ(comm,PETSC_ERR_PLIB,"Unable to create partition, check the size of the communicator and input m and n ");
308: } else if (m*n != size) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");
310: if (M < m) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
311: if (N < n) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);
313: /*
314: Determine locally owned region
315: xs is the first local node number, x is the number of local nodes
316: */
317: if (!lx) {
318: PetscMalloc1(m, &dd->lx);
319: lx = dd->lx;
320: for (i=0; i<m; i++) {
321: lx[i] = M/m + ((M % m) > i);
322: }
323: }
324: x = lx[rank % m];
325: xs = 0;
326: for (i=0; i<(rank % m); i++) {
327: xs += lx[i];
328: }
329: #if defined(PETSC_USE_DEBUG)
330: left = xs;
331: for (i=(rank % m); i<m; i++) {
332: left += lx[i];
333: }
334: if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M);
335: #endif
337: /*
338: Determine locally owned region
339: ys is the first local node number, y is the number of local nodes
340: */
341: if (!ly) {
342: PetscMalloc1(n, &dd->ly);
343: ly = dd->ly;
344: for (i=0; i<n; i++) {
345: ly[i] = N/n + ((N % n) > i);
346: }
347: }
348: y = ly[rank/m];
349: ys = 0;
350: for (i=0; i<(rank/m); i++) {
351: ys += ly[i];
352: }
353: #if defined(PETSC_USE_DEBUG)
354: left = ys;
355: for (i=(rank/m); i<n; i++) {
356: left += ly[i];
357: }
358: if (left != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N);
359: #endif
361: /*
362: check if the scatter requires more than one process neighbor or wraps around
363: the domain more than once
364: */
365: if ((x < s) && ((m > 1) || (bx == DM_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local x-width of domain x %D is smaller than stencil width s %D",x,s);
366: if ((y < s) && ((n > 1) || (by == DM_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local y-width of domain y %D is smaller than stencil width s %D",y,s);
367: xe = xs + x;
368: ye = ys + y;
370: /* determine ghost region (Xs) and region scattered into (IXs) */
371: if (xs-s > 0) {
372: Xs = xs - s; IXs = xs - s;
373: } else {
374: if (bx) {
375: Xs = xs - s;
376: } else {
377: Xs = 0;
378: }
379: IXs = 0;
380: }
381: if (xe+s <= M) {
382: Xe = xe + s; IXe = xe + s;
383: } else {
384: if (bx) {
385: Xs = xs - s; Xe = xe + s;
386: } else {
387: Xe = M;
388: }
389: IXe = M;
390: }
392: if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
393: IXs = xs - s;
394: IXe = xe + s;
395: Xs = xs - s;
396: Xe = xe + s;
397: }
399: if (ys-s > 0) {
400: Ys = ys - s; IYs = ys - s;
401: } else {
402: if (by) {
403: Ys = ys - s;
404: } else {
405: Ys = 0;
406: }
407: IYs = 0;
408: }
409: if (ye+s <= N) {
410: Ye = ye + s; IYe = ye + s;
411: } else {
412: if (by) {
413: Ye = ye + s;
414: } else {
415: Ye = N;
416: }
417: IYe = N;
418: }
420: if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) {
421: IYs = ys - s;
422: IYe = ye + s;
423: Ys = ys - s;
424: Ye = ye + s;
425: }
427: /* stencil length in each direction */
428: s_x = s;
429: s_y = s;
431: /* determine starting point of each processor */
432: nn = x*y;
433: PetscMalloc2(size+1,&bases,size,&ldims);
434: MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
435: bases[0] = 0;
436: for (i=1; i<=size; i++) {
437: bases[i] = ldims[i-1];
438: }
439: for (i=1; i<=size; i++) {
440: bases[i] += bases[i-1];
441: }
442: base = bases[rank]*dof;
444: /* allocate the base parallel and sequential vectors */
445: dd->Nlocal = x*y*dof;
446: VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,NULL,&global);
447: dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof;
448: VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,NULL,&local);
450: /* generate appropriate vector scatters */
451: /* local to global inserts non-ghost point region into global */
452: PetscMalloc1((IXe-IXs)*(IYe-IYs),&idx);
453: left = xs - Xs; right = left + x;
454: down = ys - Ys; up = down + y;
455: count = 0;
456: for (i=down; i<up; i++) {
457: for (j=left; j<right; j++) {
458: idx[count++] = i*(Xe-Xs) + j;
459: }
460: }
462: /* global to local must include ghost points within the domain,
463: but not ghost points outside the domain that aren't periodic */
464: if (stencil_type == DMDA_STENCIL_BOX) {
465: left = IXs - Xs; right = left + (IXe-IXs);
466: down = IYs - Ys; up = down + (IYe-IYs);
467: count = 0;
468: for (i=down; i<up; i++) {
469: for (j=left; j<right; j++) {
470: idx[count++] = j + i*(Xe-Xs);
471: }
472: }
473: ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
475: } else {
476: /* must drop into cross shape region */
477: /* ---------|
478: | top |
479: |--- ---| up
480: | middle |
481: | |
482: ---- ---- down
483: | bottom |
484: -----------
485: Xs xs xe Xe */
486: left = xs - Xs; right = left + x;
487: down = ys - Ys; up = down + y;
488: count = 0;
489: /* bottom */
490: for (i=(IYs-Ys); i<down; i++) {
491: for (j=left; j<right; j++) {
492: idx[count++] = j + i*(Xe-Xs);
493: }
494: }
495: /* middle */
496: for (i=down; i<up; i++) {
497: for (j=(IXs-Xs); j<(IXe-Xs); j++) {
498: idx[count++] = j + i*(Xe-Xs);
499: }
500: }
501: /* top */
502: for (i=up; i<up+IYe-ye; i++) {
503: for (j=left; j<right; j++) {
504: idx[count++] = j + i*(Xe-Xs);
505: }
506: }
507: ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
508: }
511: /* determine who lies on each side of us stored in n6 n7 n8
512: n3 n5
513: n0 n1 n2
514: */
516: /* Assume the Non-Periodic Case */
517: n1 = rank - m;
518: if (rank % m) {
519: n0 = n1 - 1;
520: } else {
521: n0 = -1;
522: }
523: if ((rank+1) % m) {
524: n2 = n1 + 1;
525: n5 = rank + 1;
526: n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
527: } else {
528: n2 = -1; n5 = -1; n8 = -1;
529: }
530: if (rank % m) {
531: n3 = rank - 1;
532: n6 = n3 + m; if (n6 >= m*n) n6 = -1;
533: } else {
534: n3 = -1; n6 = -1;
535: }
536: n7 = rank + m; if (n7 >= m*n) n7 = -1;
538: if (bx == DM_BOUNDARY_PERIODIC && by == DM_BOUNDARY_PERIODIC) {
539: /* Modify for Periodic Cases */
540: /* Handle all four corners */
541: if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
542: if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
543: if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
544: if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;
546: /* Handle Top and Bottom Sides */
547: if (n1 < 0) n1 = rank + m * (n-1);
548: if (n7 < 0) n7 = rank - m * (n-1);
549: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
550: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
551: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
552: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
554: /* Handle Left and Right Sides */
555: if (n3 < 0) n3 = rank + (m-1);
556: if (n5 < 0) n5 = rank - (m-1);
557: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
558: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
559: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
560: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
561: } else if (by == DM_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */
562: if (n1 < 0) n1 = rank + m * (n-1);
563: if (n7 < 0) n7 = rank - m * (n-1);
564: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
565: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
566: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
567: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
568: } else if (bx == DM_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */
569: if (n3 < 0) n3 = rank + (m-1);
570: if (n5 < 0) n5 = rank - (m-1);
571: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
572: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
573: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
574: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
575: }
577: PetscMalloc1(9,&dd->neighbors);
579: dd->neighbors[0] = n0;
580: dd->neighbors[1] = n1;
581: dd->neighbors[2] = n2;
582: dd->neighbors[3] = n3;
583: dd->neighbors[4] = rank;
584: dd->neighbors[5] = n5;
585: dd->neighbors[6] = n6;
586: dd->neighbors[7] = n7;
587: dd->neighbors[8] = n8;
589: if (stencil_type == DMDA_STENCIL_STAR) {
590: /* save corner processor numbers */
591: sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8;
592: n0 = n2 = n6 = n8 = -1;
593: }
595: PetscMalloc1((Xe-Xs)*(Ye-Ys),&idx);
597: nn = 0;
598: xbase = bases[rank];
599: for (i=1; i<=s_y; i++) {
600: if (n0 >= 0) { /* left below */
601: x_t = lx[n0 % m];
602: y_t = ly[(n0/m)];
603: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
604: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
605: }
607: if (n1 >= 0) { /* directly below */
608: x_t = x;
609: y_t = ly[(n1/m)];
610: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
611: for (j=0; j<x_t; j++) idx[nn++] = s_t++;
612: } else if (by == DM_BOUNDARY_MIRROR) {
613: for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(s_y - i + 1) + j;
614: }
616: if (n2 >= 0) { /* right below */
617: x_t = lx[n2 % m];
618: y_t = ly[(n2/m)];
619: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
620: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
621: }
622: }
624: for (i=0; i<y; i++) {
625: if (n3 >= 0) { /* directly left */
626: x_t = lx[n3 % m];
627: /* y_t = y; */
628: s_t = bases[n3] + (i+1)*x_t - s_x;
629: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
630: } else if (bx == DM_BOUNDARY_MIRROR) {
631: for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*i + s_x - j;
632: }
634: for (j=0; j<x; j++) idx[nn++] = xbase++; /* interior */
636: if (n5 >= 0) { /* directly right */
637: x_t = lx[n5 % m];
638: /* y_t = y; */
639: s_t = bases[n5] + (i)*x_t;
640: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
641: } else if (bx == DM_BOUNDARY_MIRROR) {
642: for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*(i + 1) - 2 - j;
643: }
644: }
646: for (i=1; i<=s_y; i++) {
647: if (n6 >= 0) { /* left above */
648: x_t = lx[n6 % m];
649: /* y_t = ly[(n6/m)]; */
650: s_t = bases[n6] + (i)*x_t - s_x;
651: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
652: }
654: if (n7 >= 0) { /* directly above */
655: x_t = x;
656: /* y_t = ly[(n7/m)]; */
657: s_t = bases[n7] + (i-1)*x_t;
658: for (j=0; j<x_t; j++) idx[nn++] = s_t++;
659: } else if (by == DM_BOUNDARY_MIRROR) {
660: for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(y - i - 1) + j;
661: }
663: if (n8 >= 0) { /* right above */
664: x_t = lx[n8 % m];
665: /* y_t = ly[(n8/m)]; */
666: s_t = bases[n8] + (i-1)*x_t;
667: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
668: }
669: }
671: ISCreateBlock(comm,dof,nn,idx,PETSC_USE_POINTER,&from);
672: VecScatterCreate(global,from,local,to,>ol);
673: PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);
674: ISDestroy(&to);
675: ISDestroy(&from);
677: if (stencil_type == DMDA_STENCIL_STAR) {
678: n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
679: }
681: if (((stencil_type == DMDA_STENCIL_STAR) || (bx && bx != DM_BOUNDARY_PERIODIC) || (by && by != DM_BOUNDARY_PERIODIC))) {
682: /*
683: Recompute the local to global mappings, this time keeping the
684: information about the cross corner processor numbers and any ghosted
685: but not periodic indices.
686: */
687: nn = 0;
688: xbase = bases[rank];
689: for (i=1; i<=s_y; i++) {
690: if (n0 >= 0) { /* left below */
691: x_t = lx[n0 % m];
692: y_t = ly[(n0/m)];
693: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
694: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
695: } else if (xs-Xs > 0 && ys-Ys > 0) {
696: for (j=0; j<s_x; j++) idx[nn++] = -1;
697: }
698: if (n1 >= 0) { /* directly below */
699: x_t = x;
700: y_t = ly[(n1/m)];
701: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
702: for (j=0; j<x_t; j++) idx[nn++] = s_t++;
703: } else if (ys-Ys > 0) {
704: if (by == DM_BOUNDARY_MIRROR) {
705: for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(s_y - i + 1) + j;
706: } else {
707: for (j=0; j<x; j++) idx[nn++] = -1;
708: }
709: }
710: if (n2 >= 0) { /* right below */
711: x_t = lx[n2 % m];
712: y_t = ly[(n2/m)];
713: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
714: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
715: } else if (Xe-xe> 0 && ys-Ys > 0) {
716: for (j=0; j<s_x; j++) idx[nn++] = -1;
717: }
718: }
720: for (i=0; i<y; i++) {
721: if (n3 >= 0) { /* directly left */
722: x_t = lx[n3 % m];
723: /* y_t = y; */
724: s_t = bases[n3] + (i+1)*x_t - s_x;
725: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
726: } else if (xs-Xs > 0) {
727: if (bx == DM_BOUNDARY_MIRROR) {
728: for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*i + s_x - j;
729: } else {
730: for (j=0; j<s_x; j++) idx[nn++] = -1;
731: }
732: }
734: for (j=0; j<x; j++) idx[nn++] = xbase++; /* interior */
736: if (n5 >= 0) { /* directly right */
737: x_t = lx[n5 % m];
738: /* y_t = y; */
739: s_t = bases[n5] + (i)*x_t;
740: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
741: } else if (Xe-xe > 0) {
742: if (bx == DM_BOUNDARY_MIRROR) {
743: for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*(i + 1) - 2 - j;
744: } else {
745: for (j=0; j<s_x; j++) idx[nn++] = -1;
746: }
747: }
748: }
750: for (i=1; i<=s_y; i++) {
751: if (n6 >= 0) { /* left above */
752: x_t = lx[n6 % m];
753: /* y_t = ly[(n6/m)]; */
754: s_t = bases[n6] + (i)*x_t - s_x;
755: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
756: } else if (xs-Xs > 0 && Ye-ye > 0) {
757: for (j=0; j<s_x; j++) idx[nn++] = -1;
758: }
759: if (n7 >= 0) { /* directly above */
760: x_t = x;
761: /* y_t = ly[(n7/m)]; */
762: s_t = bases[n7] + (i-1)*x_t;
763: for (j=0; j<x_t; j++) idx[nn++] = s_t++;
764: } else if (Ye-ye > 0) {
765: if (by == DM_BOUNDARY_MIRROR) {
766: for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(y - i - 1) + j;
767: } else {
768: for (j=0; j<x; j++) idx[nn++] = -1;
769: }
770: }
771: if (n8 >= 0) { /* right above */
772: x_t = lx[n8 % m];
773: /* y_t = ly[(n8/m)]; */
774: s_t = bases[n8] + (i-1)*x_t;
775: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
776: } else if (Xe-xe > 0 && Ye-ye > 0) {
777: for (j=0; j<s_x; j++) idx[nn++] = -1;
778: }
779: }
780: }
781: /*
782: Set the local to global ordering in the global vector, this allows use
783: of VecSetValuesLocal().
784: */
785: ISLocalToGlobalMappingCreate(comm,dof,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);
786: PetscLogObjectParent((PetscObject)da,(PetscObject)da->ltogmap);
788: PetscFree2(bases,ldims);
789: dd->m = m; dd->n = n;
790: /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
791: dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1;
792: dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1;
794: VecDestroy(&local);
795: VecDestroy(&global);
797: dd->gtol = gtol;
798: dd->base = base;
799: da->ops->view = DMView_DA_2d;
800: dd->ltol = NULL;
801: dd->ao = NULL;
802: return(0);
803: }
805: /*@C
806: DMDACreate2d - Creates an object that will manage the communication of two-dimensional
807: regular array data that is distributed across some processors.
809: Collective on MPI_Comm
811: Input Parameters:
812: + comm - MPI communicator
813: . bx,by - type of ghost nodes the array have.
814: Use one of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC.
815: . stencil_type - stencil type. Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
816: . M,N - global dimension in each direction of the array
817: . m,n - corresponding number of processors in each dimension
818: (or PETSC_DECIDE to have calculated)
819: . dof - number of degrees of freedom per node
820: . s - stencil width
821: - lx, ly - arrays containing the number of nodes in each cell along
822: the x and y coordinates, or NULL. If non-null, these
823: must be of length as m and n, and the corresponding
824: m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
825: must be M, and the sum of the ly[] entries must be N.
827: Output Parameter:
828: . da - the resulting distributed array object
830: Options Database Key:
831: + -dm_view - Calls DMView() at the conclusion of DMDACreate2d()
832: . -da_grid_x <nx> - number of grid points in x direction, if M < 0
833: . -da_grid_y <ny> - number of grid points in y direction, if N < 0
834: . -da_processors_x <nx> - number of processors in x direction
835: . -da_processors_y <ny> - number of processors in y direction
836: . -da_refine_x <rx> - refinement ratio in x direction
837: . -da_refine_y <ry> - refinement ratio in y direction
838: - -da_refine <n> - refine the DMDA n times before creating, if M or N < 0
841: Level: beginner
843: Notes:
844: The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the
845: standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes
846: the standard 9-pt stencil.
848: The array data itself is NOT stored in the DMDA, it is stored in Vec objects;
849: The appropriate vector objects can be obtained with calls to DMCreateGlobalVector()
850: and DMCreateLocalVector() and calls to VecDuplicate() if more are needed.
852: You must call DMSetUp() after this call before using this DM.
854: If you wish to use the options database to change values in the DMDA call DMSetFromOptions() after this call
855: but before DMSetUp().
857: .keywords: distributed array, create, two-dimensional
859: .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(),
860: DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDASetRefinementFactor(),
861: DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges()
863: @*/
865: PetscErrorCode DMDACreate2d(MPI_Comm comm,DMBoundaryType bx,DMBoundaryType by,DMDAStencilType stencil_type,
866: PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da)
867: {
871: DMDACreate(comm, da);
872: DMSetDimension(*da, 2);
873: DMDASetSizes(*da, M, N, 1);
874: DMDASetNumProcs(*da, m, n, PETSC_DECIDE);
875: DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE);
876: DMDASetDof(*da, dof);
877: DMDASetStencilType(*da, stencil_type);
878: DMDASetStencilWidth(*da, s);
879: DMDASetOwnershipRanges(*da, lx, ly, NULL);
880: return(0);
881: }