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