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