spatial-discretization

name: spatial-discretization description: Reference for DMStag grid layout and FlucaFD finite difference operators. Use when working on spatial discretization, stencils, operator composition, or understanding the collocated grid.

Spatial Discretization: DMStag + FlucaFD

DMStag Grid Layout

Fluca uses PETSc's DMStag for structured grid storage. DOFs live at vertices, edges, faces, and cell centers.

Stencil Locations

Fluca uses "lower-left" convention — only LEFT/DOWN/BACK faces, never RIGHT/UP/FRONT:

Face count: N+1 (non-periodic) or N (periodic) per direction. Cell centers: always N.

Coordinates

DMStag product coordinates with two slots per direction:

• DMSTAG_LEFT slot → face coordinate x_{i-1/2}
• DMSTAG_ELEMENT slot → cell center x_i

Access: DMStagGetProductCoordinateArraysRead() + DMStagGetProductCoordinateLocationSlot().

Global and Local Vectors

DMStag manages two vector types that differ in size and ghost point inclusion.

entriesPerElement packs all DOFs at the "lower-left" corner of one element (vertex + edges + faces + element center). For example in 2D with 1 DOF per stratum: entriesPerElement = 4 (vertex + left-face + down-face + element).

Global vector — distributed across all ranks in the MPI communicator, one entry per owned DOF, no ghost points. Uses global numbering (indices are unique across all ranks). Size is not simply owned_elements × entriesPerElement when non-element strata have DOFs and boundaries are non-periodic. The last rank in each direction owns extra DOFs for the upper boundary (right faces, top faces, far-right vertex, etc.) that are not covered by entriesPerElement.

Local vector — lives on a single process (PETSC_COMM_SELF), includes owned DOFs plus ghost region (neighbor data needed for stencil operations). Uses local numbering (always starts at 0, so index 0 refers to different physical elements on different ranks). Size = gxm * gym * gzm * entriesPerElement (the ghost region is always a full rectangle of elements, so entriesPerElement applies uniformly).

Distribution example — 2D, 2x2 elements, 2 ranks (split in x), stencil width 1, non-periodic, 1 DOF per stratum (entriesPerElement = 4: V vertex + D down-face + L left-face + E element):

Global vector (numbers are global indices):

    Rank 0    ||         Rank 1
              ||
  V8 -- D9 ---++--- V22 -- D23 -- V24
  |           ||     |             |
  L6    E7    ||    L18    E19    L21
  |           ||     |             |
  V4 -- D5 ---++--- V16 -- D17 -- V20
  |           ||     |             |
  L2    E3    ||    L12    E13    L15
  |           ||     |             |
  V0 -- D1 ---++--- V10 -- D11 -- V14
              ||
  10 entries  ||       15 entries
  (1x2x4 + 2) ||       (1x2x4 + 7)

Extra boundary DOFs on the last rank per direction: vertices + the face stratum perpendicular to the boundary. At x=2 (right): V + L (left-face = right boundary face). At y=2 (top): V + D (down-face = top boundary face). At (2,2): V only.

Local vectors (numbers are local indices):

Rank 0 (gxs=0, gxm=2, gys=0, gym=3):     24 = 2 x 3 x 4

  L18    E19    L22    E23      ghost(padding)
   |             |
  V16 -- D17 -- V20 -- D21    ^
   |             |            |
  L10    E11    L14    E15    |
   |             |            |
   V8 --- D9 -- V12 -- D13    | owned
   |             |            |
   L2     E3     L6     E7    |
   |             |            |
   V0 --- D1 --- V4 --- D5    v

   <------>      <------>
    owned         ghost

Rank 1 (gxs=0, gxm=3, gys=0, gym=3):     36 = 3 x 3 x 4

  L26    E27    L30    E31    L34    E35      ghost(padding)
   |             |             |
  V24 -- D25 -- V28 -- D29 -- V32 -- D33    ^
   |             |             |            |
  L14    E15    L18    E19    L22    E23    |
   |             |             |            |
  V12 -- D13 -- V16 -- D17 -- V20 -- D21    | owned
   |             |             |            |
   L2     E3     L6     E7    L10    E11    |
   |             |             |            |
   V0 --- D1 --- V4 --- D5 --- V8 --- D9    v

   <------>      <------------->
    ghost             owned         ghost
                                  (padding)

