Actual source code: da2.c
1: #define PETSCDM_DLL
2:
3: #include src/dm/da/daimpl.h
7: /*@C
8: DAGetElements - Gets an array containing the indices (in local coordinates)
9: of all the local elements
11: Not Collective
13: Input Parameter:
14: . da - the DA object
16: Output Parameters:
17: + n - number of local elements
18: - e - the indices of the elements vertices
20: Level: intermediate
22: .seealso: DAElementType, DASetElementType(), DARestoreElements()
23: @*/
24: PetscErrorCode DAGetElements(DA da,PetscInt *n,const PetscInt *e[])
25: {
29: (da->ops->getelements)(da,n,e);
30: return(0);
31: }
35: /*@C
36: DARestoreElements - Returns an array containing the indices (in local coordinates)
37: of all the local elements obtained with DAGetElements()
39: Not Collective
41: Input Parameter:
42: + da - the DA object
43: . n - number of local elements
44: - e - the indices of the elements vertices
46: Level: intermediate
48: .seealso: DAElementType, DASetElementType(), DAGetElements()
49: @*/
50: PetscErrorCode DARestoreElements(DA da,PetscInt *n,const PetscInt *e[])
51: {
55: if (da->ops->restoreelements) {
56: (da->ops->restoreelements)(da,n,e);
57: }
58: return(0);
59: }
63: PetscErrorCode DAGetOwnershipRange(DA da,PetscInt **lx,PetscInt **ly,PetscInt **lz)
64: {
67: if (lx) *lx = da->lx;
68: if (ly) *ly = da->ly;
69: if (lz) *lz = da->lz;
70: return(0);
71: }
75: PetscErrorCode DAView_2d(DA da,PetscViewer viewer)
76: {
78: PetscMPIInt rank;
79: PetscTruth iascii,isdraw;
82: MPI_Comm_rank(da->comm,&rank);
84: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
85: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
86: if (iascii) {
87: PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D m %D n %D w %D s %D\n",rank,da->M,
88: da->N,da->m,da->n,da->w,da->s);
89: PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D\n",da->xs,da->xe,da->ys,da->ye);
90: PetscViewerFlush(viewer);
91: } else if (isdraw) {
92: PetscDraw draw;
93: double ymin = -1*da->s-1,ymax = da->N+da->s;
94: double xmin = -1*da->s-1,xmax = da->M+da->s;
95: double x,y;
96: PetscInt base,*idx;
97: char node[10];
98: PetscTruth isnull;
99:
100: PetscViewerDrawGetDraw(viewer,0,&draw);
101: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
102: if (!da->coordinates) {
103: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
104: }
105: PetscDrawSynchronizedClear(draw);
107: /* first processor draw all node lines */
108: if (!rank) {
109: ymin = 0.0; ymax = da->N - 1;
110: for (xmin=0; xmin<da->M; xmin++) {
111: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
112: }
113: xmin = 0.0; xmax = da->M - 1;
114: for (ymin=0; ymin<da->N; ymin++) {
115: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
116: }
117: }
118: PetscDrawSynchronizedFlush(draw);
119: PetscDrawPause(draw);
121: /* draw my box */
122: ymin = da->ys; ymax = da->ye - 1; xmin = da->xs/da->w;
123: xmax =(da->xe-1)/da->w;
124: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
125: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
126: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
127: PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);
129: /* put in numbers */
130: base = (da->base)/da->w;
131: for (y=ymin; y<=ymax; y++) {
132: for (x=xmin; x<=xmax; x++) {
133: sprintf(node,"%d",(int)base++);
134: PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
135: }
136: }
138: PetscDrawSynchronizedFlush(draw);
139: PetscDrawPause(draw);
140: /* overlay ghost numbers, useful for error checking */
141: /* put in numbers */
143: base = 0; idx = da->idx;
144: ymin = da->Ys; ymax = da->Ye; xmin = da->Xs; xmax = da->Xe;
145: for (y=ymin; y<ymax; y++) {
146: for (x=xmin; x<xmax; x++) {
147: if ((base % da->w) == 0) {
148: sprintf(node,"%d",(int)(idx[base]/da->w));
149: PetscDrawString(draw,x/da->w,y,PETSC_DRAW_BLUE,node);
150: }
151: base++;
152: }
153: }
154: PetscDrawSynchronizedFlush(draw);
155: PetscDrawPause(draw);
156: } else {
157: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for DA2d",((PetscObject)viewer)->type_name);
158: }
159: return(0);
160: }
164: PetscErrorCode DAPublish_Petsc(PetscObject obj)
165: {
167: return(0);
168: }
172: PetscErrorCode DAGetElements_2d_P1(DA da,PetscInt *n,const PetscInt *e[])
173: {
175: PetscInt i,j,cnt,xs,xe = da->xe,ys,ye = da->ye,Xs = da->Xs, Xe = da->Xe, Ys = da->Ys;
178: if (!da->e) {
179: if (da->xs == Xs) xs = da->xs; else xs = da->xs - 1;
180: if (da->ys == Ys) ys = da->ys; else ys = da->ys - 1;
181: da->ne = 2*(xe - xs - 1)*(ye - ys - 1);
182: PetscMalloc((1 + 3*da->ne)*sizeof(PetscInt),&da->e);
183: cnt = 0;
184: for (j=ys; j<ye-1; j++) {
185: for (i=xs; i<xe-1; i++) {
186: da->e[cnt] = i - Xs + (j - Ys)*(Xe - Xs);
187: da->e[cnt+1] = i - Xs + 1 + (j - Ys)*(Xe - Xs);
188: da->e[cnt+2] = i - Xs + (j - Ys + 1)*(Xe - Xs);
190: da->e[cnt+3] = i - Xs + 1 + (j - Ys + 1)*(Xe - Xs);
191: da->e[cnt+4] = i - Xs + (j - Ys + 1)*(Xe - Xs);
192: da->e[cnt+5] = i - Xs + 1 + (j - Ys)*(Xe - Xs);
193: cnt += 6;
194: }
195: }
196: }
197: *n = da->ne;
198: *e = da->e;
199: return(0);
200: }
205: /*@C
206: DACreate2d - Creates an object that will manage the communication of two-dimensional
207: regular array data that is distributed across some processors.
209: Collective on MPI_Comm
211: Input Parameters:
212: + comm - MPI communicator
213: . wrap - type of periodicity should the array have.
214: Use one of DA_NONPERIODIC, DA_XPERIODIC, DA_YPERIODIC, or DA_XYPERIODIC.
215: . stencil_type - stencil type. Use either DA_STENCIL_BOX or DA_STENCIL_STAR.
216: . 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
217: from the command line with -da_grid_x <M> -da_grid_y <N>)
218: . m,n - corresponding number of processors in each dimension
219: (or PETSC_DECIDE to have calculated)
220: . dof - number of degrees of freedom per node
221: . s - stencil width
222: - lx, ly - arrays containing the number of nodes in each cell along
223: the x and y coordinates, or PETSC_NULL. If non-null, these
224: must be of length as m and n, and the corresponding
225: m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
226: must be M, and the sum of the ly[] entries must be N.
228: Output Parameter:
229: . inra - the resulting distributed array object
231: Options Database Key:
232: + -da_view - Calls DAView() at the conclusion of DACreate2d()
233: . -da_grid_x <nx> - number of grid points in x direction, if M < 0
234: . -da_grid_y <ny> - number of grid points in y direction, if N < 0
235: . -da_processors_x <nx> - number of processors in x direction
236: . -da_processors_y <ny> - number of processors in y direction
237: . -da_refine_x - refinement ratio in x direction
238: - -da_refine_y - refinement ratio in y direction
240: Level: beginner
242: Notes:
243: The stencil type DA_STENCIL_STAR with width 1 corresponds to the
244: standard 5-pt stencil, while DA_STENCIL_BOX with width 1 denotes
245: the standard 9-pt stencil.
247: The array data itself is NOT stored in the DA, it is stored in Vec objects;
248: The appropriate vector objects can be obtained with calls to DACreateGlobalVector()
249: and DACreateLocalVector() and calls to VecDuplicate() if more are needed.
251: .keywords: distributed array, create, two-dimensional
253: .seealso: DADestroy(), DAView(), DACreate1d(), DACreate3d(), DAGlobalToLocalBegin(), DAGetRefinementFactor(),
254: DAGlobalToLocalEnd(), DALocalToGlobal(), DALocalToLocalBegin(), DALocalToLocalEnd(), DASetRefinementFactor(),
255: DAGetInfo(), DACreateGlobalVector(), DACreateLocalVector(), DACreateNaturalVector(), DALoad(), DAView()
257: @*/
258: PetscErrorCode DACreate2d(MPI_Comm comm,DAPeriodicType wrap,DAStencilType stencil_type,
259: PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,PetscInt *lx,PetscInt *ly,DA *inra)
260: {
262: PetscMPIInt rank,size;
263: PetscInt xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end;
264: PetscInt up,down,left,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn;
265: PetscInt xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count;
266: PetscInt s_x,s_y; /* s proportionalized to w */
267: PetscInt *flx = 0,*fly = 0;
268: PetscInt sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0,refine_x = 2, refine_y = 2,tM = M,tN = N;
269: DA da;
270: Vec local,global;
271: VecScatter ltog,gtol;
272: IS to,from;
276: *inra = 0;
277: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
278: DMInitializePackage(PETSC_NULL);
279: #endif
281: if (dof < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
282: if (s < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);
284: PetscOptionsBegin(comm,PETSC_NULL,"2d DA Options","DA");
285: if (M < 0){
286: tM = -M;
287: PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DACreate2d",tM,&tM,PETSC_NULL);
288: }
289: if (N < 0){
290: tN = -N;
291: PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DACreate2d",tN,&tN,PETSC_NULL);
292: }
293: PetscOptionsInt("-da_processors_x","Number of processors in x direction","DACreate2d",m,&m,PETSC_NULL);
294: PetscOptionsInt("-da_processors_y","Number of processors in y direction","DACreate2d",n,&n,PETSC_NULL);
295: PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DASetRefinementFactor",refine_x,&refine_x,PETSC_NULL);
296: PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DASetRefinementFactor",refine_y,&refine_y,PETSC_NULL);
297: PetscOptionsEnd();
298: M = tM; N = tN;
300: PetscHeaderCreate(da,_p_DA,struct _DAOps,DA_COOKIE,0,"DA",comm,DADestroy,DAView);
301: da->bops->publish = DAPublish_Petsc;
302: da->ops->createglobalvector = DACreateGlobalVector;
303: da->ops->getinterpolation = DAGetInterpolation;
304: da->ops->getcoloring = DAGetColoring;
305: da->ops->getmatrix = DAGetMatrix;
306: da->ops->refine = DARefine;
307: da->ops->getinjection = DAGetInjection;
308: da->ops->getelements = DAGetElements_2d_P1;
309: da->elementtype = DA_ELEMENT_P1;
311: PetscLogObjectMemory(da,sizeof(struct _p_DA));
312: da->dim = 2;
313: da->interptype = DA_Q1;
314: da->refine_x = refine_x;
315: da->refine_y = refine_y;
316: PetscMalloc(dof*sizeof(char*),&da->fieldname);
317: PetscMemzero(da->fieldname,dof*sizeof(char*));
319: MPI_Comm_size(comm,&size);
320: MPI_Comm_rank(comm,&rank);
322: if (m != PETSC_DECIDE) {
323: if (m < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);}
324: else if (m > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);}
325: }
326: if (n != PETSC_DECIDE) {
327: if (n < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);}
328: else if (n > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);}
329: }
331: if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
332: if (n != PETSC_DECIDE) {
333: m = size/n;
334: } else if (m != PETSC_DECIDE) {
335: n = size/m;
336: } else {
337: /* try for squarish distribution */
338: m = (PetscInt)(0.5 + sqrt(((double)M)*((double)size)/((double)N)));
339: if (!m) m = 1;
340: while (m > 0) {
341: n = size/m;
342: if (m*n == size) break;
343: m--;
344: }
345: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
346: }
347: if (m*n != size) SETERRQ(PETSC_ERR_PLIB,"Unable to create partition, check the size of the communicator and input m and n ");
348: } else if (m*n != size) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");
350: if (M < m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
351: if (N < n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);
353: /*
354: Determine locally owned region
355: xs is the first local node number, x is the number of local nodes
356: */
357: if (!lx) { /* user sets distribution */
358: PetscMalloc(m*sizeof(PetscInt),&lx);
359: flx = lx;
360: for (i=0; i<m; i++) {
361: lx[i] = M/m + ((M % m) > i);
362: }
363: }
364: x = lx[rank % m];
365: xs = 0;
366: for (i=0; i<(rank % m); i++) {
367: xs += lx[i];
368: }
369: #if defined(PETSC_USE_DEBUG)
370: left = xs;
371: for (i=(rank % m); i<m; i++) {
372: left += lx[i];
373: }
374: if (left != M) {
375: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M);
376: }
377: #endif
379: /*
380: Determine locally owned region
381: ys is the first local node number, y is the number of local nodes
382: */
383: if (!ly) { /* user sets distribution */
384: PetscMalloc(n*sizeof(PetscInt),&ly);
385: fly = ly;
386: for (i=0; i<n; i++) {
387: ly[i] = N/n + ((N % n) > i);
388: }
389: }
390: y = ly[rank/m];
391: ys = 0;
392: for (i=0; i<(rank/m); i++) {
393: ys += ly[i];
394: }
395: #if defined(PETSC_USE_DEBUG)
396: left = ys;
397: for (i=(rank/m); i<n; i++) {
398: left += ly[i];
399: }
400: if (left != N) {
401: SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N);
402: }
403: #endif
405: if (x < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Local x-width of domain x %D is smaller than stencil width s %D",x,s);
406: if (y < s) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Local y-width of domain y %D is smaller than stencil width s %D",y,s);
407: xe = xs + x;
408: ye = ys + y;
410: /* determine ghost region */
411: /* Assume No Periodicity */
412: if (xs-s > 0) Xs = xs - s; else Xs = 0;
413: if (ys-s > 0) Ys = ys - s; else Ys = 0;
414: if (xe+s <= M) Xe = xe + s; else Xe = M;
415: if (ye+s <= N) Ye = ye + s; else Ye = N;
417: /* X Periodic */
418: if (DAXPeriodic(wrap)){
419: Xs = xs - s;
420: Xe = xe + s;
421: }
423: /* Y Periodic */
424: if (DAYPeriodic(wrap)){
425: Ys = ys - s;
426: Ye = ye + s;
427: }
429: /* Resize all X parameters to reflect w */
430: x *= dof;
431: xs *= dof;
432: xe *= dof;
433: Xs *= dof;
434: Xe *= dof;
435: s_x = s*dof;
436: s_y = s;
438: /* determine starting point of each processor */
439: nn = x*y;
440: PetscMalloc((2*size+1)*sizeof(PetscInt),&bases);
441: ldims = bases+size+1;
442: MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
443: bases[0] = 0;
444: for (i=1; i<=size; i++) {
445: bases[i] = ldims[i-1];
446: }
447: for (i=1; i<=size; i++) {
448: bases[i] += bases[i-1];
449: }
451: /* allocate the base parallel and sequential vectors */
452: da->Nlocal = x*y;
453: VecCreateMPIWithArray(comm,da->Nlocal,PETSC_DECIDE,0,&global);
454: VecSetBlockSize(global,dof);
455: da->nlocal = (Xe-Xs)*(Ye-Ys);
456: VecCreateSeqWithArray(PETSC_COMM_SELF,da->nlocal,0,&local);
457: VecSetBlockSize(local,dof);
460: /* generate appropriate vector scatters */
461: /* local to global inserts non-ghost point region into global */
462: VecGetOwnershipRange(global,&start,&end);
463: ISCreateStride(comm,x*y,start,1,&to);
465: left = xs - Xs; down = ys - Ys; up = down + y;
466: PetscMalloc(x*(up - down)*sizeof(PetscInt),&idx);
467: count = 0;
468: for (i=down; i<up; i++) {
469: for (j=0; j<x/dof; j++) {
470: idx[count++] = left + i*(Xe-Xs) + j*dof;
471: }
472: }
473: ISCreateBlock(comm,dof,count,idx,&from);
474: PetscFree(idx);
476: VecScatterCreate(local,from,global,to,<og);
477: PetscLogObjectParent(da,to);
478: PetscLogObjectParent(da,from);
479: PetscLogObjectParent(da,ltog);
480: ISDestroy(from);
481: ISDestroy(to);
483: /* global to local must include ghost points */
484: if (stencil_type == DA_STENCIL_BOX) {
485: ISCreateStride(comm,(Xe-Xs)*(Ye-Ys),0,1,&to);
486: } else {
487: /* must drop into cross shape region */
488: /* ---------|
489: | top |
490: |--- ---|
491: | middle |
492: | |
493: ---- ----
494: | bottom |
495: -----------
496: Xs xs xe Xe */
497: /* bottom */
498: left = xs - Xs; down = ys - Ys; up = down + y;
499: count = down*(xe-xs) + (up-down)*(Xe-Xs) + (Ye-Ys-up)*(xe-xs);
500: PetscMalloc(count*sizeof(PetscInt)/dof,&idx);
501: count = 0;
502: for (i=0; i<down; i++) {
503: for (j=0; j<xe-xs; j += dof) {
504: idx[count++] = left + i*(Xe-Xs) + j;
505: }
506: }
507: /* middle */
508: for (i=down; i<up; i++) {
509: for (j=0; j<Xe-Xs; j += dof) {
510: idx[count++] = i*(Xe-Xs) + j;
511: }
512: }
513: /* top */
514: for (i=up; i<Ye-Ys; i++) {
515: for (j=0; j<xe-xs; j += dof) {
516: idx[count++] = left + i*(Xe-Xs) + j;
517: }
518: }
519: ISCreateBlock(comm,dof,count,idx,&to);
520: PetscFree(idx);
521: }
524: /* determine who lies on each side of us stored in n6 n7 n8
525: n3 n5
526: n0 n1 n2
527: */
529: /* Assume the Non-Periodic Case */
530: n1 = rank - m;
531: if (rank % m) {
532: n0 = n1 - 1;
533: } else {
534: n0 = -1;
535: }
536: if ((rank+1) % m) {
537: n2 = n1 + 1;
538: n5 = rank + 1;
539: n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
540: } else {
541: n2 = -1; n5 = -1; n8 = -1;
542: }
543: if (rank % m) {
544: n3 = rank - 1;
545: n6 = n3 + m; if (n6 >= m*n) n6 = -1;
546: } else {
547: n3 = -1; n6 = -1;
548: }
549: n7 = rank + m; if (n7 >= m*n) n7 = -1;
552: /* Modify for Periodic Cases */
553: if (wrap == DA_YPERIODIC) { /* Handle Top and Bottom Sides */
554: if (n1 < 0) n1 = rank + m * (n-1);
555: if (n7 < 0) n7 = rank - m * (n-1);
556: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
557: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
558: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
559: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
560: } else if (wrap == DA_XPERIODIC) { /* Handle Left and Right Sides */
561: if (n3 < 0) n3 = rank + (m-1);
562: if (n5 < 0) n5 = rank - (m-1);
563: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
564: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
565: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
566: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
567: } else if (wrap == DA_XYPERIODIC) {
569: /* Handle all four corners */
570: if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
571: if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
572: if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
573: if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;
575: /* Handle Top and Bottom Sides */
576: if (n1 < 0) n1 = rank + m * (n-1);
577: if (n7 < 0) n7 = rank - m * (n-1);
578: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
579: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
580: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
581: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
583: /* Handle Left and Right Sides */
584: if (n3 < 0) n3 = rank + (m-1);
585: if (n5 < 0) n5 = rank - (m-1);
586: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
587: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
588: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
589: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
590: }
592: if (stencil_type == DA_STENCIL_STAR) {
593: /* save corner processor numbers */
594: sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8;
595: n0 = n2 = n6 = n8 = -1;
596: }
598: PetscMalloc((x+2*s_x)*(y+2*s_y)*sizeof(PetscInt),&idx);
599: PetscLogObjectMemory(da,(x+2*s_x)*(y+2*s_y)*sizeof(PetscInt));
600: nn = 0;
602: xbase = bases[rank];
603: for (i=1; i<=s_y; i++) {
604: if (n0 >= 0) { /* left below */
605: x_t = lx[n0 % m]*dof;
606: y_t = ly[(n0/m)];
607: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
608: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
609: }
610: if (n1 >= 0) { /* directly below */
611: x_t = x;
612: y_t = ly[(n1/m)];
613: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
614: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
615: }
616: if (n2 >= 0) { /* right below */
617: x_t = lx[n2 % m]*dof;
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]*dof;
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: }
632: for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */
634: if (n5 >= 0) { /* directly right */
635: x_t = lx[n5 % m]*dof;
636: /* y_t = y; */
637: s_t = bases[n5] + (i)*x_t;
638: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
639: }
640: }
642: for (i=1; i<=s_y; i++) {
643: if (n6 >= 0) { /* left above */
644: x_t = lx[n6 % m]*dof;
645: /* y_t = ly[(n6/m)]; */
646: s_t = bases[n6] + (i)*x_t - s_x;
647: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
648: }
649: if (n7 >= 0) { /* directly above */
650: x_t = x;
651: /* y_t = ly[(n7/m)]; */
652: s_t = bases[n7] + (i-1)*x_t;
653: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
654: }
655: if (n8 >= 0) { /* right above */
656: x_t = lx[n8 % m]*dof;
657: /* y_t = ly[(n8/m)]; */
658: s_t = bases[n8] + (i-1)*x_t;
659: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
660: }
661: }
663: base = bases[rank];
664: {
665: PetscInt nnn = nn/dof,*iidx;
666: PetscMalloc(nnn*sizeof(PetscInt),&iidx);
667: for (i=0; i<nnn; i++) {
668: iidx[i] = idx[dof*i];
669: }
670: ISCreateBlock(comm,dof,nnn,iidx,&from);
671: PetscFree(iidx);
672: }
673: VecScatterCreate(global,from,local,to,>ol);
674: PetscLogObjectParent(da,to);
675: PetscLogObjectParent(da,from);
676: PetscLogObjectParent(da,gtol);
677: ISDestroy(to);
678: ISDestroy(from);
680: if (stencil_type == DA_STENCIL_STAR) {
681: /*
682: Recompute the local to global mappings, this time keeping the
683: information about the cross corner processor numbers.
684: */
685: n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
686: nn = 0;
687: xbase = bases[rank];
688: for (i=1; i<=s_y; i++) {
689: if (n0 >= 0) { /* left below */
690: x_t = lx[n0 % m]*dof;
691: y_t = ly[(n0/m)];
692: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
693: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
694: }
695: if (n1 >= 0) { /* directly below */
696: x_t = x;
697: y_t = ly[(n1/m)];
698: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
699: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
700: }
701: if (n2 >= 0) { /* right below */
702: x_t = lx[n2 % m]*dof;
703: y_t = ly[(n2/m)];
704: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
705: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
706: }
707: }
709: for (i=0; i<y; i++) {
710: if (n3 >= 0) { /* directly left */
711: x_t = lx[n3 % m]*dof;
712: /* y_t = y; */
713: s_t = bases[n3] + (i+1)*x_t - s_x;
714: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
715: }
717: for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */
719: if (n5 >= 0) { /* directly right */
720: x_t = lx[n5 % m]*dof;
721: /* y_t = y; */
722: s_t = bases[n5] + (i)*x_t;
723: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
724: }
725: }
727: for (i=1; i<=s_y; i++) {
728: if (n6 >= 0) { /* left above */
729: x_t = lx[n6 % m]*dof;
730: /* y_t = ly[(n6/m)]; */
731: s_t = bases[n6] + (i)*x_t - s_x;
732: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
733: }
734: if (n7 >= 0) { /* directly above */
735: x_t = x;
736: /* y_t = ly[(n7/m)]; */
737: s_t = bases[n7] + (i-1)*x_t;
738: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
739: }
740: if (n8 >= 0) { /* right above */
741: x_t = lx[n8 % m]*dof;
742: /* y_t = ly[(n8/m)]; */
743: s_t = bases[n8] + (i-1)*x_t;
744: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
745: }
746: }
747: }
748: PetscFree(bases);
750: da->M = M; da->N = N; da->m = m; da->n = n; da->w = dof; da->s = s;
751: da->xs = xs; da->xe = xe; da->ys = ys; da->ye = ye; da->zs = 0; da->ze = 1;
752: da->Xs = Xs; da->Xe = Xe; da->Ys = Ys; da->Ye = Ye; da->Zs = 0; da->Ze = 1;
753: da->P = 1; da->p = 1;
755: VecDestroy(local);
756: VecDestroy(global);
758: da->gtol = gtol;
759: da->ltog = ltog;
760: da->idx = idx;
761: da->Nl = nn;
762: da->base = base;
763: da->wrap = wrap;
764: da->ops->view = DAView_2d;
765: da->stencil_type = stencil_type;
767: /*
768: Set the local to global ordering in the global vector, this allows use
769: of VecSetValuesLocal().
770: */
771: ISLocalToGlobalMappingCreateNC(comm,nn,idx,&da->ltogmap);
772: ISLocalToGlobalMappingBlock(da->ltogmap,da->w,&da->ltogmapb);
773: PetscLogObjectParent(da,da->ltogmap);
775: *inra = da;
777: da->ltol = PETSC_NULL;
778: da->ao = PETSC_NULL;
781: if (!flx) {
782: PetscMalloc(m*sizeof(PetscInt),&flx);
783: PetscMemcpy(flx,lx,m*sizeof(PetscInt));
784: }
785: if (!fly) {
786: PetscMalloc(n*sizeof(PetscInt),&fly);
787: PetscMemcpy(fly,ly,n*sizeof(PetscInt));
788: }
789: da->lx = flx;
790: da->ly = fly;
791: DAView_Private(da);
792: PetscPublishAll(da);
793: return(0);
794: }
798: /*@
799: DARefine - Creates a new distributed array that is a refinement of a given
800: distributed array.
802: Collective on DA
804: Input Parameter:
805: + da - initial distributed array
806: - comm - communicator to contain refined DA, must be either same as the da communicator or include the
807: da communicator and be 2, 4, or 8 times larger. Currently ignored
809: Output Parameter:
810: . daref - refined distributed array
812: Level: advanced
814: Note:
815: Currently, refinement consists of just doubling the number of grid spaces
816: in each dimension of the DA.
818: .keywords: distributed array, refine
820: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy()
821: @*/
822: PetscErrorCode DARefine(DA da,MPI_Comm comm,DA *daref)
823: {
825: PetscInt M,N,P;
826: DA da2;
832: if (DAXPeriodic(da->wrap) || da->interptype == DA_Q0){
833: M = da->refine_x*da->M;
834: } else {
835: M = 1 + da->refine_x*(da->M - 1);
836: }
837: if (DAYPeriodic(da->wrap) || da->interptype == DA_Q0){
838: N = da->refine_y*da->N;
839: } else {
840: N = 1 + da->refine_y*(da->N - 1);
841: }
842: if (DAZPeriodic(da->wrap) || da->interptype == DA_Q0){
843: P = da->refine_z*da->P;
844: } else {
845: P = 1 + da->refine_z*(da->P - 1);
846: }
847: DACreate(da->comm,da->dim,da->wrap,da->stencil_type,M,N,P,da->m,da->n,da->p,da->w,da->s,0,0,0,&da2);
849: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
850: da2->ops->getmatrix = da->ops->getmatrix;
851: da2->ops->getinterpolation = da->ops->getinterpolation;
852: da2->ops->getcoloring = da->ops->getcoloring;
853: da2->interptype = da->interptype;
854:
855: /* copy fill information if given */
856: if (da->dfill) {
857: PetscMalloc((da->dfill[da->w]+da->w+1)*sizeof(PetscInt),&da2->dfill);
858: PetscMemcpy(da2->dfill,da->dfill,(da->dfill[da->w]+da->w+1)*sizeof(PetscInt));
859: }
860: if (da->ofill) {
861: PetscMalloc((da->ofill[da->w]+da->w+1)*sizeof(PetscInt),&da2->ofill);
862: PetscMemcpy(da2->ofill,da->ofill,(da->ofill[da->w]+da->w+1)*sizeof(PetscInt));
863: }
864: /* copy the refine information */
865: da2->refine_x = da->refine_x;
866: da2->refine_y = da->refine_y;
867: da2->refine_z = da->refine_z;
868: *daref = da2;
869: return(0);
870: }
872: /*@C
873: DASetRefinementFactor - Set the ratios that the DA grid is refined
875: Collective on DA
877: Input Parameters:
878: + da - the DA object
879: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
880: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
881: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
883: Options Database:
884: + -da_refine_x - refinement ratio in x direction
885: . -da_refine_y - refinement ratio in y direction
886: - -da_refine_y - refinement ratio in z direction
888: Level: intermediate
890: Notes: Pass PETSC_IGNORE to leave a value unchanged
892: .seealso: DARefine(), DAGetRefinementFactor()
893: @*/
894: PetscErrorCode DASetRefinementFactor(DA da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
895: {
897: if (refine_x > 0) da->refine_x = refine_x;
898: if (refine_y > 0) da->refine_y = refine_y;
899: if (refine_z > 0) da->refine_z = refine_z;
900: return(0);
901: }
903: /*@C
904: DAGetRefinementFactor - Gets the ratios that the DA grid is refined
906: Not Collective
908: Input Parameter:
909: . da - the DA object
911: Output Parameters:
912: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
913: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
914: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
916: Level: intermediate
918: Notes: Pass PETSC_NULL for values you do not need
920: .seealso: DARefine(), DASetRefinementFactor()
921: @*/
922: PetscErrorCode DAGetRefinementFactor(DA da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
923: {
925: if (refine_x) *refine_x = da->refine_x;
926: if (refine_y) *refine_y = da->refine_y;
927: if (refine_z) *refine_z = da->refine_z;
928: return(0);
929: }
931: /*@C
932: DASetGetMatrix - Sets the routine used by the DA to allocate a matrix.
934: Collective on DA
936: Input Parameters:
937: + da - the DA object
938: - f - the function that allocates the matrix for that specific DA
940: Level: developer
942: Notes: See DASetBlockFills() that provides a simple way to provide the nonzero structure for
943: the diagonal and off-diagonal blocks of the matrix
945: .seealso: DAGetMatrix(), DASetBlockFills()
946: @*/
947: PetscErrorCode DASetGetMatrix(DA da,PetscErrorCode (*f)(DA, MatType,Mat*))
948: {
950: da->ops->getmatrix = f;
951: return(0);
952: }
954: /*
955: M is number of grid points
956: m is number of processors
958: */
961: PetscErrorCode DASplitComm2d(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt sw,MPI_Comm *outcomm)
962: {
964: PetscInt m,n = 0,x = 0,y = 0;
965: PetscMPIInt size,csize,rank;
968: MPI_Comm_size(comm,&size);
969: MPI_Comm_rank(comm,&rank);
971: csize = 4*size;
972: do {
973: if (csize % 4) SETERRQ4(PETSC_ERR_ARG_INCOMP,"Cannot split communicator of size %d tried %d %D %D",size,csize,x,y);
974: csize = csize/4;
975:
976: m = (PetscInt)(0.5 + sqrt(((double)M)*((double)csize)/((double)N)));
977: if (!m) m = 1;
978: while (m > 0) {
979: n = csize/m;
980: if (m*n == csize) break;
981: m--;
982: }
983: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
985: x = M/m + ((M % m) > ((csize-1) % m));
986: y = (N + (csize-1)/m)/n;
987: } while ((x < 4 || y < 4) && csize > 1);
988: if (size != csize) {
989: MPI_Group entire_group,sub_group;
990: PetscMPIInt i,*groupies;
992: MPI_Comm_group(comm,&entire_group);
993: PetscMalloc(csize*sizeof(PetscInt),&groupies);
994: for (i=0; i<csize; i++) {
995: groupies[i] = (rank/csize)*csize + i;
996: }
997: MPI_Group_incl(entire_group,csize,groupies,&sub_group);
998: PetscFree(groupies);
999: MPI_Comm_create(comm,sub_group,outcomm);
1000: MPI_Group_free(&entire_group);
1001: MPI_Group_free(&sub_group);
1002: PetscInfo1(0,"DASplitComm2d:Creating redundant coarse problems of size %d\n",csize);
1003: } else {
1004: *outcomm = comm;
1005: }
1006: return(0);
1007: }
1011: /*@C
1012: DASetLocalFunction - Caches in a DA a local function.
1014: Collective on DA
1016: Input Parameter:
1017: + da - initial distributed array
1018: - lf - the local function
1020: Level: intermediate
1022: Notes: The routine SNESDAFormFunction() uses this the cached function to evaluate the user provided function.
1024: .keywords: distributed array, refine
1026: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunctioni()
1027: @*/
1028: PetscErrorCode DASetLocalFunction(DA da,DALocalFunction1 lf)
1029: {
1032: da->lf = lf;
1033: return(0);
1034: }
1038: /*@C
1039: DASetLocalFunctioni - Caches in a DA a local function that evaluates a single component
1041: Collective on DA
1043: Input Parameter:
1044: + da - initial distributed array
1045: - lfi - the local function
1047: Level: intermediate
1049: .keywords: distributed array, refine
1051: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction()
1052: @*/
1053: PetscErrorCode DASetLocalFunctioni(DA da,PetscErrorCode (*lfi)(DALocalInfo*,MatStencil*,void*,PetscScalar*,void*))
1054: {
1057: da->lfi = lfi;
1058: return(0);
1059: }
1063: /*@C
1064: DASetLocalFunctionib - Caches in a DA a block local function that evaluates a single component
1066: Collective on DA
1068: Input Parameter:
1069: + da - initial distributed array
1070: - lfi - the local function
1072: Level: intermediate
1074: .keywords: distributed array, refine
1076: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction()
1077: @*/
1078: PetscErrorCode DASetLocalFunctionib(DA da,PetscErrorCode (*lfi)(DALocalInfo*,MatStencil*,void*,PetscScalar*,void*))
1079: {
1082: da->lfib = lfi;
1083: return(0);
1084: }
1088: PetscErrorCode DASetLocalAdicFunction_Private(DA da,DALocalFunction1 ad_lf)
1089: {
1092: da->adic_lf = ad_lf;
1093: return(0);
1094: }
1096: /*MC
1097: DASetLocalAdicFunctioni - Caches in a DA a local functioni computed by ADIC/ADIFOR
1099: Collective on DA
1101: Synopsis:
1102: PetscErrorCode DASetLocalAdicFunctioni(DA da,PetscInt (ad_lf*)(DALocalInfo*,MatStencil*,void*,void*,void*)
1103:
1104: Input Parameter:
1105: + da - initial distributed array
1106: - ad_lfi - the local function as computed by ADIC/ADIFOR
1108: Level: intermediate
1110: .keywords: distributed array, refine
1112: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1113: DASetLocalJacobian(), DASetLocalFunctioni()
1114: M*/
1118: PetscErrorCode DASetLocalAdicFunctioni_Private(DA da,PetscErrorCode (*ad_lfi)(DALocalInfo*,MatStencil*,void*,void*,void*))
1119: {
1122: da->adic_lfi = ad_lfi;
1123: return(0);
1124: }
1126: /*MC
1127: DASetLocalAdicMFFunctioni - Caches in a DA a local functioni computed by ADIC/ADIFOR
1129: Collective on DA
1131: Synopsis:
1132: PetscErrorCode DASetLocalAdicFunctioni(DA da,int (ad_lf*)(DALocalInfo*,MatStencil*,void*,void*,void*)
1133:
1134: Input Parameter:
1135: + da - initial distributed array
1136: - admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR
1138: Level: intermediate
1140: .keywords: distributed array, refine
1142: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1143: DASetLocalJacobian(), DASetLocalFunctioni()
1144: M*/
1148: PetscErrorCode DASetLocalAdicMFFunctioni_Private(DA da,PetscErrorCode (*admf_lfi)(DALocalInfo*,MatStencil*,void*,void*,void*))
1149: {
1152: da->adicmf_lfi = admf_lfi;
1153: return(0);
1154: }
1156: /*MC
1157: DASetLocalAdicFunctionib - Caches in a DA a block local functioni computed by ADIC/ADIFOR
1159: Collective on DA
1161: Synopsis:
1162: PetscErrorCode DASetLocalAdicFunctionib(DA da,PetscInt (ad_lf*)(DALocalInfo*,MatStencil*,void*,void*,void*)
1163:
1164: Input Parameter:
1165: + da - initial distributed array
1166: - ad_lfi - the local function as computed by ADIC/ADIFOR
1168: Level: intermediate
1170: .keywords: distributed array, refine
1172: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1173: DASetLocalJacobian(), DASetLocalFunctionib()
1174: M*/
1178: PetscErrorCode DASetLocalAdicFunctionib_Private(DA da,PetscErrorCode (*ad_lfi)(DALocalInfo*,MatStencil*,void*,void*,void*))
1179: {
1182: da->adic_lfib = ad_lfi;
1183: return(0);
1184: }
1186: /*MC
1187: DASetLocalAdicMFFunctionib - Caches in a DA a block local functioni computed by ADIC/ADIFOR
1189: Collective on DA
1191: Synopsis:
1192: PetscErrorCode DASetLocalAdicFunctionib(DA da,int (ad_lf*)(DALocalInfo*,MatStencil*,void*,void*,void*)
1193:
1194: Input Parameter:
1195: + da - initial distributed array
1196: - admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR
1198: Level: intermediate
1200: .keywords: distributed array, refine
1202: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1203: DASetLocalJacobian(), DASetLocalFunctionib()
1204: M*/
1208: PetscErrorCode DASetLocalAdicMFFunctionib_Private(DA da,PetscErrorCode (*admf_lfi)(DALocalInfo*,MatStencil*,void*,void*,void*))
1209: {
1212: da->adicmf_lfib = admf_lfi;
1213: return(0);
1214: }
1216: /*MC
1217: DASetLocalAdicMFFunction - Caches in a DA a local function computed by ADIC/ADIFOR
1219: Collective on DA
1221: Synopsis:
1222: PetscErrorCode DASetLocalAdicMFFunction(DA da,DALocalFunction1 ad_lf)
1223:
1224: Input Parameter:
1225: + da - initial distributed array
1226: - ad_lf - the local function as computed by ADIC/ADIFOR
1228: Level: intermediate
1230: .keywords: distributed array, refine
1232: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
1233: DASetLocalJacobian()
1234: M*/
1238: PetscErrorCode DASetLocalAdicMFFunction_Private(DA da,DALocalFunction1 ad_lf)
1239: {
1242: da->adicmf_lf = ad_lf;
1243: return(0);
1244: }
1246: /*@C
1247: DASetLocalJacobian - Caches in a DA a local Jacobian
1249: Collective on DA
1251:
1252: Input Parameter:
1253: + da - initial distributed array
1254: - lj - the local Jacobian
1256: Level: intermediate
1258: Notes: The routine SNESDAFormFunction() uses this the cached function to evaluate the user provided function.
1260: .keywords: distributed array, refine
1262: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction()
1263: @*/
1266: PetscErrorCode DASetLocalJacobian(DA da,DALocalFunction1 lj)
1267: {
1270: da->lj = lj;
1271: return(0);
1272: }
1276: /*@C
1277: DAGetLocalFunction - Gets from a DA a local function and its ADIC/ADIFOR Jacobian
1279: Collective on DA
1281: Input Parameter:
1282: . da - initial distributed array
1284: Output Parameters:
1285: . lf - the local function
1287: Level: intermediate
1289: .keywords: distributed array, refine
1291: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DASetLocalFunction()
1292: @*/
1293: PetscErrorCode