Actual source code: da2.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/
3: #include <petscdraw.h>
7: 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: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
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: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
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,*idx;
47: char node[10];
48: PetscBool isnull;
50: PetscViewerDrawGetDraw(viewer,0,&draw);
51: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
52: if (!da->coordinates) {
53: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
54: }
55: PetscDrawSynchronizedClear(draw);
57: /* first processor draw all node lines */
58: if (!rank) {
59: ymin = 0.0; ymax = dd->N - 1;
60: for (xmin=0; xmin<dd->M; xmin++) {
61: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
62: }
63: xmin = 0.0; xmax = dd->M - 1;
64: for (ymin=0; ymin<dd->N; ymin++) {
65: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
66: }
67: }
68: PetscDrawSynchronizedFlush(draw);
69: PetscDrawPause(draw);
71: /* draw my box */
72: ymin = dd->ys; ymax = dd->ye - 1; xmin = dd->xs/dd->w;
73: xmax =(dd->xe-1)/dd->w;
74: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
75: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
76: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
77: PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);
79: /* put in numbers */
80: base = (dd->base)/dd->w;
81: for (y=ymin; y<=ymax; y++) {
82: for (x=xmin; x<=xmax; x++) {
83: sprintf(node,"%d",(int)base++);
84: PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
85: }
86: }
88: PetscDrawSynchronizedFlush(draw);
89: PetscDrawPause(draw);
90: /* overlay ghost numbers, useful for error checking */
91: /* put in numbers */
93: base = 0; idx = dd->idx;
94: ymin = dd->Ys; ymax = dd->Ye; xmin = dd->Xs; xmax = dd->Xe;
95: for (y=ymin; y<ymax; y++) {
96: for (x=xmin; x<xmax; x++) {
97: if ((base % dd->w) == 0) {
98: sprintf(node,"%d",(int)(idx[base]/dd->w));
99: PetscDrawString(draw,x/dd->w,y,PETSC_DRAW_BLUE,node);
100: }
101: base++;
102: }
103: }
104: PetscDrawSynchronizedFlush(draw);
105: PetscDrawPause(draw);
106: } else if (isbinary) {
107: DMView_DA_Binary(da,viewer);
108: #if defined(PETSC_HAVE_MATLAB_ENGINE)
109: } else if (ismatlab) {
110: DMView_DA_Matlab(da,viewer);
111: #endif
112: }
113: return(0);
114: }
116: /*
117: M is number of grid points
118: m is number of processors
120: */
123: PetscErrorCode DMDASplitComm2d(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt sw,MPI_Comm *outcomm)
124: {
126: PetscInt m,n = 0,x = 0,y = 0;
127: PetscMPIInt size,csize,rank;
130: MPI_Comm_size(comm,&size);
131: MPI_Comm_rank(comm,&rank);
133: csize = 4*size;
134: do {
135: 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);
136: csize = csize/4;
138: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)csize)/((PetscReal)N)));
139: if (!m) m = 1;
140: while (m > 0) {
141: n = csize/m;
142: if (m*n == csize) break;
143: m--;
144: }
145: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
147: x = M/m + ((M % m) > ((csize-1) % m));
148: y = (N + (csize-1)/m)/n;
149: } while ((x < 4 || y < 4) && csize > 1);
150: if (size != csize) {
151: MPI_Group entire_group,sub_group;
152: PetscMPIInt i,*groupies;
154: MPI_Comm_group(comm,&entire_group);
155: PetscMalloc(csize*sizeof(PetscInt),&groupies);
156: for (i=0; i<csize; i++) {
157: groupies[i] = (rank/csize)*csize + i;
158: }
159: MPI_Group_incl(entire_group,csize,groupies,&sub_group);
160: PetscFree(groupies);
161: MPI_Comm_create(comm,sub_group,outcomm);
162: MPI_Group_free(&entire_group);
163: MPI_Group_free(&sub_group);
164: PetscInfo1(0,"DMDASplitComm2d:Creating redundant coarse problems of size %d\n",csize);
165: } else {
166: *outcomm = comm;
167: }
168: return(0);
169: }
171: #if defined(new)
174: /*
175: DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local
176: function lives on a DMDA
178: y ~= (F(u + ha) - F(u))/h,
179: where F = nonlinear function, as set by SNESSetFunction()
180: u = current iterate
181: h = difference interval
182: */
183: PetscErrorCode DMDAGetDiagonal_MFFD(DM da,Vec U,Vec a)
184: {
185: PetscScalar h,*aa,*ww,v;
186: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
188: PetscInt gI,nI;
189: MatStencil stencil;
190: DMDALocalInfo info;
193: (*ctx->func)(0,U,a,ctx->funcctx);
194: (*ctx->funcisetbase)(U,ctx->funcctx);
196: VecGetArray(U,&ww);
197: VecGetArray(a,&aa);
199: nI = 0;
200: h = ww[gI];
201: if (h == 0.0) h = 1.0;
202: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
203: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
204: h *= epsilon;
206: ww[gI] += h;
207: (*ctx->funci)(i,w,&v,ctx->funcctx);
208: aa[nI] = (v - aa[nI])/h;
209: ww[gI] -= h;
210: nI++;
212: VecRestoreArray(U,&ww);
213: VecRestoreArray(a,&aa);
214: return(0);
215: }
216: #endif
220: PetscErrorCode DMSetUp_DA_2D(DM da)
221: {
222: DM_DA *dd = (DM_DA*)da->data;
223: const PetscInt M = dd->M;
224: const PetscInt N = dd->N;
225: PetscInt m = dd->m;
226: PetscInt n = dd->n;
227: const PetscInt dof = dd->w;
228: const PetscInt s = dd->s;
229: DMDABoundaryType bx = dd->bx;
230: DMDABoundaryType by = dd->by;
231: DMDAStencilType stencil_type = dd->stencil_type;
232: PetscInt *lx = dd->lx;
233: PetscInt *ly = dd->ly;
234: MPI_Comm comm;
235: PetscMPIInt rank,size;
236: PetscInt xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end,IXs,IXe,IYs,IYe;
237: PetscInt up,down,left,right,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn,*idx_cpy;
238: const PetscInt *idx_full;
239: PetscInt xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count;
240: PetscInt s_x,s_y; /* s proportionalized to w */
241: PetscInt sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0;
242: Vec local,global;
243: VecScatter ltog,gtol;
244: IS to,from,ltogis;
245: PetscErrorCode ierr;
248: if (stencil_type == DMDA_STENCIL_BOX && (bx == DMDA_BOUNDARY_MIRROR || by == DMDA_BOUNDARY_MIRROR)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Mirror boundary and box stencil");
249: PetscObjectGetComm((PetscObject)da,&comm);
250: #if !defined(PETSC_USE_64BIT_INDICES)
251: 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);
252: #endif
254: if (dof < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
255: if (s < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);
257: MPI_Comm_size(comm,&size);
258: MPI_Comm_rank(comm,&rank);
260: if (m != PETSC_DECIDE) {
261: if (m < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
262: else if (m > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
263: }
264: if (n != PETSC_DECIDE) {
265: if (n < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
266: else if (n > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
267: }
269: if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
270: if (n != PETSC_DECIDE) {
271: m = size/n;
272: } else if (m != PETSC_DECIDE) {
273: n = size/m;
274: } else {
275: /* try for squarish distribution */
276: m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M)*((PetscReal)size)/((PetscReal)N)));
277: if (!m) m = 1;
278: while (m > 0) {
279: n = size/m;
280: if (m*n == size) break;
281: m--;
282: }
283: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
284: }
285: if (m*n != size) SETERRQ(comm,PETSC_ERR_PLIB,"Unable to create partition, check the size of the communicator and input m and n ");
286: } else if (m*n != size) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");
288: if (M < m) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
289: if (N < n) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);
291: /*
292: Determine locally owned region
293: xs is the first local node number, x is the number of local nodes
294: */
295: if (!lx) {
296: PetscMalloc(m*sizeof(PetscInt), &dd->lx);
297: lx = dd->lx;
298: for (i=0; i<m; i++) {
299: lx[i] = M/m + ((M % m) > i);
300: }
301: }
302: x = lx[rank % m];
303: xs = 0;
304: for (i=0; i<(rank % m); i++) {
305: xs += lx[i];
306: }
307: #if defined(PETSC_USE_DEBUG)
308: left = xs;
309: for (i=(rank % m); i<m; i++) {
310: left += lx[i];
311: }
312: if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M);
313: #endif
315: /*
316: Determine locally owned region
317: ys is the first local node number, y is the number of local nodes
318: */
319: if (!ly) {
320: PetscMalloc(n*sizeof(PetscInt), &dd->ly);
321: ly = dd->ly;
322: for (i=0; i<n; i++) {
323: ly[i] = N/n + ((N % n) > i);
324: }
325: }
326: y = ly[rank/m];
327: ys = 0;
328: for (i=0; i<(rank/m); i++) {
329: ys += ly[i];
330: }
331: #if defined(PETSC_USE_DEBUG)
332: left = ys;
333: for (i=(rank/m); i<n; i++) {
334: left += ly[i];
335: }
336: if (left != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N);
337: #endif
339: /*
340: check if the scatter requires more than one process neighbor or wraps around
341: the domain more than once
342: */
343: if ((x < s) && ((m > 1) || (bx == DMDA_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);
344: if ((y < s) && ((n > 1) || (by == DMDA_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);
345: xe = xs + x;
346: ye = ys + y;
348: /* determine ghost region (Xs) and region scattered into (IXs) */
349: if (xs-s > 0) {
350: Xs = xs - s; IXs = xs - s;
351: } else {
352: if (bx) {
353: Xs = xs - s;
354: } else {
355: Xs = 0;
356: }
357: IXs = 0;
358: }
359: if (xe+s <= M) {
360: Xe = xe + s; IXe = xe + s;
361: } else {
362: if (bx) {
363: Xs = xs - s; Xe = xe + s;
364: } else {
365: Xe = M;
366: }
367: IXe = M;
368: }
370: if (bx == DMDA_BOUNDARY_PERIODIC || bx == DMDA_BOUNDARY_MIRROR) {
371: IXs = xs - s;
372: IXe = xe + s;
373: Xs = xs - s;
374: Xe = xe + s;
375: }
377: if (ys-s > 0) {
378: Ys = ys - s; IYs = ys - s;
379: } else {
380: if (by) {
381: Ys = ys - s;
382: } else {
383: Ys = 0;
384: }
385: IYs = 0;
386: }
387: if (ye+s <= N) {
388: Ye = ye + s; IYe = ye + s;
389: } else {
390: if (by) {
391: Ye = ye + s;
392: } else {
393: Ye = N;
394: }
395: IYe = N;
396: }
398: if (by == DMDA_BOUNDARY_PERIODIC || by == DMDA_BOUNDARY_MIRROR) {
399: IYs = ys - s;
400: IYe = ye + s;
401: Ys = ys - s;
402: Ye = ye + s;
403: }
405: /* stencil length in each direction */
406: s_x = s;
407: s_y = s;
409: /* determine starting point of each processor */
410: nn = x*y;
411: PetscMalloc2(size+1,PetscInt,&bases,size,PetscInt,&ldims);
412: MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
413: bases[0] = 0;
414: for (i=1; i<=size; i++) {
415: bases[i] = ldims[i-1];
416: }
417: for (i=1; i<=size; i++) {
418: bases[i] += bases[i-1];
419: }
420: base = bases[rank]*dof;
422: /* allocate the base parallel and sequential vectors */
423: dd->Nlocal = x*y*dof;
424: VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,0,&global);
425: dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof;
426: VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,0,&local);
428: /* generate appropriate vector scatters */
429: /* local to global inserts non-ghost point region into global */
430: VecGetOwnershipRange(global,&start,&end);
431: ISCreateStride(comm,x*y*dof,start,1,&to);
433: PetscMalloc(x*y*sizeof(PetscInt),&idx);
434: left = xs - Xs; right = left + x;
435: down = ys - Ys; up = down + y;
436: count = 0;
437: for (i=down; i<up; i++) {
438: for (j=left; j<right; j++) {
439: idx[count++] = i*(Xe-Xs) + j;
440: }
441: }
443: ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);
444: VecScatterCreate(local,from,global,to,<og);
445: PetscLogObjectParent(dd,ltog);
446: ISDestroy(&from);
447: ISDestroy(&to);
449: /* global to local must include ghost points within the domain,
450: but not ghost points outside the domain that aren't periodic */
451: if (stencil_type == DMDA_STENCIL_BOX) {
452: count = (IXe-IXs)*(IYe-IYs);
453: PetscMalloc(count*sizeof(PetscInt),&idx);
455: left = IXs - Xs; right = left + (IXe-IXs);
456: down = IYs - Ys; up = down + (IYe-IYs);
457: count = 0;
458: for (i=down; i<up; i++) {
459: for (j=left; j<right; j++) {
460: idx[count++] = j + i*(Xe-Xs);
461: }
462: }
463: ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
465: } else {
466: /* must drop into cross shape region */
467: /* ---------|
468: | top |
469: |--- ---| up
470: | middle |
471: | |
472: ---- ---- down
473: | bottom |
474: -----------
475: Xs xs xe Xe */
476: count = (ys-IYs)*x + y*(IXe-IXs) + (IYe-ye)*x;
477: PetscMalloc(count*sizeof(PetscInt),&idx);
479: left = xs - Xs; right = left + x;
480: down = ys - Ys; up = down + y;
481: count = 0;
482: /* bottom */
483: for (i=(IYs-Ys); i<down; i++) {
484: for (j=left; j<right; j++) {
485: idx[count++] = j + i*(Xe-Xs);
486: }
487: }
488: /* middle */
489: for (i=down; i<up; i++) {
490: for (j=(IXs-Xs); j<(IXe-Xs); j++) {
491: idx[count++] = j + i*(Xe-Xs);
492: }
493: }
494: /* top */
495: for (i=up; i<up+IYe-ye; i++) {
496: for (j=left; j<right; j++) {
497: idx[count++] = j + i*(Xe-Xs);
498: }
499: }
500: ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
501: }
504: /* determine who lies on each side of us stored in n6 n7 n8
505: n3 n5
506: n0 n1 n2
507: */
509: /* Assume the Non-Periodic Case */
510: n1 = rank - m;
511: if (rank % m) {
512: n0 = n1 - 1;
513: } else {
514: n0 = -1;
515: }
516: if ((rank+1) % m) {
517: n2 = n1 + 1;
518: n5 = rank + 1;
519: n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
520: } else {
521: n2 = -1; n5 = -1; n8 = -1;
522: }
523: if (rank % m) {
524: n3 = rank - 1;
525: n6 = n3 + m; if (n6 >= m*n) n6 = -1;
526: } else {
527: n3 = -1; n6 = -1;
528: }
529: n7 = rank + m; if (n7 >= m*n) n7 = -1;
531: if (bx == DMDA_BOUNDARY_PERIODIC && by == DMDA_BOUNDARY_PERIODIC) {
532: /* Modify for Periodic Cases */
533: /* Handle all four corners */
534: if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
535: if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
536: if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
537: if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;
539: /* Handle Top and Bottom Sides */
540: if (n1 < 0) n1 = rank + m * (n-1);
541: if (n7 < 0) n7 = rank - m * (n-1);
542: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
543: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
544: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
545: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
547: /* Handle Left and Right Sides */
548: if (n3 < 0) n3 = rank + (m-1);
549: if (n5 < 0) n5 = rank - (m-1);
550: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
551: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
552: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
553: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
554: } else if (by == DMDA_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */
555: if (n1 < 0) n1 = rank + m * (n-1);
556: if (n7 < 0) n7 = rank - m * (n-1);
557: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
558: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
559: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
560: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
561: } else if (bx == DMDA_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */
562: if (n3 < 0) n3 = rank + (m-1);
563: if (n5 < 0) n5 = rank - (m-1);
564: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
565: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
566: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
567: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
568: }
570: PetscMalloc(9*sizeof(PetscInt),&dd->neighbors);
572: dd->neighbors[0] = n0;
573: dd->neighbors[1] = n1;
574: dd->neighbors[2] = n2;
575: dd->neighbors[3] = n3;
576: dd->neighbors[4] = rank;
577: dd->neighbors[5] = n5;
578: dd->neighbors[6] = n6;
579: dd->neighbors[7] = n7;
580: dd->neighbors[8] = n8;
582: if (stencil_type == DMDA_STENCIL_STAR) {
583: /* save corner processor numbers */
584: sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8;
585: n0 = n2 = n6 = n8 = -1;
586: }
588: PetscMalloc((Xe-Xs)*(Ye-Ys)*sizeof(PetscInt),&idx);
589: PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*sizeof(PetscInt));
591: nn = 0;
592: xbase = bases[rank];
593: for (i=1; i<=s_y; i++) {
594: if (n0 >= 0) { /* left below */
595: x_t = lx[n0 % m];
596: y_t = ly[(n0/m)];
597: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
598: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
599: }
601: if (n1 >= 0) { /* directly below */
602: x_t = x;
603: y_t = ly[(n1/m)];
604: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
605: for (j=0; j<x_t; j++) idx[nn++] = s_t++;
606: } else if (by == DMDA_BOUNDARY_MIRROR) {
607: for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(s_y - i + 1) + j;
608: }
610: if (n2 >= 0) { /* right below */
611: x_t = lx[n2 % m];
612: y_t = ly[(n2/m)];
613: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
614: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
615: }
616: }
618: for (i=0; i<y; i++) {
619: if (n3 >= 0) { /* directly left */
620: x_t = lx[n3 % m];
621: /* y_t = y; */
622: s_t = bases[n3] + (i+1)*x_t - s_x;
623: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
624: } else if (bx == DMDA_BOUNDARY_MIRROR) {
625: for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*i + s_x - j;
626: }
628: for (j=0; j<x; j++) idx[nn++] = xbase++; /* interior */
630: if (n5 >= 0) { /* directly right */
631: x_t = lx[n5 % m];
632: /* y_t = y; */
633: s_t = bases[n5] + (i)*x_t;
634: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
635: } else if (bx == DMDA_BOUNDARY_MIRROR) {
636: for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*(i + 1) - 2 - j;
637: }
638: }
640: for (i=1; i<=s_y; i++) {
641: if (n6 >= 0) { /* left above */
642: x_t = lx[n6 % m];
643: /* y_t = ly[(n6/m)]; */
644: s_t = bases[n6] + (i)*x_t - s_x;
645: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
646: }
648: if (n7 >= 0) { /* directly above */
649: x_t = x;
650: /* y_t = ly[(n7/m)]; */
651: s_t = bases[n7] + (i-1)*x_t;
652: for (j=0; j<x_t; j++) idx[nn++] = s_t++;
653: } else if (by == DMDA_BOUNDARY_MIRROR) {
654: for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(y - i - 1) + j;
655: }
657: if (n8 >= 0) { /* right above */
658: x_t = lx[n8 % m];
659: /* y_t = ly[(n8/m)]; */
660: s_t = bases[n8] + (i-1)*x_t;
661: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
662: }
663: }
665: ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);
666: VecScatterCreate(global,from,local,to,>ol);
667: PetscLogObjectParent(da,gtol);
668: ISDestroy(&to);
669: ISDestroy(&from);
671: if (stencil_type == DMDA_STENCIL_STAR) {
672: n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
673: }
675: if (((stencil_type == DMDA_STENCIL_STAR) ||
676: (bx && bx != DMDA_BOUNDARY_PERIODIC) ||
677: (by && by != DMDA_BOUNDARY_PERIODIC))) {
678: /*
679: Recompute the local to global mappings, this time keeping the
680: information about the cross corner processor numbers and any ghosted
681: but not periodic indices.
682: */
683: nn = 0;
684: xbase = bases[rank];
685: for (i=1; i<=s_y; i++) {
686: if (n0 >= 0) { /* left below */
687: x_t = lx[n0 % m];
688: y_t = ly[(n0/m)];
689: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
690: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
691: } else if (xs-Xs > 0 && ys-Ys > 0) {
692: for (j=0; j<s_x; j++) idx[nn++] = -1;
693: }
694: if (n1 >= 0) { /* directly below */
695: x_t = x;
696: y_t = ly[(n1/m)];
697: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
698: for (j=0; j<x_t; j++) idx[nn++] = s_t++;
699: } else if (ys-Ys > 0) {
700: if (by == DMDA_BOUNDARY_MIRROR) {
701: for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(s_y - i + 1) + j;
702: } else {
703: for (j=0; j<x; j++) idx[nn++] = -1;
704: }
705: }
706: if (n2 >= 0) { /* right below */
707: x_t = lx[n2 % m];
708: y_t = ly[(n2/m)];
709: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
710: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
711: } else if (Xe-xe> 0 && ys-Ys > 0) {
712: for (j=0; j<s_x; j++) idx[nn++] = -1;
713: }
714: }
716: for (i=0; i<y; i++) {
717: if (n3 >= 0) { /* directly left */
718: x_t = lx[n3 % m];
719: /* y_t = y; */
720: s_t = bases[n3] + (i+1)*x_t - s_x;
721: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
722: } else if (xs-Xs > 0) {
723: if (bx == DMDA_BOUNDARY_MIRROR) {
724: for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*i + s_x - j;
725: } else {
726: for (j=0; j<s_x; j++) idx[nn++] = -1;
727: }
728: }
730: for (j=0; j<x; j++) idx[nn++] = xbase++; /* interior */
732: if (n5 >= 0) { /* directly right */
733: x_t = lx[n5 % m];
734: /* y_t = y; */
735: s_t = bases[n5] + (i)*x_t;
736: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
737: } else if (Xe-xe > 0) {
738: if (bx == DMDA_BOUNDARY_MIRROR) {
739: for (j=0; j<s_x; j++) idx[nn++] = bases[rank] + x*(i + 1) - 2 - j;
740: } else {
741: for (j=0; j<s_x; j++) idx[nn++] = -1;
742: }
743: }
744: }
746: for (i=1; i<=s_y; i++) {
747: if (n6 >= 0) { /* left above */
748: x_t = lx[n6 % m];
749: /* y_t = ly[(n6/m)]; */
750: s_t = bases[n6] + (i)*x_t - s_x;
751: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
752: } else if (xs-Xs > 0 && Ye-ye > 0) {
753: for (j=0; j<s_x; j++) idx[nn++] = -1;
754: }
755: if (n7 >= 0) { /* directly above */
756: x_t = x;
757: /* y_t = ly[(n7/m)]; */
758: s_t = bases[n7] + (i-1)*x_t;
759: for (j=0; j<x_t; j++) idx[nn++] = s_t++;
760: } else if (Ye-ye > 0) {
761: if (by == DMDA_BOUNDARY_MIRROR) {
762: for (j=0; j<x; j++) idx[nn++] = bases[rank] + x*(y - i - 1) + j;
763: } else {
764: for (j=0; j<x; j++) idx[nn++] = -1;
765: }
766: }
767: if (n8 >= 0) { /* right above */
768: x_t = lx[n8 % m];
769: /* y_t = ly[(n8/m)]; */
770: s_t = bases[n8] + (i-1)*x_t;
771: for (j=0; j<s_x; j++) idx[nn++] = s_t++;
772: } else if (Xe-xe > 0 && Ye-ye > 0) {
773: for (j=0; j<s_x; j++) idx[nn++] = -1;
774: }
775: }
776: }
777: /*
778: Set the local to global ordering in the global vector, this allows use
779: of VecSetValuesLocal().
780: */
781: ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);
782: PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);
783: PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));
784: ISGetIndices(ltogis, &idx_full);
785: PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));
786: ISRestoreIndices(ltogis, &idx_full);
787: ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);
788: PetscLogObjectParent(da,da->ltogmap);
789: ISDestroy(<ogis);
790: ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);
791: PetscLogObjectParent(da,da->ltogmap);
793: PetscFree2(bases,ldims);
794: dd->m = m; dd->n = n;
795: /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
796: dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1;
797: dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1;
799: VecDestroy(&local);
800: VecDestroy(&global);
802: dd->gtol = gtol;
803: dd->ltog = ltog;
804: dd->idx = idx_cpy;
805: dd->Nl = nn*dof;
806: dd->base = base;
807: da->ops->view = DMView_DA_2d;
808: dd->ltol = NULL;
809: dd->ao = NULL;
810: return(0);
811: }
815: /*@C
816: DMDACreate2d - Creates an object that will manage the communication of two-dimensional
817: regular array data that is distributed across some processors.
819: Collective on MPI_Comm
821: Input Parameters:
822: + comm - MPI communicator
823: . bx,by - type of ghost nodes the array have.
824: Use one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC.
825: . stencil_type - stencil type. Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
826: . 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
827: from the command line with -da_grid_x <M> -da_grid_y <N>)
828: . m,n - corresponding number of processors in each dimension
829: (or PETSC_DECIDE to have calculated)
830: . dof - number of degrees of freedom per node
831: . s - stencil width
832: - lx, ly - arrays containing the number of nodes in each cell along
833: the x and y coordinates, or NULL. If non-null, these
834: must be of length as m and n, and the corresponding
835: m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
836: must be M, and the sum of the ly[] entries must be N.
838: Output Parameter:
839: . da - the resulting distributed array object
841: Options Database Key:
842: + -dm_view - Calls DMView() at the conclusion of DMDACreate2d()
843: . -da_grid_x <nx> - number of grid points in x direction, if M < 0
844: . -da_grid_y <ny> - number of grid points in y direction, if N < 0
845: . -da_processors_x <nx> - number of processors in x direction
846: . -da_processors_y <ny> - number of processors in y direction
847: . -da_refine_x <rx> - refinement ratio in x direction
848: . -da_refine_y <ry> - refinement ratio in y direction
849: - -da_refine <n> - refine the DMDA n times before creating, if M or N < 0
852: Level: beginner
854: Notes:
855: The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the
856: standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes
857: the standard 9-pt stencil.
859: The array data itself is NOT stored in the DMDA, it is stored in Vec objects;
860: The appropriate vector objects can be obtained with calls to DMCreateGlobalVector()
861: and DMCreateLocalVector() and calls to VecDuplicate() if more are needed.
863: .keywords: distributed array, create, two-dimensional
865: .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(),
866: DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDASetRefinementFactor(),
867: DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges()
869: @*/
871: PetscErrorCode DMDACreate2d(MPI_Comm comm,DMDABoundaryType bx,DMDABoundaryType by,DMDAStencilType stencil_type,
872: PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da)
873: {
877: DMDACreate(comm, da);
878: DMDASetDim(*da, 2);
879: DMDASetSizes(*da, M, N, 1);
880: DMDASetNumProcs(*da, m, n, PETSC_DECIDE);
881: DMDASetBoundaryType(*da, bx, by, DMDA_BOUNDARY_NONE);
882: DMDASetDof(*da, dof);
883: DMDASetStencilType(*da, stencil_type);
884: DMDASetStencilWidth(*da, s);
885: DMDASetOwnershipRanges(*da, lx, ly, NULL);
886: /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */
887: DMSetFromOptions(*da);
888: DMSetUp(*da);
889: DMViewFromOptions(*da,NULL,"-dm_view");
890: return(0);
891: }