Actual source code: da.c
petsc-3.6.4 2016-04-12
1: #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
5: /*@
6: DMDASetSizes - Sets the global sizes
8: Logically Collective on DMDA
10: Input Parameters:
11: + da - the DMDA
12: . M - the global X size (or PETSC_DECIDE)
13: . N - the global Y size (or PETSC_DECIDE)
14: - P - the global Z size (or PETSC_DECIDE)
16: Level: intermediate
18: .seealso: DMDAGetSize(), PetscSplitOwnership()
19: @*/
20: PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
21: {
22: DM_DA *dd = (DM_DA*)da->data;
29: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
31: dd->M = M;
32: dd->N = N;
33: dd->P = P;
34: return(0);
35: }
39: /*@
40: DMDASetNumProcs - Sets the number of processes in each dimension
42: Logically Collective on DMDA
44: Input Parameters:
45: + da - the DMDA
46: . m - the number of X procs (or PETSC_DECIDE)
47: . n - the number of Y procs (or PETSC_DECIDE)
48: - p - the number of Z procs (or PETSC_DECIDE)
50: Level: intermediate
52: .seealso: DMDASetSizes(), DMDAGetSize(), PetscSplitOwnership()
53: @*/
54: PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
55: {
56: DM_DA *dd = (DM_DA*)da->data;
64: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
65: dd->m = m;
66: dd->n = n;
67: dd->p = p;
68: if (da->dim == 2) {
69: PetscMPIInt size;
70: MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);
71: if ((dd->m > 0) && (dd->n < 0)) {
72: dd->n = size/dd->m;
73: if (dd->n*dd->m != size) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"%D processes in X direction not divisible into comm size %d",m,size);
74: }
75: if ((dd->n > 0) && (dd->m < 0)) {
76: dd->m = size/dd->n;
77: if (dd->n*dd->m != size) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"%D processes in Y direction not divisible into comm size %d",n,size);
78: }
79: }
80: return(0);
81: }
85: /*@
86: DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.
88: Not collective
90: Input Parameter:
91: + da - The DMDA
92: - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC
94: Level: intermediate
96: .keywords: distributed array, periodicity
97: .seealso: DMDACreate(), DMDestroy(), DMDA, DMBoundaryType
98: @*/
99: PetscErrorCode DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)
100: {
101: DM_DA *dd = (DM_DA*)da->data;
108: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
109: dd->bx = bx;
110: dd->by = by;
111: dd->bz = bz;
112: return(0);
113: }
117: /*@
118: DMDASetDof - Sets the number of degrees of freedom per vertex
120: Not collective
122: Input Parameter:
123: + da - The DMDA
124: - dof - Number of degrees of freedom
126: Level: intermediate
128: .keywords: distributed array, degrees of freedom
129: .seealso: DMDACreate(), DMDestroy(), DMDA
130: @*/
131: PetscErrorCode DMDASetDof(DM da, PetscInt dof)
132: {
133: DM_DA *dd = (DM_DA*)da->data;
138: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
139: dd->w = dof;
140: da->bs = dof;
141: return(0);
142: }
146: /*@
147: DMDAGetOverlap - Gets the size of the per-processor overlap.
149: Not collective
151: Input Parameters:
152: . da - The DMDA
154: Output Parameters:
155: + x - Overlap in the x direction
156: . y - Overlap in the y direction
157: - z - Overlap in the z direction
159: Level: intermediate
161: .keywords: distributed array, overlap, domain decomposition
162: .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
163: @*/
164: PetscErrorCode DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
165: {
166: DM_DA *dd = (DM_DA*)da->data;
170: if (x) *x = dd->xol;
171: if (y) *y = dd->yol;
172: if (z) *z = dd->zol;
173: return(0);
174: }
178: /*@
179: DMDASetOverlap - Sets the size of the per-processor overlap.
181: Not collective
183: Input Parameters:
184: + da - The DMDA
185: . x - Overlap in the x direction
186: . y - Overlap in the y direction
187: - z - Overlap in the z direction
189: Level: intermediate
191: .keywords: distributed array, overlap, domain decomposition
192: .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
193: @*/
194: PetscErrorCode DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
195: {
196: DM_DA *dd = (DM_DA*)da->data;
203: dd->xol = x;
204: dd->yol = y;
205: dd->zol = z;
206: return(0);
207: }
212: /*@
213: DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.
215: Not collective
217: Input Parameters:
218: . da - The DMDA
220: Output Parameters:
221: + Nsub - Number of local subdomains created upon decomposition
223: Level: intermediate
225: .keywords: distributed array, domain decomposition
226: .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
227: @*/
228: PetscErrorCode DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
229: {
230: DM_DA *dd = (DM_DA*)da->data;
234: if (Nsub) *Nsub = dd->Nsub;
235: return(0);
236: }
240: /*@
241: DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.
243: Not collective
245: Input Parameters:
246: + da - The DMDA
247: - Nsub - The number of local subdomains requested
249: Level: intermediate
251: .keywords: distributed array, domain decomposition
252: .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
253: @*/
254: PetscErrorCode DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
255: {
256: DM_DA *dd = (DM_DA*)da->data;
261: dd->Nsub = Nsub;
262: return(0);
263: }
267: /*@
268: DMDASetOffset - Sets the index offset of the DA.
270: Collective on DA
272: Input Parameter:
273: + da - The DMDA
274: . xo - The offset in the x direction
275: . yo - The offset in the y direction
276: - zo - The offset in the z direction
278: Level: intermediate
280: Notes: This is used primarily to overlap a computation on a local DA with that on a global DA without
281: changing boundary conditions or subdomain features that depend upon the global offsets.
283: .keywords: distributed array, degrees of freedom
284: .seealso: DMDAGetOffset(), DMDAVecGetArray()
285: @*/
286: PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
287: {
289: DM_DA *dd = (DM_DA*)da->data;
299: dd->xo = xo;
300: dd->yo = yo;
301: dd->zo = zo;
302: dd->Mo = Mo;
303: dd->No = No;
304: dd->Po = Po;
306: if (da->coordinateDM) {
307: DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);
308: }
309: return(0);
310: }
314: /*@
315: DMDAGetOffset - Gets the index offset of the DA.
317: Not collective
319: Input Parameter:
320: . da - The DMDA
322: Output Parameters:
323: + xo - The offset in the x direction
324: . yo - The offset in the y direction
325: . zo - The offset in the z direction
326: . Mo - The global size in the x direction
327: . No - The global size in the y direction
328: - Po - The global size in the z direction
330: Level: intermediate
332: .keywords: distributed array, degrees of freedom
333: .seealso: DMDASetOffset(), DMDAVecGetArray()
334: @*/
335: PetscErrorCode DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
336: {
337: DM_DA *dd = (DM_DA*)da->data;
341: if (xo) *xo = dd->xo;
342: if (yo) *yo = dd->yo;
343: if (zo) *zo = dd->zo;
344: if (Mo) *Mo = dd->Mo;
345: if (No) *No = dd->No;
346: if (Po) *Po = dd->Po;
347: return(0);
348: }
352: /*@
353: DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.
355: Not collective
357: Input Parameter:
358: . da - The DMDA
360: Output Parameters:
361: + xs - The start of the region in x
362: . ys - The start of the region in y
363: . zs - The start of the region in z
364: . xs - The size of the region in x
365: . ys - The size of the region in y
366: . zs - The size of the region in z
368: Level: intermediate
370: .keywords: distributed array, degrees of freedom
371: .seealso: DMDAGetOffset(), DMDAVecGetArray()
372: @*/
373: PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
374: {
375: DM_DA *dd = (DM_DA*)da->data;
379: if (xs) *xs = dd->nonxs;
380: if (ys) *ys = dd->nonys;
381: if (zs) *zs = dd->nonzs;
382: if (xm) *xm = dd->nonxm;
383: if (ym) *ym = dd->nonym;
384: if (zm) *zm = dd->nonzm;
385: return(0);
386: }
391: /*@
392: DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.
394: Collective on DA
396: Input Parameter:
397: + da - The DMDA
398: . xs - The start of the region in x
399: . ys - The start of the region in y
400: . zs - The start of the region in z
401: . xs - The size of the region in x
402: . ys - The size of the region in y
403: . zs - The size of the region in z
405: Level: intermediate
407: .keywords: distributed array, degrees of freedom
408: .seealso: DMDAGetOffset(), DMDAVecGetArray()
409: @*/
410: PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
411: {
412: DM_DA *dd = (DM_DA*)da->data;
422: dd->nonxs = xs;
423: dd->nonys = ys;
424: dd->nonzs = zs;
425: dd->nonxm = xm;
426: dd->nonym = ym;
427: dd->nonzm = zm;
429: return(0);
430: }
434: /*@
435: DMDASetStencilType - Sets the type of the communication stencil
437: Logically Collective on DMDA
439: Input Parameter:
440: + da - The DMDA
441: - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
443: Level: intermediate
445: .keywords: distributed array, stencil
446: .seealso: DMDACreate(), DMDestroy(), DMDA
447: @*/
448: PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
449: {
450: DM_DA *dd = (DM_DA*)da->data;
455: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
456: dd->stencil_type = stype;
457: return(0);
458: }
462: /*@
463: DMDASetStencilWidth - Sets the width of the communication stencil
465: Logically Collective on DMDA
467: Input Parameter:
468: + da - The DMDA
469: - width - The stencil width
471: Level: intermediate
473: .keywords: distributed array, stencil
474: .seealso: DMDACreate(), DMDestroy(), DMDA
475: @*/
476: PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
477: {
478: DM_DA *dd = (DM_DA*)da->data;
483: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
484: dd->s = width;
485: return(0);
486: }
490: static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
491: {
492: PetscInt i,sum;
495: if (M < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
496: for (i=sum=0; i<m; i++) sum += lx[i];
497: if (sum != M) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
498: return(0);
499: }
503: /*@
504: DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
506: Logically Collective on DMDA
508: Input Parameter:
509: + da - The DMDA
510: . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
511: . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
512: - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.
514: Level: intermediate
516: Note: these numbers are NOT multiplied by the number of dof per node.
518: .keywords: distributed array
519: .seealso: DMDACreate(), DMDestroy(), DMDA
520: @*/
521: PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
522: {
524: DM_DA *dd = (DM_DA*)da->data;
528: if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
529: if (lx) {
530: if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
531: DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);
532: if (!dd->lx) {
533: PetscMalloc1(dd->m, &dd->lx);
534: }
535: PetscMemcpy(dd->lx, lx, dd->m*sizeof(PetscInt));
536: }
537: if (ly) {
538: if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
539: DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);
540: if (!dd->ly) {
541: PetscMalloc1(dd->n, &dd->ly);
542: }
543: PetscMemcpy(dd->ly, ly, dd->n*sizeof(PetscInt));
544: }
545: if (lz) {
546: if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
547: DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);
548: if (!dd->lz) {
549: PetscMalloc1(dd->p, &dd->lz);
550: }
551: PetscMemcpy(dd->lz, lz, dd->p*sizeof(PetscInt));
552: }
553: return(0);
554: }
558: /*@
559: DMDASetInterpolationType - Sets the type of interpolation that will be
560: returned by DMCreateInterpolation()
562: Logically Collective on DMDA
564: Input Parameter:
565: + da - initial distributed array
566: . ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
568: Level: intermediate
570: Notes: you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
572: .keywords: distributed array, interpolation
574: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
575: @*/
576: PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
577: {
578: DM_DA *dd = (DM_DA*)da->data;
583: dd->interptype = ctype;
584: return(0);
585: }
589: /*@
590: DMDAGetInterpolationType - Gets the type of interpolation that will be
591: used by DMCreateInterpolation()
593: Not Collective
595: Input Parameter:
596: . da - distributed array
598: Output Parameter:
599: . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
601: Level: intermediate
603: .keywords: distributed array, interpolation
605: .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
606: @*/
607: PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
608: {
609: DM_DA *dd = (DM_DA*)da->data;
614: *ctype = dd->interptype;
615: return(0);
616: }
620: /*@C
621: DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
622: processes neighbors.
624: Not Collective
626: Input Parameter:
627: . da - the DMDA object
629: Output Parameters:
630: . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
631: this process itself is in the list
633: Notes: In 2d the array is of length 9, in 3d of length 27
634: Not supported in 1d
635: Do not free the array, it is freed when the DMDA is destroyed.
637: Fortran Notes: In fortran you must pass in an array of the appropriate length.
639: Level: intermediate
641: @*/
642: PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
643: {
644: DM_DA *dd = (DM_DA*)da->data;
648: *ranks = dd->neighbors;
649: return(0);
650: }
654: /*@C
655: DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
657: Not Collective
659: Input Parameter:
660: . da - the DMDA object
662: Output Parameter:
663: + lx - ownership along x direction (optional)
664: . ly - ownership along y direction (optional)
665: - lz - ownership along z direction (optional)
667: Level: intermediate
669: Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
671: In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
672: eighth arguments from DMDAGetInfo()
674: In C you should not free these arrays, nor change the values in them. They will only have valid values while the
675: DMDA they came from still exists (has not been destroyed).
677: These numbers are NOT multiplied by the number of dof per node.
679: .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
680: @*/
681: PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
682: {
683: DM_DA *dd = (DM_DA*)da->data;
687: if (lx) *lx = dd->lx;
688: if (ly) *ly = dd->ly;
689: if (lz) *lz = dd->lz;
690: return(0);
691: }
695: /*@
696: DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
698: Logically Collective on DMDA
700: Input Parameters:
701: + da - the DMDA object
702: . refine_x - ratio of fine grid to coarse in x direction (2 by default)
703: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
704: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
706: Options Database:
707: + -da_refine_x - refinement ratio in x direction
708: . -da_refine_y - refinement ratio in y direction
709: - -da_refine_z - refinement ratio in z direction
711: Level: intermediate
713: Notes: Pass PETSC_IGNORE to leave a value unchanged
715: .seealso: DMRefine(), DMDAGetRefinementFactor()
716: @*/
717: PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
718: {
719: DM_DA *dd = (DM_DA*)da->data;
727: if (refine_x > 0) dd->refine_x = refine_x;
728: if (refine_y > 0) dd->refine_y = refine_y;
729: if (refine_z > 0) dd->refine_z = refine_z;
730: return(0);
731: }
735: /*@C
736: DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
738: Not Collective
740: Input Parameter:
741: . da - the DMDA object
743: Output Parameters:
744: + refine_x - ratio of fine grid to coarse in x direction (2 by default)
745: . refine_y - ratio of fine grid to coarse in y direction (2 by default)
746: - refine_z - ratio of fine grid to coarse in z direction (2 by default)
748: Level: intermediate
750: Notes: Pass NULL for values you do not need
752: .seealso: DMRefine(), DMDASetRefinementFactor()
753: @*/
754: PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
755: {
756: DM_DA *dd = (DM_DA*)da->data;
760: if (refine_x) *refine_x = dd->refine_x;
761: if (refine_y) *refine_y = dd->refine_y;
762: if (refine_z) *refine_z = dd->refine_z;
763: return(0);
764: }
768: /*@C
769: DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
771: Logically Collective on DMDA
773: Input Parameters:
774: + da - the DMDA object
775: - f - the function that allocates the matrix for that specific DMDA
777: Level: developer
779: Notes: See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
780: the diagonal and off-diagonal blocks of the matrix
782: .seealso: DMCreateMatrix(), DMDASetBlockFills()
783: @*/
784: PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
785: {
788: da->ops->creatematrix = f;
789: return(0);
790: }
794: /*
795: Creates "balanced" ownership ranges after refinement, constrained by the need for the
796: fine grid boundaries to fall within one stencil width of the coarse partition.
798: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
799: */
800: static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
801: {
802: PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
806: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
807: if (ratio == 1) {
808: PetscMemcpy(lf,lc,m*sizeof(lc[0]));
809: return(0);
810: }
811: for (i=0; i<m; i++) totalc += lc[i];
812: remaining = (!periodic) + ratio * (totalc - (!periodic));
813: for (i=0; i<m; i++) {
814: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
815: if (i == m-1) lf[i] = want;
816: else {
817: const PetscInt nextc = startc + lc[i];
818: /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
819: * coarse stencil width of the first coarse node in the next subdomain. */
820: while ((startf+want)/ratio < nextc - stencil_width) want++;
821: /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
822: * coarse stencil width of the last coarse node in the current subdomain. */
823: while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
824: /* Make sure all constraints are satisfied */
825: if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
826: || ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible refined ownership range");
827: }
828: lf[i] = want;
829: startc += lc[i];
830: startf += lf[i];
831: remaining -= lf[i];
832: }
833: return(0);
834: }
838: /*
839: Creates "balanced" ownership ranges after coarsening, constrained by the need for the
840: fine grid boundaries to fall within one stencil width of the coarse partition.
842: Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
843: */
844: static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
845: {
846: PetscInt i,totalf,remaining,startc,startf;
850: if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
851: if (ratio == 1) {
852: PetscMemcpy(lc,lf,m*sizeof(lf[0]));
853: return(0);
854: }
855: for (i=0,totalf=0; i<m; i++) totalf += lf[i];
856: remaining = (!periodic) + (totalf - (!periodic)) / ratio;
857: for (i=0,startc=0,startf=0; i<m; i++) {
858: PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
859: if (i == m-1) lc[i] = want;
860: else {
861: const PetscInt nextf = startf+lf[i];
862: /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
863: * node is within one stencil width. */
864: while (nextf/ratio < startc+want-stencil_width) want--;
865: /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
866: * fine node is within one stencil width. */
867: while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
868: if (want < 0 || want > remaining
869: || (nextf/ratio < startc+want-stencil_width) || ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible coarsened ownership range");
870: }
871: lc[i] = want;
872: startc += lc[i];
873: startf += lf[i];
874: remaining -= lc[i];
875: }
876: return(0);
877: }
881: PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
882: {
884: PetscInt M,N,P,i,dim;
885: DM da2;
886: DM_DA *dd = (DM_DA*)da->data,*dd2;
892: DMGetDimension(da, &dim);
893: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
894: M = dd->refine_x*dd->M;
895: } else {
896: M = 1 + dd->refine_x*(dd->M - 1);
897: }
898: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
899: if (dim > 1) {
900: N = dd->refine_y*dd->N;
901: } else {
902: N = 1;
903: }
904: } else {
905: N = 1 + dd->refine_y*(dd->N - 1);
906: }
907: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
908: if (dim > 2) {
909: P = dd->refine_z*dd->P;
910: } else {
911: P = 1;
912: }
913: } else {
914: P = 1 + dd->refine_z*(dd->P - 1);
915: }
916: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
917: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
918: DMSetDimension(da2,dim);
919: DMDASetSizes(da2,M,N,P);
920: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
921: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
922: DMDASetDof(da2,dd->w);
923: DMDASetStencilType(da2,dd->stencil_type);
924: DMDASetStencilWidth(da2,dd->s);
925: if (dim == 3) {
926: PetscInt *lx,*ly,*lz;
927: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
928: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
929: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
930: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);
931: DMDASetOwnershipRanges(da2,lx,ly,lz);
932: PetscFree3(lx,ly,lz);
933: } else if (dim == 2) {
934: PetscInt *lx,*ly;
935: PetscMalloc2(dd->m,&lx,dd->n,&ly);
936: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
937: DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);
938: DMDASetOwnershipRanges(da2,lx,ly,NULL);
939: PetscFree2(lx,ly);
940: } else if (dim == 1) {
941: PetscInt *lx;
942: PetscMalloc1(dd->m,&lx);
943: DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);
944: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
945: PetscFree(lx);
946: }
947: dd2 = (DM_DA*)da2->data;
949: /* allow overloaded (user replaced) operations to be inherited by refinement clones */
950: da2->ops->creatematrix = da->ops->creatematrix;
951: /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
952: da2->ops->getcoloring = da->ops->getcoloring;
953: dd2->interptype = dd->interptype;
955: /* copy fill information if given */
956: if (dd->dfill) {
957: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
958: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
959: }
960: if (dd->ofill) {
961: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
962: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
963: }
964: /* copy the refine information */
965: dd2->coarsen_x = dd2->refine_x = dd->refine_x;
966: dd2->coarsen_y = dd2->refine_y = dd->refine_y;
967: dd2->coarsen_z = dd2->refine_z = dd->refine_z;
969: /* copy vector type information */
970: DMSetVecType(da2,da->vectype);
972: dd2->lf = dd->lf;
973: dd2->lj = dd->lj;
975: da2->leveldown = da->leveldown;
976: da2->levelup = da->levelup + 1;
978: DMSetFromOptions(da2);
979: DMSetUp(da2);
981: /* interpolate coordinates if they are set on the coarse grid */
982: if (da->coordinates) {
983: DM cdaf,cdac;
984: Vec coordsc,coordsf;
985: Mat II;
987: DMGetCoordinateDM(da,&cdac);
988: DMGetCoordinates(da,&coordsc);
989: DMGetCoordinateDM(da2,&cdaf);
990: /* force creation of the coordinate vector */
991: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
992: DMGetCoordinates(da2,&coordsf);
993: DMCreateInterpolation(cdac,cdaf,&II,NULL);
994: MatInterpolate(II,coordsc,coordsf);
995: MatDestroy(&II);
996: }
998: for (i=0; i<da->bs; i++) {
999: const char *fieldname;
1000: DMDAGetFieldName(da,i,&fieldname);
1001: DMDASetFieldName(da2,i,fieldname);
1002: }
1004: *daref = da2;
1005: return(0);
1006: }
1011: PetscErrorCode DMCoarsen_DA(DM da, MPI_Comm comm,DM *daref)
1012: {
1014: PetscInt M,N,P,i,dim;
1015: DM da2;
1016: DM_DA *dd = (DM_DA*)da->data,*dd2;
1022: DMGetDimension(da, &dim);
1023: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1024: M = dd->M / dd->coarsen_x;
1025: } else {
1026: M = 1 + (dd->M - 1) / dd->coarsen_x;
1027: }
1028: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1029: if (dim > 1) {
1030: N = dd->N / dd->coarsen_y;
1031: } else {
1032: N = 1;
1033: }
1034: } else {
1035: N = 1 + (dd->N - 1) / dd->coarsen_y;
1036: }
1037: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1038: if (dim > 2) {
1039: P = dd->P / dd->coarsen_z;
1040: } else {
1041: P = 1;
1042: }
1043: } else {
1044: P = 1 + (dd->P - 1) / dd->coarsen_z;
1045: }
1046: DMDACreate(PetscObjectComm((PetscObject)da),&da2);
1047: DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);
1048: DMSetDimension(da2,dim);
1049: DMDASetSizes(da2,M,N,P);
1050: DMDASetNumProcs(da2,dd->m,dd->n,dd->p);
1051: DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);
1052: DMDASetDof(da2,dd->w);
1053: DMDASetStencilType(da2,dd->stencil_type);
1054: DMDASetStencilWidth(da2,dd->s);
1055: if (dim == 3) {
1056: PetscInt *lx,*ly,*lz;
1057: PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);
1058: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1059: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1060: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);
1061: DMDASetOwnershipRanges(da2,lx,ly,lz);
1062: PetscFree3(lx,ly,lz);
1063: } else if (dim == 2) {
1064: PetscInt *lx,*ly;
1065: PetscMalloc2(dd->m,&lx,dd->n,&ly);
1066: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1067: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);
1068: DMDASetOwnershipRanges(da2,lx,ly,NULL);
1069: PetscFree2(lx,ly);
1070: } else if (dim == 1) {
1071: PetscInt *lx;
1072: PetscMalloc1(dd->m,&lx);
1073: DMDACoarsenOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);
1074: DMDASetOwnershipRanges(da2,lx,NULL,NULL);
1075: PetscFree(lx);
1076: }
1077: dd2 = (DM_DA*)da2->data;
1079: /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1080: /* da2->ops->createinterpolation = da->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1081: da2->ops->creatematrix = da->ops->creatematrix;
1082: da2->ops->getcoloring = da->ops->getcoloring;
1083: dd2->interptype = dd->interptype;
1085: /* copy fill information if given */
1086: if (dd->dfill) {
1087: PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);
1088: PetscMemcpy(dd2->dfill,dd->dfill,(dd->dfill[dd->w]+dd->w+1)*sizeof(PetscInt));
1089: }
1090: if (dd->ofill) {
1091: PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);
1092: PetscMemcpy(dd2->ofill,dd->ofill,(dd->ofill[dd->w]+dd->w+1)*sizeof(PetscInt));
1093: }
1094: /* copy the refine information */
1095: dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1096: dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1097: dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1099: /* copy vector type information */
1100: DMSetVecType(da2,da->vectype);
1102: dd2->lf = dd->lf;
1103: dd2->lj = dd->lj;
1105: da2->leveldown = da->leveldown + 1;
1106: da2->levelup = da->levelup;
1108: DMSetFromOptions(da2);
1109: DMSetUp(da2);
1111: /* inject coordinates if they are set on the fine grid */
1112: if (da->coordinates) {
1113: DM cdaf,cdac;
1114: Vec coordsc,coordsf;
1115: Mat inject;
1116: VecScatter vscat;
1118: DMGetCoordinateDM(da,&cdaf);
1119: DMGetCoordinates(da,&coordsf);
1120: DMGetCoordinateDM(da2,&cdac);
1121: /* force creation of the coordinate vector */
1122: DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);
1123: DMGetCoordinates(da2,&coordsc);
1125: DMCreateInjection(cdac,cdaf,&inject);
1126: MatScatterGetVecScatter(inject,&vscat);
1127: VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1128: VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);
1129: MatDestroy(&inject);
1130: }
1132: for (i=0; i<da->bs; i++) {
1133: const char *fieldname;
1134: DMDAGetFieldName(da,i,&fieldname);
1135: DMDASetFieldName(da2,i,fieldname);
1136: }
1138: *daref = da2;
1139: return(0);
1140: }
1144: PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1145: {
1147: PetscInt i,n,*refx,*refy,*refz;
1151: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1152: if (nlevels == 0) return(0);
1155: /* Get refinement factors, defaults taken from the coarse DMDA */
1156: PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);
1157: for (i=0; i<nlevels; i++) {
1158: DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);
1159: }
1160: n = nlevels;
1161: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);
1162: n = nlevels;
1163: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);
1164: n = nlevels;
1165: PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);
1167: DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);
1168: DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);
1169: for (i=1; i<nlevels; i++) {
1170: DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);
1171: DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);
1172: }
1173: PetscFree3(refx,refy,refz);
1174: return(0);
1175: }
1179: PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1180: {
1182: PetscInt i;
1186: if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1187: if (nlevels == 0) return(0);
1189: DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);
1190: for (i=1; i<nlevels; i++) {
1191: DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);
1192: }
1193: return(0);
1194: }