Each element occupies 4 sequential slots [V, D, L, E]. At boundary phantom elements, only boundary-facing slots carry real data (V+D at top, V+L at right, V only at corner). The remaining slots are padding — they have no physical meaning and exist solely to keep the local vector size exactly gxm * gym * gzm * entriesPerElement, so that arr[j][i][slot] indexing works uniformly without special-casing boundary elements.

Typical workflow:

/* 1. Get a temporary local vector (includes ghost region) */
DMGetLocalVector(dm, &local);

/* 2. Scatter global → local, filling ghost points from neighbor ranks */
DMGlobalToLocal(dm, global, INSERT_VALUES, local);

/* 3. Access as multi-dimensional array: arr[k][j][i][slot] */
DMStagVecGetArrayRead(dm, local, &arr);

/* 4. Loop over owned region (DMStagGetCorners gives start + count + extra) */
DMStagGetCorners(dm, &xs, &ys, &zs, &xm, &ym, &zm, &xe, &ye, &ze);
for (j = ys; j < ys + ym + ye; j++)
  for (i = xs; i < xs + xm + xe; i++)
    val = arr[j][i][slot];  /* ghost points accessible via neighbor indices */

/* 5. Restore and return */
DMStagVecRestoreArrayRead(dm, local, &arr);
DMRestoreLocalVector(dm, &local);

Key functions:

Corner variable naming:

Array indexing — DMStagVecGetArray returns an array indexed by global element indices (not zero-based local). In 2D: arr[j][i][slot] where i,j range over ghost corners (gxs..gxs+gxm-1). The last dimension (slot) packs all DOFs within one element. Use DMStagGetLocationSlot() to find which slot corresponds to a given (location, component).

xe, ye, ze — 0 or 1, indicating whether the local rank owns the extra boundary face/vertex in that direction. Loop bounds for face-located DOFs extend by these values.

IMPORTANT: DMStagVecGetArray and DMStagVecGetArrayRead require a local vector. Never pass a global vector to these functions — it will silently produce wrong index mappings because global vectors lack the uniform gxm × gym × gzm × entriesPerElement layout that the array view assumes. Always scatter to a local vector first (DMGlobalToLocal) or obtain one via DMGetLocalVector.

FlucaFD Operator System

Polymorphic FD operators that compute stencils mapping input DOFs to output DOFs on DMStag.

Subtypes

Pre-computed subtypes (Derivative, SecondOrderTVD) build stencil coefficients at SetUp time and store them per grid point. Lazy subtypes (Composition, Scale, Sum) have no SetUp stencil storage — they compute stencils on-the-fly by recursively calling GetStencilRaw on their children and combining the results.

Lifecycle

Create → [SetType] → SetFromOptions → SetBoundaryConditions → SetUp → Apply/GetOperator → Destroy

Convenience constructors (FlucaFDDerivativeCreate, etc.) handle Create + SetType + config.

Stencil Representation

A stencil = array of (DMStagStencil col, PetscScalar v) pairs. The col.c field encodes:

• >= 0: Interior grid point (component index)
• FLUCAFD_BOUNDARY_LEFT (-1) through _FRONT (-6): Boundary value marker
• FLUCAFD_CONSTANT (-7): Additive constant term

Core API

FlucaFDGetStencil(fd, i, j, k, &n, col[], v[])   /* Stencil with BC handling */
FlucaFDApply(fd, input_dm, output_dm, x, y)        /* y = fd(x) */
FlucaFDGetOperator(fd, input_dm, output_dm, mat)    /* Sparse matrix (interior only) */

Apply Internal Flow

FlucaFDApply(fd, input_dm, output_dm, x, y) computes y = fd(x):

