Actual source code: dacreate.c
1: #include <petsc/private/dmdaimpl.h>
3: static PetscErrorCode DMSetFromOptions_DA(DM da, PetscOptionItems *PetscOptionsObject)
4: {
5: DM_DA *dd = (DM_DA *)da->data;
6: PetscInt refine = 0, dim = da->dim, maxnlevels = 100, refx[100], refy[100], refz[100], n, i;
7: PetscBool flg;
9: PetscFunctionBegin;
10: PetscCheck(dd->M >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Dimension must be non-negative, call DMSetFromOptions() if you want to change the value at runtime");
11: PetscCheck(dd->N >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Dimension must be non-negative, call DMSetFromOptions() if you want to change the value at runtime");
12: PetscCheck(dd->P >= 0, PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_OUTOFRANGE, "Dimension must be non-negative, call DMSetFromOptions() if you want to change the value at runtime");
14: PetscOptionsHeadBegin(PetscOptionsObject, "DMDA Options");
15: PetscCall(PetscOptionsBoundedInt("-da_grid_x", "Number of grid points in x direction", "DMDASetSizes", dd->M, &dd->M, NULL, 1));
16: if (dim > 1) PetscCall(PetscOptionsBoundedInt("-da_grid_y", "Number of grid points in y direction", "DMDASetSizes", dd->N, &dd->N, NULL, 1));
17: if (dim > 2) PetscCall(PetscOptionsBoundedInt("-da_grid_z", "Number of grid points in z direction", "DMDASetSizes", dd->P, &dd->P, NULL, 1));
19: PetscCall(PetscOptionsBoundedInt("-da_overlap", "Decomposition overlap in all directions", "DMDASetOverlap", dd->xol, &dd->xol, &flg, 0));
20: if (flg) PetscCall(DMDASetOverlap(da, dd->xol, dd->xol, dd->xol));
21: PetscCall(PetscOptionsBoundedInt("-da_overlap_x", "Decomposition overlap in x direction", "DMDASetOverlap", dd->xol, &dd->xol, NULL, 0));
22: if (dim > 1) PetscCall(PetscOptionsBoundedInt("-da_overlap_y", "Decomposition overlap in y direction", "DMDASetOverlap", dd->yol, &dd->yol, NULL, 0));
23: if (dim > 2) PetscCall(PetscOptionsBoundedInt("-da_overlap_z", "Decomposition overlap in z direction", "DMDASetOverlap", dd->zol, &dd->zol, NULL, 0));
25: PetscCall(PetscOptionsBoundedInt("-da_local_subdomains", "", "DMDASetNumLocalSubdomains", dd->Nsub, &dd->Nsub, &flg, PETSC_DECIDE));
26: if (flg) PetscCall(DMDASetNumLocalSubDomains(da, dd->Nsub));
28: /* Handle DMDA parallel distribution */
29: PetscCall(PetscOptionsBoundedInt("-da_processors_x", "Number of processors in x direction", "DMDASetNumProcs", dd->m, &dd->m, NULL, PETSC_DECIDE));
30: if (dim > 1) PetscCall(PetscOptionsBoundedInt("-da_processors_y", "Number of processors in y direction", "DMDASetNumProcs", dd->n, &dd->n, NULL, PETSC_DECIDE));
31: if (dim > 2) PetscCall(PetscOptionsBoundedInt("-da_processors_z", "Number of processors in z direction", "DMDASetNumProcs", dd->p, &dd->p, NULL, PETSC_DECIDE));
32: /* Handle DMDA refinement */
33: PetscCall(PetscOptionsBoundedInt("-da_refine_x", "Refinement ratio in x direction", "DMDASetRefinementFactor", dd->refine_x, &dd->refine_x, NULL, 1));
34: if (dim > 1) PetscCall(PetscOptionsBoundedInt("-da_refine_y", "Refinement ratio in y direction", "DMDASetRefinementFactor", dd->refine_y, &dd->refine_y, NULL, 1));
35: if (dim > 2) PetscCall(PetscOptionsBoundedInt("-da_refine_z", "Refinement ratio in z direction", "DMDASetRefinementFactor", dd->refine_z, &dd->refine_z, NULL, 1));
36: dd->coarsen_x = dd->refine_x;
37: dd->coarsen_y = dd->refine_y;
38: dd->coarsen_z = dd->refine_z;
40: /* Get refinement factors, defaults taken from the coarse DMDA */
41: PetscCall(DMDAGetRefinementFactor(da, &refx[0], &refy[0], &refz[0]));
42: for (i = 1; i < maxnlevels; i++) {
43: refx[i] = refx[0];
44: refy[i] = refy[0];
45: refz[i] = refz[0];
46: }
47: n = maxnlevels;
48: PetscCall(PetscOptionsIntArray("-da_refine_hierarchy_x", "Refinement factor for each level", "None", refx, &n, &flg));
49: if (flg) {
50: dd->refine_x = refx[0];
51: dd->refine_x_hier_n = n;
52: PetscCall(PetscMalloc1(n, &dd->refine_x_hier));
53: PetscCall(PetscArraycpy(dd->refine_x_hier, refx, n));
54: }
55: if (dim > 1) {
56: n = maxnlevels;
57: PetscCall(PetscOptionsIntArray("-da_refine_hierarchy_y", "Refinement factor for each level", "None", refy, &n, &flg));
58: if (flg) {
59: dd->refine_y = refy[0];
60: dd->refine_y_hier_n = n;
61: PetscCall(PetscMalloc1(n, &dd->refine_y_hier));
62: PetscCall(PetscArraycpy(dd->refine_y_hier, refy, n));
63: }
64: }
65: if (dim > 2) {
66: n = maxnlevels;
67: PetscCall(PetscOptionsIntArray("-da_refine_hierarchy_z", "Refinement factor for each level", "None", refz, &n, &flg));
68: if (flg) {
69: dd->refine_z = refz[0];
70: dd->refine_z_hier_n = n;
71: PetscCall(PetscMalloc1(n, &dd->refine_z_hier));
72: PetscCall(PetscArraycpy(dd->refine_z_hier, refz, n));
73: }
74: }
76: PetscCall(PetscOptionsBoundedInt("-da_refine", "Uniformly refine DA one or more times", "None", refine, &refine, NULL, 0));
77: PetscOptionsHeadEnd();
79: while (refine--) {
80: if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
81: PetscCall(PetscIntMultError(dd->refine_x, dd->M, &dd->M));
82: } else {
83: PetscCall(PetscIntMultError(dd->refine_x, dd->M - 1, &dd->M));
84: dd->M += 1;
85: }
86: if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
87: PetscCall(PetscIntMultError(dd->refine_y, dd->N, &dd->N));
88: } else {
89: PetscCall(PetscIntMultError(dd->refine_y, dd->N - 1, &dd->N));
90: dd->N += 1;
91: }
92: if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
93: PetscCall(PetscIntMultError(dd->refine_z, dd->P, &dd->P));
94: } else {
95: PetscCall(PetscIntMultError(dd->refine_z, dd->P - 1, &dd->P));
96: dd->P += 1;
97: }
98: da->levelup++;
99: if (da->levelup - da->leveldown >= 0) {
100: dd->refine_x = refx[da->levelup - da->leveldown];
101: dd->refine_y = refy[da->levelup - da->leveldown];
102: dd->refine_z = refz[da->levelup - da->leveldown];
103: }
104: if (da->levelup - da->leveldown >= 1) {
105: dd->coarsen_x = refx[da->levelup - da->leveldown - 1];
106: dd->coarsen_y = refy[da->levelup - da->leveldown - 1];
107: dd->coarsen_z = refz[da->levelup - da->leveldown - 1];
108: }
109: }
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: static PetscErrorCode DMLoad_DA(DM da, PetscViewer viewer)
114: {
115: PetscInt dim, m, n, p, dof, swidth;
116: DMDAStencilType stencil;
117: DMBoundaryType bx, by, bz;
118: PetscBool coors;
119: DM dac;
120: Vec c;
122: PetscFunctionBegin;
123: PetscCall(PetscViewerBinaryRead(viewer, &dim, 1, NULL, PETSC_INT));
124: PetscCall(PetscViewerBinaryRead(viewer, &m, 1, NULL, PETSC_INT));
125: PetscCall(PetscViewerBinaryRead(viewer, &n, 1, NULL, PETSC_INT));
126: PetscCall(PetscViewerBinaryRead(viewer, &p, 1, NULL, PETSC_INT));
127: PetscCall(PetscViewerBinaryRead(viewer, &dof, 1, NULL, PETSC_INT));
128: PetscCall(PetscViewerBinaryRead(viewer, &swidth, 1, NULL, PETSC_INT));
129: PetscCall(PetscViewerBinaryRead(viewer, &bx, 1, NULL, PETSC_ENUM));
130: PetscCall(PetscViewerBinaryRead(viewer, &by, 1, NULL, PETSC_ENUM));
131: PetscCall(PetscViewerBinaryRead(viewer, &bz, 1, NULL, PETSC_ENUM));
132: PetscCall(PetscViewerBinaryRead(viewer, &stencil, 1, NULL, PETSC_ENUM));
134: PetscCall(DMSetDimension(da, dim));
135: PetscCall(DMDASetSizes(da, m, n, p));
136: PetscCall(DMDASetBoundaryType(da, bx, by, bz));
137: PetscCall(DMDASetDof(da, dof));
138: PetscCall(DMDASetStencilType(da, stencil));
139: PetscCall(DMDASetStencilWidth(da, swidth));
140: PetscCall(DMSetUp(da));
141: PetscCall(PetscViewerBinaryRead(viewer, &coors, 1, NULL, PETSC_ENUM));
142: if (coors) {
143: PetscCall(DMGetCoordinateDM(da, &dac));
144: PetscCall(DMCreateGlobalVector(dac, &c));
145: PetscCall(VecLoad(c, viewer));
146: PetscCall(DMSetCoordinates(da, c));
147: PetscCall(VecDestroy(&c));
148: }
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
153: {
154: DM_DA *da = (DM_DA *)dm->data;
156: PetscFunctionBegin;
157: if (subdm) {
158: PetscSF sf;
159: Vec coords;
160: void *ctx;
161: /* Cannot use DMClone since the dof stuff is mixed in. Ugh
162: PetscCall(DMClone(dm, subdm)); */
163: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), subdm));
164: PetscCall(DMGetPointSF(dm, &sf));
165: PetscCall(DMSetPointSF(*subdm, sf));
166: PetscCall(DMGetApplicationContext(dm, &ctx));
167: PetscCall(DMSetApplicationContext(*subdm, ctx));
168: PetscCall(DMGetCoordinatesLocal(dm, &coords));
169: if (coords) {
170: PetscCall(DMSetCoordinatesLocal(*subdm, coords));
171: } else {
172: PetscCall(DMGetCoordinates(dm, &coords));
173: if (coords) PetscCall(DMSetCoordinates(*subdm, coords));
174: }
176: PetscCall(DMSetType(*subdm, DMDA));
177: PetscCall(DMSetDimension(*subdm, dm->dim));
178: PetscCall(DMDASetSizes(*subdm, da->M, da->N, da->P));
179: PetscCall(DMDASetNumProcs(*subdm, da->m, da->n, da->p));
180: PetscCall(DMDASetBoundaryType(*subdm, da->bx, da->by, da->bz));
181: PetscCall(DMDASetDof(*subdm, numFields));
182: PetscCall(DMDASetStencilType(*subdm, da->stencil_type));
183: PetscCall(DMDASetStencilWidth(*subdm, da->s));
184: PetscCall(DMDASetOwnershipRanges(*subdm, da->lx, da->ly, da->lz));
185: }
186: if (is) {
187: PetscInt *indices, cnt = 0, dof = da->w, i, j;
189: PetscCall(PetscMalloc1(da->Nlocal * numFields / dof, &indices));
190: for (i = da->base / dof; i < (da->base + da->Nlocal) / dof; ++i) {
191: for (j = 0; j < numFields; ++j) indices[cnt++] = dof * i + fields[j];
192: }
193: PetscCheck(cnt == da->Nlocal * numFields / dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count %" PetscInt_FMT " does not equal expected value %" PetscInt_FMT, cnt, da->Nlocal * numFields / dof);
194: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), cnt, indices, PETSC_OWN_POINTER, is));
195: }
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: static PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
200: {
201: PetscInt i;
202: DM_DA *dd = (DM_DA *)dm->data;
203: PetscInt dof = dd->w;
205: PetscFunctionBegin;
206: if (len) *len = dof;
207: if (islist) {
208: Vec v;
209: PetscInt rstart, n;
211: PetscCall(DMGetGlobalVector(dm, &v));
212: PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
213: PetscCall(VecGetLocalSize(v, &n));
214: PetscCall(DMRestoreGlobalVector(dm, &v));
215: PetscCall(PetscMalloc1(dof, islist));
216: for (i = 0; i < dof; i++) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)dm), n / dof, rstart + i, dof, &(*islist)[i]));
217: }
218: if (namelist) {
219: PetscCall(PetscMalloc1(dof, namelist));
220: if (dd->fieldname) {
221: for (i = 0; i < dof; i++) PetscCall(PetscStrallocpy(dd->fieldname[i], &(*namelist)[i]));
222: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Currently DMDA must have fieldnames");
223: }
224: if (dmlist) {
225: DM da;
227: PetscCall(DMDACreate(PetscObjectComm((PetscObject)dm), &da));
228: PetscCall(DMSetDimension(da, dm->dim));
229: PetscCall(DMDASetSizes(da, dd->M, dd->N, dd->P));
230: PetscCall(DMDASetNumProcs(da, dd->m, dd->n, dd->p));
231: PetscCall(DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz));
232: PetscCall(DMDASetDof(da, 1));
233: PetscCall(DMDASetStencilType(da, dd->stencil_type));
234: PetscCall(DMDASetStencilWidth(da, dd->s));
235: PetscCall(DMSetUp(da));
236: PetscCall(PetscMalloc1(dof, dmlist));
237: for (i = 0; i < dof - 1; i++) PetscCall(PetscObjectReference((PetscObject)da));
238: for (i = 0; i < dof; i++) (*dmlist)[i] = da;
239: }
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode DMClone_DA(DM dm, DM *newdm)
244: {
245: DM_DA *da = (DM_DA *)dm->data;
247: PetscFunctionBegin;
248: PetscCall(DMSetType(*newdm, DMDA));
249: PetscCall(DMSetDimension(*newdm, dm->dim));
250: PetscCall(DMDASetSizes(*newdm, da->M, da->N, da->P));
251: PetscCall(DMDASetNumProcs(*newdm, da->m, da->n, da->p));
252: PetscCall(DMDASetBoundaryType(*newdm, da->bx, da->by, da->bz));
253: PetscCall(DMDASetDof(*newdm, da->w));
254: PetscCall(DMDASetStencilType(*newdm, da->stencil_type));
255: PetscCall(DMDASetStencilWidth(*newdm, da->s));
256: PetscCall(DMDASetOwnershipRanges(*newdm, da->lx, da->ly, da->lz));
257: PetscCall(DMSetUp(*newdm));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: static PetscErrorCode DMHasCreateInjection_DA(DM dm, PetscBool *flg)
262: {
263: DM_DA *da = (DM_DA *)dm->data;
265: PetscFunctionBegin;
267: PetscAssertPointer(flg, 2);
268: *flg = da->interptype == DMDA_Q1 ? PETSC_TRUE : PETSC_FALSE;
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode DMGetDimPoints_DA(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
273: {
274: PetscFunctionBegin;
275: PetscCall(DMDAGetDepthStratum(dm, dim, pStart, pEnd));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: static PetscErrorCode DMGetNeighbors_DA(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[])
280: {
281: PetscInt dim;
282: DMDAStencilType st;
284: PetscFunctionBegin;
285: PetscCall(DMDAGetNeighbors(dm, ranks));
286: PetscCall(DMDAGetInfo(dm, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &st));
288: switch (dim) {
289: case 1:
290: *nranks = 3;
291: /* if (st == DMDA_STENCIL_STAR) *nranks = 3; */
292: break;
293: case 2:
294: *nranks = 9;
295: /* if (st == DMDA_STENCIL_STAR) *nranks = 5; */
296: break;
297: case 3:
298: *nranks = 27;
299: /* if (st == DMDA_STENCIL_STAR) *nranks = 7; */
300: break;
301: default:
302: break;
303: }
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: /*MC
308: DMDA = "da" - A `DM` object that is used to manage data for a structured grid in 1, 2, or 3 dimensions.
309: In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points.
310: In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width set with `DMDASetStencilWidth()`.
312: The vectors can be thought of as either cell centered or vertex centered on the mesh. But some variables cannot be cell centered and others
313: vertex centered; see the documentation for `DMSTAG`, a similar `DM` implementation which supports more general staggered grids.
315: Level: intermediate
317: .seealso: [](sec_struct), `DMType`, `DMCOMPOSITE`, `DMSTAG`, `DMDACreate()`, `DMCreate()`, `DMSetType()`, `DMDASetStencilWidth()`, `DMDASetStencilType()`,
318: `DMDAStencilType`
319: M*/
321: PETSC_EXTERN PetscErrorCode DMCreate_DA(DM da)
322: {
323: DM_DA *dd;
325: PetscFunctionBegin;
326: PetscAssertPointer(da, 1);
327: PetscCall(PetscNew(&dd));
328: da->data = dd;
330: da->dim = -1;
331: dd->interptype = DMDA_Q1;
332: dd->refine_x = 2;
333: dd->refine_y = 2;
334: dd->refine_z = 2;
335: dd->coarsen_x = 2;
336: dd->coarsen_y = 2;
337: dd->coarsen_z = 2;
338: dd->fieldname = NULL;
339: dd->nlocal = -1;
340: dd->Nlocal = -1;
341: dd->M = -1;
342: dd->N = -1;
343: dd->P = -1;
344: dd->m = -1;
345: dd->n = -1;
346: dd->p = -1;
347: dd->w = -1;
348: dd->s = -1;
350: dd->xs = -1;
351: dd->xe = -1;
352: dd->ys = -1;
353: dd->ye = -1;
354: dd->zs = -1;
355: dd->ze = -1;
356: dd->Xs = -1;
357: dd->Xe = -1;
358: dd->Ys = -1;
359: dd->Ye = -1;
360: dd->Zs = -1;
361: dd->Ze = -1;
363: dd->Nsub = 1;
364: dd->xol = 0;
365: dd->yol = 0;
366: dd->zol = 0;
367: dd->xo = 0;
368: dd->yo = 0;
369: dd->zo = 0;
370: dd->Mo = -1;
371: dd->No = -1;
372: dd->Po = -1;
374: dd->gtol = NULL;
375: dd->ltol = NULL;
376: dd->ao = NULL;
377: PetscCall(PetscStrallocpy(AOBASIC, (char **)&dd->aotype));
378: dd->base = -1;
379: dd->bx = DM_BOUNDARY_NONE;
380: dd->by = DM_BOUNDARY_NONE;
381: dd->bz = DM_BOUNDARY_NONE;
382: dd->stencil_type = DMDA_STENCIL_BOX;
383: dd->interptype = DMDA_Q1;
384: dd->lx = NULL;
385: dd->ly = NULL;
386: dd->lz = NULL;
388: dd->elementtype = DMDA_ELEMENT_Q1;
390: da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA;
391: da->ops->globaltolocalend = DMGlobalToLocalEnd_DA;
392: da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA;
393: da->ops->localtoglobalend = DMLocalToGlobalEnd_DA;
394: da->ops->localtolocalbegin = DMLocalToLocalBegin_DA;
395: da->ops->localtolocalend = DMLocalToLocalEnd_DA;
396: da->ops->createglobalvector = DMCreateGlobalVector_DA;
397: da->ops->createlocalvector = DMCreateLocalVector_DA;
398: da->ops->createinterpolation = DMCreateInterpolation_DA;
399: da->ops->getcoloring = DMCreateColoring_DA;
400: da->ops->creatematrix = DMCreateMatrix_DA;
401: da->ops->refine = DMRefine_DA;
402: da->ops->coarsen = DMCoarsen_DA;
403: da->ops->refinehierarchy = DMRefineHierarchy_DA;
404: da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA;
405: da->ops->createinjection = DMCreateInjection_DA;
406: da->ops->hascreateinjection = DMHasCreateInjection_DA;
407: da->ops->destroy = DMDestroy_DA;
408: da->ops->view = NULL;
409: da->ops->setfromoptions = DMSetFromOptions_DA;
410: da->ops->setup = DMSetUp_DA;
411: da->ops->clone = DMClone_DA;
412: da->ops->load = DMLoad_DA;
413: da->ops->createcoordinatedm = DMCreateCoordinateDM_DA;
414: da->ops->createsubdm = DMCreateSubDM_DA;
415: da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA;
416: da->ops->createdomaindecomposition = DMCreateDomainDecomposition_DA;
417: da->ops->createddscatters = DMCreateDomainDecompositionScatters_DA;
418: da->ops->getdimpoints = DMGetDimPoints_DA;
419: da->ops->getneighbors = DMGetNeighbors_DA;
420: da->ops->locatepoints = DMLocatePoints_DA_Regular;
421: da->ops->getcompatibility = DMGetCompatibility_DA;
422: PetscCall(PetscObjectComposeFunction((PetscObject)da, "DMSetUpGLVisViewer_C", DMSetUpGLVisViewer_DMDA));
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: /*@
427: DMDACreate - Creates a `DMDA` object for managing structured grids.
429: Collective
431: Input Parameter:
432: . comm - The communicator for the `DMDA` object
434: Output Parameter:
435: . da - The DMDA object
437: Level: advanced
439: Developer Note:
440: Since there exists DMDACreate1/2/3d() should this routine even exist?
442: .seealso: [](sec_struct), `DMDASetSizes()`, `DMClone()`, `DMDACreate1d()`, `DMDACreate2d()`, `DMDACreate3d()`
443: @*/
444: PetscErrorCode DMDACreate(MPI_Comm comm, DM *da)
445: {
446: PetscFunctionBegin;
447: PetscAssertPointer(da, 2);
448: PetscCall(DMCreate(comm, da));
449: PetscCall(DMSetType(*da, DMDA));
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }