Actual source code: da1.c
petsc-3.7.5 2017-01-01
2: /*
3: Code for manipulating distributed regular 1d arrays in parallel.
4: This file was created by Peter Mell 6/30/95
5: */
7: #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
9: #include <petscdraw.h>
12: static PetscErrorCode DMView_DA_1d(DM da,PetscViewer viewer)
13: {
15: PetscMPIInt rank;
16: PetscBool iascii,isdraw,isbinary;
17: DM_DA *dd = (DM_DA*)da->data;
18: #if defined(PETSC_HAVE_MATLAB_ENGINE)
19: PetscBool ismatlab;
20: #endif
23: MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);
25: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
26: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
27: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
28: #if defined(PETSC_HAVE_MATLAB_ENGINE)
29: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
30: #endif
31: if (iascii) {
32: PetscViewerFormat format;
34: PetscViewerGetFormat(viewer, &format);
35: if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL) {
36: DMDALocalInfo info;
37: DMDAGetLocalInfo(da,&info);
38: PetscViewerASCIIPushSynchronized(viewer);
39: PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D m %D w %D s %D\n",rank,dd->M,dd->m,dd->w,dd->s);
40: PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D\n",info.xs,info.xs+info.xm);
41: PetscViewerFlush(viewer);
42: PetscViewerASCIIPopSynchronized(viewer);
43: } else {
44: DMView_DA_VTK(da, viewer);
45: }
46: } else if (isdraw) {
47: PetscDraw draw;
48: double ymin = -1,ymax = 1,xmin = -1,xmax = dd->M,x;
49: PetscInt base;
50: char node[10];
51: PetscBool isnull;
53: PetscViewerDrawGetDraw(viewer,0,&draw);
54: PetscDrawIsNull(draw,&isnull);
55: if (isnull) return(0);
57: PetscDrawCheckResizedWindow(draw);
58: PetscDrawClear(draw);
59: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
61: PetscDrawCollectiveBegin(draw);
62: /* first processor draws all node lines */
63: if (!rank) {
64: PetscInt xmin_tmp;
65: ymin = 0.0; ymax = 0.3;
66: for (xmin_tmp=0; xmin_tmp < dd->M; xmin_tmp++) {
67: PetscDrawLine(draw,(double)xmin_tmp,ymin,(double)xmin_tmp,ymax,PETSC_DRAW_BLACK);
68: }
69: xmin = 0.0; xmax = dd->M - 1;
70: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
71: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_BLACK);
72: }
73: PetscDrawCollectiveEnd(draw);
74: PetscDrawFlush(draw);
75: PetscDrawPause(draw);
77: PetscDrawCollectiveBegin(draw);
78: /* draw my box */
79: ymin = 0; ymax = 0.3; xmin = dd->xs / dd->w; xmax = (dd->xe / dd->w) - 1;
80: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
81: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
82: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
83: PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);
84: /* Put in index numbers */
85: base = dd->base / dd->w;
86: for (x=xmin; x<=xmax; x++) {
87: PetscSNPrintf(node,sizeof(node),"%d",(int)base++);
88: PetscDrawString(draw,x,ymin,PETSC_DRAW_RED,node);
89: }
90: PetscDrawCollectiveEnd(draw);
91: PetscDrawFlush(draw);
92: PetscDrawPause(draw);
93: PetscDrawSave(draw);
94: } else if (isbinary) {
95: DMView_DA_Binary(da,viewer);
96: #if defined(PETSC_HAVE_MATLAB_ENGINE)
97: } else if (ismatlab) {
98: DMView_DA_Matlab(da,viewer);
99: #endif
100: }
101: return(0);
102: }
107: PetscErrorCode DMSetUp_DA_1D(DM da)
108: {
109: DM_DA *dd = (DM_DA*)da->data;
110: const PetscInt M = dd->M;
111: const PetscInt dof = dd->w;
112: const PetscInt s = dd->s;
113: const PetscInt sDist = s; /* stencil distance in points */
114: const PetscInt *lx = dd->lx;
115: DMBoundaryType bx = dd->bx;
116: MPI_Comm comm;
117: Vec local, global;
118: VecScatter gtol;
119: IS to, from;
120: PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE;
121: PetscMPIInt rank, size;
122: PetscInt i,*idx,nn,left,xs,xe,x,Xs,Xe,start,m,IXs,IXe;
123: PetscErrorCode ierr;
126: PetscObjectGetComm((PetscObject) da, &comm);
127: MPI_Comm_size(comm,&size);
128: MPI_Comm_rank(comm,&rank);
130: dd->p = 1;
131: dd->n = 1;
132: dd->m = size;
133: m = dd->m;
135: if (s > 0) {
136: /* if not communicating data then should be ok to have nothing on some processes */
137: if (M < m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"More processes than data points! %D %D",m,M);
138: if ((M-1) < s && size > 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Array is too small for stencil! %D %D",M-1,s);
139: }
141: /*
142: Determine locally owned region
143: xs is the first local node number, x is the number of local nodes
144: */
145: if (!lx) {
146: PetscMalloc1(m, &dd->lx);
147: PetscOptionsGetBool(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_partition_blockcomm",&flg1,NULL);
148: PetscOptionsGetBool(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_partition_nodes_at_end",&flg2,NULL);
149: if (flg1) { /* Block Comm type Distribution */
150: xs = rank*M/m;
151: x = (rank + 1)*M/m - xs;
152: } else if (flg2) { /* The odd nodes are evenly distributed across last nodes */
153: x = (M + rank)/m;
154: if (M/m == x) xs = rank*x;
155: else xs = rank*(x-1) + (M+rank)%(x*m);
156: } else { /* The odd nodes are evenly distributed across the first k nodes */
157: /* Regular PETSc Distribution */
158: x = M/m + ((M % m) > rank);
159: if (rank >= (M % m)) xs = (rank * (PetscInt)(M/m) + M % m);
160: else xs = rank * (PetscInt)(M/m) + rank;
161: }
162: MPI_Allgather(&xs,1,MPIU_INT,dd->lx,1,MPIU_INT,comm);
163: for (i=0; i<m-1; i++) dd->lx[i] = dd->lx[i+1] - dd->lx[i];
164: dd->lx[m-1] = M - dd->lx[m-1];
165: } else {
166: x = lx[rank];
167: xs = 0;
168: for (i=0; i<rank; i++) xs += lx[i];
169: /* verify that data user provided is consistent */
170: left = xs;
171: for (i=rank; i<size; i++) left += lx[i];
172: if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M %D %D",left,M);
173: }
175: /*
176: check if the scatter requires more than one process neighbor or wraps around
177: the domain more than once
178: */
179: 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);
181: xe = xs + x;
183: /* determine ghost region (Xs) and region scattered into (IXs) */
184: if (xs-sDist > 0) {
185: Xs = xs - sDist;
186: IXs = xs - sDist;
187: } else {
188: if (bx) Xs = xs - sDist;
189: else Xs = 0;
190: IXs = 0;
191: }
192: if (xe+sDist <= M) {
193: Xe = xe + sDist;
194: IXe = xe + sDist;
195: } else {
196: if (bx) Xe = xe + sDist;
197: else Xe = M;
198: IXe = M;
199: }
201: if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) {
202: Xs = xs - sDist;
203: Xe = xe + sDist;
204: IXs = xs - sDist;
205: IXe = xe + sDist;
206: }
208: /* allocate the base parallel and sequential vectors */
209: dd->Nlocal = dof*x;
210: VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,NULL,&global);
211: dd->nlocal = dof*(Xe-Xs);
212: VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,NULL,&local);
214: VecGetOwnershipRange(global,&start,NULL);
216: /* Create Global to Local Vector Scatter Context */
217: /* global to local must retrieve ghost points */
218: ISCreateStride(comm,dof*(IXe-IXs),dof*(IXs-Xs),1,&to);
220: PetscMalloc1(x+2*sDist,&idx);
221: PetscLogObjectMemory((PetscObject)da,(x+2*(sDist))*sizeof(PetscInt));
223: for (i=0; i<IXs-Xs; i++) idx[i] = -1; /* prepend with -1s if needed for ghosted case*/
225: nn = IXs-Xs;
226: if (bx == DM_BOUNDARY_PERIODIC) { /* Handle all cases with periodic first */
227: for (i=0; i<sDist; i++) { /* Left ghost points */
228: if ((xs-sDist+i)>=0) idx[nn++] = xs-sDist+i;
229: else idx[nn++] = M+(xs-sDist+i);
230: }
232: for (i=0; i<x; i++) idx [nn++] = xs + i; /* Non-ghost points */
234: for (i=0; i<sDist; i++) { /* Right ghost points */
235: if ((xe+i)<M) idx [nn++] = xe+i;
236: else idx [nn++] = (xe+i) - M;
237: }
238: } else if (bx == DM_BOUNDARY_MIRROR) { /* Handle all cases with periodic first */
239: for (i=0; i<(sDist); i++) { /* Left ghost points */
240: if ((xs-sDist+i)>=0) idx[nn++] = xs-sDist+i;
241: else idx[nn++] = sDist - i;
242: }
244: for (i=0; i<x; i++) idx [nn++] = xs + i; /* Non-ghost points */
246: for (i=0; i<(sDist); i++) { /* Right ghost points */
247: if ((xe+i)<M) idx[nn++] = xe+i;
248: else idx[nn++] = M - (i + 1);
249: }
250: } else { /* Now do all cases with no periodicity */
251: if (0 <= xs-sDist) {
252: for (i=0; i<sDist; i++) idx[nn++] = xs - sDist + i;
253: } else {
254: for (i=0; i<xs; i++) idx[nn++] = i;
255: }
257: for (i=0; i<x; i++) idx [nn++] = xs + i;
259: if ((xe+sDist)<=M) {
260: for (i=0; i<sDist; i++) idx[nn++]=xe+i;
261: } else {
262: for (i=xe; i<M; i++) idx[nn++]=i;
263: }
264: }
266: ISCreateBlock(comm,dof,nn-IXs+Xs,&idx[IXs-Xs],PETSC_USE_POINTER,&from);
267: VecScatterCreate(global,from,local,to,>ol);
268: PetscLogObjectParent((PetscObject)da,(PetscObject)gtol);
269: ISDestroy(&to);
270: ISDestroy(&from);
271: VecDestroy(&local);
272: VecDestroy(&global);
274: dd->xs = dof*xs; dd->xe = dof*xe; dd->ys = 0; dd->ye = 1; dd->zs = 0; dd->ze = 1;
275: dd->Xs = dof*Xs; dd->Xe = dof*Xe; dd->Ys = 0; dd->Ye = 1; dd->Zs = 0; dd->Ze = 1;
277: dd->gtol = gtol;
278: dd->base = dof*xs;
279: da->ops->view = DMView_DA_1d;
281: /*
282: Set the local to global ordering in the global vector, this allows use
283: of VecSetValuesLocal().
284: */
285: for (i=0; i<Xe-IXe; i++) idx[nn++] = -1; /* pad with -1s if needed for ghosted case*/
287: ISLocalToGlobalMappingCreate(comm,dof,nn,idx,PETSC_OWN_POINTER,&da->ltogmap);
288: PetscLogObjectParent((PetscObject)da,(PetscObject)da->ltogmap);
290: return(0);
291: }
296: /*@C
297: DMDACreate1d - Creates an object that will manage the communication of one-dimensional
298: regular array data that is distributed across some processors.
300: Collective on MPI_Comm
302: Input Parameters:
303: + comm - MPI communicator
304: . bx - type of ghost cells at the boundary the array should have, if any. Use
305: DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, or DM_BOUNDARY_PERIODIC.
306: . M - global dimension of the array (use -M to indicate that it may be set to a different value
307: from the command line with -da_grid_x <M>)
308: . dof - number of degrees of freedom per node
309: . s - stencil width
310: - lx - array containing number of nodes in the X direction on each processor,
311: or NULL. If non-null, must be of length as the number of processes in the MPI_Comm.
313: Output Parameter:
314: . da - the resulting distributed array object
316: Options Database Key:
317: + -dm_view - Calls DMView() at the conclusion of DMDACreate1d()
318: . -da_grid_x <nx> - number of grid points in x direction; can set if M < 0
319: . -da_refine_x <rx> - refinement factor
320: - -da_refine <n> - refine the DMDA n times before creating it, if M < 0
322: Level: beginner
324: Notes:
325: The array data itself is NOT stored in the DMDA, it is stored in Vec objects;
326: The appropriate vector objects can be obtained with calls to DMCreateGlobalVector()
327: and DMCreateLocalVector() and calls to VecDuplicate() if more are needed.
329: .keywords: distributed array, create, one-dimensional
331: .seealso: DMDestroy(), DMView(), DMDACreate2d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDASetRefinementFactor(),
332: DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToLocalBegin(), DMLocalToLocalEnd(), DMDAGetRefinementFactor(),
333: DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges()
335: @*/
336: PetscErrorCode DMDACreate1d(MPI_Comm comm, DMBoundaryType bx, PetscInt M, PetscInt dof, PetscInt s, const PetscInt lx[], DM *da)
337: {
339: PetscMPIInt size;
342: DMDACreate(comm, da);
343: DMSetDimension(*da, 1);
344: DMDASetSizes(*da, M, 1, 1);
345: MPI_Comm_size(comm, &size);
346: DMDASetNumProcs(*da, size, PETSC_DECIDE, PETSC_DECIDE);
347: DMDASetBoundaryType(*da, bx, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE);
348: DMDASetDof(*da, dof);
349: DMDASetStencilWidth(*da, s);
350: DMDASetOwnershipRanges(*da, lx, NULL, NULL);
351: /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */
352: DMSetFromOptions(*da);
353: DMSetUp(*da);
354: return(0);
355: }