1(ch_stag)= 2 3# DMSTAG: Staggered, Structured Grid 4 5For structured (aka "regular") grids with staggered data (degrees of freedom potentially living on elements, faces, edges, 6and/or vertices), the `DMSTAG` object is available. This can 7be useful for problems in many domains, including fluid flow, MHD, and seismology. 8 9It is possible, though extremely cumbersome, to implement a staggered-grid code using multiple `DMDA` objects, or a single 10multi-component `DMDA` object where some degrees of freedom are unused. 11`DMSTAG` was developed for two main purposes: 12 131. To help manage some of the burden of choosing and adhering to the complex indexing conventions needed for staggered grids (in parallel) 142. To provide a uniform abstraction for which scalable solvers and preconditioners may be developed (in particular, using `PCFIELDSPLIT` and `PCMG`). 15 16`DMSTAG` is design to behave much 17like {ref}`DMDA <sec_struct>`, with a couple of important distinctions, and borrows some terminology 18from {doc}`DMPLEX <dmplex>`. 19 20## Terminology 21 22Like a `DMPLEX` object, a `DMSTAG` represents a [cell complex](https://en.wikipedia.org/wiki/CW_complex), 23distributed in parallel over the ranks of an `MPI_Comm`. It is, however, 24a very regular complex, consisting of a structured grid of $d$-dimensional cells, with $d \in \{1,2,3\}$, 25which are referred to as *elements*, $d-1$ dimensional cells defining boundaries between these elements, 26and the boundaries of the domain, and in 2 or more dimensions, boundaries of *these* cells, 27all the way down to 0 dimensional cells referred to as *vertices*. In 2 dimensions, the 1-dimensional 28element boundaries are referred to as *edges* or *faces*. In 3 dimensions, the 2-dimensional element boundaries 29are referred to as *faces* and the 1-dimensional boundaries between faces are referred to as *edges* 30The set of cells of a given dimension is referred to as a *stratum* (which one can think of as a level in DAG representation of the mesh); a `DMSTAG` object of dimension $d$ 31represents a complete cell complex with $d+1$ *strata* (levels). 32 33In the description of any:`ch_unstructured` the cells at each level are referred to as *points*. Thus we adopt that terminology uniformly in PETSc and 34so furthermore in this document, 35point will refer to a cell. 36 37Each stratum has a constant number of unknowns (which may be zero) associated with each point (cell) on that level. 38The distinct unknowns associated with each point are referred to as *components*. For a `DMPLEX` there may be a different 39number of unknowns for each point on the same level. 40 41The structured grid, is like with `DMDA`, decomposed via a Cartesian product of decompositions in each dimension, 42giving a rectangular local subdomain on each rank. This is extended by an element-wise stencil width 43of *ghost* elements to create an atlas of overlapping patches. 44 45## Working with vectors and operators (matrices) 46 47`DMSTAG` allows the user to reason almost entirely about a global indexing of elements. 48Element indices are simply 1-3 `PetscInt` values, starting at $0$, in the 49back, bottom, left corner of the domain. For instance, element $(1,2,3)$, in 3D, 50is the element second from the left, third from the bottom, and fourth from the back 51(regardless of how many MPI ranks are used). 52 53To refer to points (elements, faces, edges, and vertices), a value of 54`DMStagStencilLocation` is used, relative to the element index. The element 55itself is referred to with `DMSTAG_ELEMENT`, the top right vertex (in 2D) 56or the top right edge (in 3D) with `DMSTAG_UP_RIGHT`, the back bottom left 57corner in 3D with `DMSTAG_BACK_DOWN_LEFT`, and so on. 58 59{numref}`figure_dmstag_indexing` gives a few examples in 2D. 60 61:::{figure} /images/manual/dmstag_indexing.svg 62:name: figure_dmstag_indexing 63 64Locations in `DMSTAG` are indexed according to global element indices (here, two in 2D) and a location name. Elements have unique names but other locations can be referred to in more than one way. Element colors correspond to a parallel decomposition, but locations on the grid have names which are invariant to this. Note that the face on the top right can be referred to as being to the left of a "dummy" element $(3,3)$ outside the physical domain. 65::: 66 67Crucially, this global indexing scheme does not include any "ghost" or "padding" unknowns outside the physical domain. 68This is useful for higher-level operations such as computing norms or developing physics-based solvers. However 69(unlike `DMDA`), this implies that the global `Vec` do not have a natural block structure, as different 70strata have different numbers of points (e.g. in 1D there is an "extra" vertex on the right). This regular block 71structure is, however, very useful for the *local* representation of the data, so in that case *dummy* DOF 72are included, drawn as grey in {numref}`figure_dmstag_local_global`. 73 74:::{figure} /images/manual/dmstag_local_global.svg 75:name: figure_dmstag_local_global 76 77Local and global representations for a 2D `DMSTAG` object, 3 by 4 elements, with one degree of freedom on each of the three strata: element (squares), faces (triangles), and vertices (circles). The cell complex is parallelized across 4 MPI ranks. In the global representation, the colors correspond to which rank holds the native representation of the unknown. The 4 local representations are shown, with an (elementwise) stencil "box" stencil width of 1. Unknowns are colored by their native rank. Dummy unknowns, which correspond to no global degree of freedom, are colored grey. Note that the local representations have a natural block size of 4, and the global representation has no natural block size. 78::: 79 80For working with `Vec` data, this approach is used to allow direct access to a multi-dimensional, regular-blocked 81array. To avoid the user having to know about the {ref}`internal numbering conventions used <sec_dmstag_numbering>`, 82helper functions are used to produce the proper final integer index for a given location and component, referred to as a "slot". 83Similarly to `DMDAVecGetArrayDOF()`, this uses a $d+1$ dimensional array in $d$ dimensions. 84The following snippet give an example of this usage. 85 86```{literalinclude} /../src/dm/impls/stag/tests/ex51.c 87:end-at: PetscCall(DMStagVecRestoreArray 88:start-at: /* Set 89``` 90 91`DMSTAG` provides a stencil-based method for getting and setting entries of `Mat` and `Vec` objects. 92The follow excerpt from <a href="PETSC_DOC_OUT_ROOT_PLACEHOLDER/src/dm/impls/stag/tutorials/ex1.c.html">DMSTAG Tutorial ex1</a> demonstrates 93the idea. For more, see the manual page for `DMStagMatSetValuesStencil()`. 94 95```{literalinclude} /../src/dm/impls/stag/tutorials/ex1.c 96:end-at: PetscCall(DMStagMatSetValuesStencil 97:start-at: /* Velocity is either a BC 98``` 99 100The array-based approach for `Vec` is likely to be more efficient than the stencil-based method just introduced above. 101 102## Coordinates 103 104`DMSTAG`, unlike `DMDA`, supports two approaches to defining coordinates. This is captured by which type of `DM` 105is used to represent the coordinates. No default is imposed, so the user must directly or indirectly call 106`DMStagSetCoordinateDMType()`. 107 108If a second `DMSTAG` object is used to represent coordinates in "explicit" form, behavior is much like with `DMDA` - the coordinate `DM` 109has $d$ DOF on each stratum corresponding to coordinates associated with each point. 110 111If `DMPRODUCT` is used instead, coordinates are represented by a `DMPRODUCT` object referring to a 112Cartesian product of 1D `DMSTAG` objects, each of which features explicit coordinates as just mentioned. 113 114Navigating these nested `DM` in `DMPRODUCT` can be tedious, but note the existence of helper functions like 115`DMStagSetUniformCoordinatesProduct()` and `DMStagGetProductCoordinateArrays()`. 116 117(sec_dmstag_numbering)= 118 119## Numberings and internal data layout 120 121While `DMSTAG` aims to hide the details of its internal data layout, for debugging, optimization, and 122customization purposes, it can be important to know how `DMSTAG` internally numbers unknowns. 123 124Internally, each point is canonically associated with an element (top-level point (cell)). For purposes of local, 125regular-blocked storage, an element is grouped with lower-dimensional points left of, below ("down"), and behind ("back") it. 126This means that "canonical" values of `DMStagStencilLocation` are `DMSTAG_ELEMENT`, plus all entries consisting only of "LEFT", "DOWN", and "BACK". In general, these are the most efficient values to use, unless convenience dictates otherwise, as they are the ones used internally. 127 128When creating the decomposition of the domain to local ranks, and extending these local domains to handle overlapping halo regions and boundary ghost unknowns, this same per-element association is used. This has the advantage of maintaining a regular blocking, but may not be optimal in some situations in terms of data movement. 129 130Numberings are, like `DMDA`, based on a local "x-fastest, z-slowest" or "PETSc" ordering of elements (see {ref}`sec_ao`), with ordering of locations canonically associated with each element decided by considering unknowns on each point 131to be located at the center of their point, and using a nested ordering of the same style. Thus, in 3-D, the ordering of the 8 canonical `DMStagStencilLocation` values associated 132with an element is 133 134``` 135DMSTAG_BACK_DOWN_LEFT 136DMSTAG_BACK_DOWN 137DMSTAG_BACK_LEFT 138DMSTAG_BACK 139DMSTAG_DOWN_LEFT 140DMSTAG_DOWN 141DMSTAG_LEFT 142DMSTAG_ELEMENT 143``` 144 145Multiple DOF associated with a given point are stored sequentially (as with `DMDA`). 146 147For local `Vec`s, this gives a regular-blocked numbering, with the same number of unknowns associated with each element, including some "dummy" unknowns which to not correspond to any (local or global) unknown in the global representation. See {numref}`figure_dmstag_numbering_local` for an example. 148 149In the global representation, only physical unknowns are numbered (using the same "Z" ordering for unknowns which are present), giving irregular numbers of unknowns, depending on whether a domain boundary is present. See {numref}`figure_dmstag_numbering_global` for an example. 150 151:::{figure} /images/manual/dmstag_numbering_global.svg 152:name: figure_dmstag_numbering_global 153 154Global numbering scheme for a 2D `DMSTAG` object with one DOF per stratum. Note that the numbering depends on the parallel decomposition (over 4 ranks, here). 155::: 156 157:::{figure} /images/manual/dmstag_numbering_local.svg 158:name: figure_dmstag_numbering_local 159 160Local numbering scheme on rank 1 (Cf. {numref}`figure_dmstag_local_global`) for a 2D `DMSTAG` object with one DOF per stratum. Note that dummy locations (grey) are used to give a regular block size (here, 4). 161::: 162 163It should be noted that this is an *interlaced* (AoS) representation. If a segregated (SoA) representation is required, 164one should use `DMCOMPOSITE` collecting several `DMSTAG` objects, perhaps using `DMStagCreateCompatibleDMStag()` to 165quickly create additional `DMSTAG` objects from an initial one. 166