1. Scatter input: DMGlobalToLocal(input_dm, x, INSERT_VALUES, x_local) — fills ghost points
2. Zero local output: VecZeroEntries(y_local)
3. Loop over owned output cells (range from DMStagGetCorners + extra):
   • Call FlucaFDGetStencil(fd, i, j, k, &ncols, col, v) — returns stencil with BCs resolved
   • Dot product: for each stencil entry, col[c].c >= 0 reads from x_local, col[c].c == FLUCAFD_CONSTANT adds constant, otherwise applies bcs[].value
   • Write result to y_local
4. Scatter output: DMLocalToGlobal(output_dm, y_local, INSERT_VALUES, y)

FlucaFDGetOperator(fd, input_dm, output_dm, mat) assembles the sparse matrix by the same loop but only inserts interior stencil points (col[c].c >= 0) via MatSetValuesLocal with ADD_VALUES. Boundary and constant terms are excluded. Caller must call MatAssemblyBegin/End afterward.

Boundary Conditions

FlucaFDBoundaryCondition bcs[2 * dim];  /* [left, right, down, up, back, front] */
bcs[0].type = FLUCAFD_BC_DIRICHLET;    /* or _NEUMANN, _NONE */
bcs[0].value = 0.0;
FlucaFDSetBoundaryConditions(fd, bcs);

Off-grid stencil points are resolved by FlucaFDRemoveOffGridPoints_Internal():

• BC_NONE: Extrapolate from interior (Vandermonde)
• BC_DIRICHLET: Include boundary value in polynomial
• BC_NEUMANN: Infer from prescribed derivative

Operator Composition Patterns

Derivative (element ↔ face transition):

FlucaFDDerivativeCreate(dm, FLUCAFD_X, 1, 2,
    DMSTAG_ELEMENT, 0, DMSTAG_LEFT, 0, &ddx);  /* d/dx: elem→face */

Convection (ρu · ∇φ):

TVD(ELEMENT→face) → Scale(ρ) → d/dx(face→ELEMENT)

FlucaFDSecondOrderTVDCreate(dm, FLUCAFD_X, 0, 0, &tvd);
FlucaFDScaleCreateConstant(tvd, rho, &scaled);
FlucaFDDerivativeCreate(dm, FLUCAFD_X, 1, 2, DMSTAG_LEFT, 0, DMSTAG_ELEMENT, 0, &ddx);
FlucaFDCompositionCreate(scaled, ddx, &conv_x);

Diffusion (∇·(μ∇φ)):

d/dx(ELEMENT→face) → Scale(μ) → d/dx(face→ELEMENT)

FlucaFDDerivativeCreate(dm, FLUCAFD_X, 1, 2, DMSTAG_ELEMENT, 0, DMSTAG_LEFT, 0, &inner);
FlucaFDScaleCreateConstant(inner, mu, &scaled);
FlucaFDDerivativeCreate(dm, FLUCAFD_X, 1, 2, DMSTAG_LEFT, 0, DMSTAG_ELEMENT, 0, &outer);
FlucaFDCompositionCreate(scaled, outer, &diff_x);

Multi-direction sum (Laplacian, full convection-diffusion):

FlucaFD ops[] = {conv_x, conv_y, neg_diff_x, neg_diff_y};
FlucaFDSumCreate(4, ops, &full);

SecondOrderTVD Usage

Nonlinear — stencil depends on solution. Must update before each Apply:

FlucaFDSecondOrderTVDSetVelocity(tvd, vel, vel_c);
FlucaFDSecondOrderTVDSetCurrentSolution(tvd, phi);

11 limiters: superbee, minmod, MC, vanLeer, vanAlbada, Barth-Jespersen, Venkatakrishnan, Koren, upwind, SOU, QUICK.

Options Prefixes

Multiple operators share the same option names. Use prefixes to disambiguate:

FlucaFDSetOptionsPrefix(fd_x, "x_");  /* -x_flucafd_deriv_order */
FlucaFDSetOptionsPrefix(fd_y, "y_");  /* -y_flucafd_deriv_order */

Key Source Files