Unified Stencil Operations for Lattice Field Theory

This module provides a single, unified stencil type that works identically across all backends (OpenMP, OpenCL, SYCL) and handles distributed lattices with ghost regions automatically.

Key Design Principles:

One Type: LatticeStencil[D] is the only stencil type you need
Always Distributed: Lattices are always distributed, so stencils always understand ghost regions and padded layouts
Backend Agnostic: Same API inside each loops for all backends
Clean Operators: Use view[stencil.shift(n, dir)] for neighbor access

Usage Example:

# Create stencil from lattice - D is inferred automatically
let stencil = newLatticeStencil(nearestNeighborStencil(lat), lat)
# Or even simpler (defaults to nearest-neighbor):
let stencil = newLatticeStencil(lat)

# Inside any backend's each loop:
for n in each 0..<view.numSites:
  for dir in 0..<4:
    # Forward neighbor
    let fwd = view[stencil.shift(n, +1, dir)]
    # Backward neighbor
    let bwd = view[stencil.shift(n, -1, dir)]
    result[n] = fwd + bwd - 2.0 * view[n]

The stencil handles:

Periodic boundary conditions within local domain
Ghost region detection and padded coordinate mapping
SIMD-aware indexing for vectorized backends
Device buffer offsets for GPU backends

Imports

latticeconcept, ../globalarrays/gatypes

Types

Direction = distinct int

Type-safe direction index (0 = x, 1 = y, 2 = z, 3 = t for 4D)

LatticeStencil[D] = object
  pattern*: StencilPattern[D]
  localGeom*: array[D, int]
  ghostWidth*: array[D, int]
  paddedGeom*: array[D, int]
  nLocalSites*: int
  nPaddedSites*: int
  nPoints*: int
  backend*: StencilBackend
  offsets*: seq[int32]
  isGhost*: seq[bool]
  simdGrid*: array[D, int]
  nSitesInner*: int
  nSitesOuter*: int
  nbrOuter*: seq[int32]
  nbrLane*: seq[int32]

Unified stencil for all backends and distributed lattices

Provides a single, clean API for neighbor access that works identically whether you're using OpenMP, OpenCL, or SYCL.

PathStep = object
  dir*: int
  sign*: int

SignedDirection = object
  dir*: Direction
  sign*: int

Direction with sign for forward (+1) or backward (-1) shifts

Stencil[D] = StencilPattern[D]

StencilBackend = enum
  sbScalar,                 ## Scalar CPU (no SIMD)
  sbSimd,                   ## SIMD-vectorized CPU (OpenMP)
  sbOpenCL,                 ## OpenCL GPU
  sbSycl                     ## SYCL GPU

StencilEntry = object
  offset*: int
  isLocal*: bool
  permute*: uint8
  wrapAround*: bool

StencilPattern[D] = object
  points*: seq[StencilPoint[D]]
  name*: string

A stencil pattern defining neighbor offsets

StencilPoint[D] = object
  offset*: array[D, int]
  id*: int

A single point in a stencil pattern (offset from center)

StencilShift[D] = object
  stencil*: ptr LatticeStencil[D]
  site*: int
  point*: int
  neighborIdx*: int

Result of stencil.shift() - used with view for neighbor access

StencilView[D] = object
  stencil*: StencilPattern[D]
  localGeom*: array[D, int]
  nSites*: int
  nPoints*: int
  entries*: seq[StencilEntry]

Consts

UseOpenCL {.booldefine.} = false
UseOpenMP {.booldefine.} = false
UseSYCL {.booldefine.} = false
VectorWidth {.intdefine.} = 1

Procs

proc `$`(d: Direction): string {....raises: [], tags: [], forbids: [].}
proc `$`(sd: SignedDirection): string {....raises: [], tags: [], forbids: [].}
proc `$`[D: static int](p: StencilPoint[D]): string
proc `$`[D: static int](s: StencilPattern[D]): string
proc `==`(a, b: Direction): bool {.borrow, ...raises: [], tags: [], forbids: [].}
proc addPoint[D: static int](s: var StencilPattern[D]; offset: array[D, int])
proc addPoint[D: static int](s: var StencilPattern[D]; sd: SignedDirection)
proc addPoint[D: static int](s: var StencilPattern[D]; sign: int; dir: int)
proc backward(d: Direction): SignedDirection {.inline, ...raises: [], tags: [], forbids: [].}
proc backwardStencil[D: static int; L: Lattice[D]](lat: L): StencilPattern[D]: Create a backward-only stencil, inferring D from the lattice
proc backwardStencil[D: static int](): StencilPattern[D]: Backward-only stencil (-x, -y, -z, -t)
proc bwd(dir: int): PathStep {....raises: [], tags: [], forbids: [].}
proc bwd[D: static int](s: LatticeStencil[D]; site: int; dir: int): StencilShift[ D] {.inline.}: Shorthand for backward shift
proc forward(d: Direction): SignedDirection {.inline, ...raises: [], tags: [], forbids: [].}
proc forwardStencil[D: static int; L: Lattice[D]](lat: L): StencilPattern[D]: Create a forward-only stencil, inferring D from the lattice
proc forwardStencil[D: static int](): StencilPattern[D]: Forward-only stencil (+x, +y, +z, +t)
proc fwd(dir: int): PathStep {....raises: [], tags: [], forbids: [].}
proc fwd[D: static int](s: LatticeStencil[D]; site: int; dir: int): StencilShift[ D] {.inline.}: Shorthand for forward shift
proc getEntry[D: static int](sv: StencilView[D]; site, point: int): StencilEntry {. inline.}
proc getOffsetBuffer[D: static int](s: LatticeStencil[D]): ptr int32 {.inline.}
proc hash(d: Direction): int {....raises: [], tags: [], forbids: [].}
proc idx(sh: StencilShift): int {.inline.}
proc isGhostNeighbor[D: static int](s: LatticeStencil[D]; site, point: int): bool {. inline.}: Check if neighbor is in a ghost region
proc laplacianStencil[D: static int; L: Lattice[D]](lat: L): StencilPattern[D]: Create a laplacian stencil, inferring D from the lattice
proc laplacianStencil[D: static int](): StencilPattern[D]
proc localToPadded[D: static int](localSite: int; localGeom: array[D, int]; ghostWidth: array[D, int]): int
proc localToPadded[D: static int](s: LatticeStencil[D]; localSite: int): int: Convert local site index to padded site index
proc nearestNeighborStencil[D: static int; L: Lattice[D]](lat: L): StencilPattern[ D]: Create a nearest-neighbor stencil, inferring D from the lattice

Example: let lat = newSimpleCubicLattice(8, 8, 8, 16) let pattern = nearestNeighborStencil(lat) # D=4 inferred
proc nearestNeighborStencil[D: static int](): StencilPattern[D]: Create a nearest-neighbor stencil with 2D points (forward and backward in each dimension).
proc neighbor[D: static int](s: LatticeStencil[D]; site, point: int): int {. inline.}: Get neighbor index in PADDED layout

This is the primary API for reading neighbors. The returned index is valid for accessing a padded tensor that includes ghost regions.
proc neighborOffset[D: static int](sv: StencilView[D]; site, point: int): int {. inline.}
proc neighborSimd[D: static int](s: LatticeStencil[D]; outer, lane, point: int): tuple[ outer, lane: int] {.inline.}: SIMD-aware neighbor lookup for vectorized backends
proc newLatticeStencil[D: static int; L: Lattice[D]](lat: L): LatticeStencil[D]: Create a nearest-neighbor stencil from a Lattice (most convenient constructor)

Infers D and uses nearestNeighborStencil automatically.

Example:

let lat = newSimpleCubicLattice([8, 8, 8, 16]) let stencil = newLatticeStencil(lat) # 4D nearest-neighbor stencil
proc newLatticeStencil[D: static int; L: Lattice[D]](pattern: StencilPattern[D]; lat: L): LatticeStencil[D]: Create stencil from a Lattice type (preferred constructor)

Automatically extracts local geometry and ghost width from the lattice. When mpiGrid uses auto-detect sentinels (values ≤ 0), queries GlobalArrays for the actual MPI decomposition.

Example:

let lat = newSimpleCubicLattice([16, 16, 16, 32], [2, 2, 2, 4], [1, 1, 1, 1]) let stencil = newLatticeStencil(nearestNeighborStencil[4](), lat)
proc newLatticeStencil[D: static int](pattern: StencilPattern[D]; localGeom: array[D, int]): LatticeStencil[ D]: Create stencil with zero ghost width (local-only periodic BC)
proc newLatticeStencil[D: static int](pattern: StencilPattern[D]; localGeom: array[D, int]; ghostWidth: array[D, int]): LatticeStencil[ D]: Create stencil with default SIMD grid
proc newLatticeStencil[D: static int](pattern: StencilPattern[D]; localGeom: array[D, int]; ghostWidth: array[D, int]; simdGrid: array[D, int]): LatticeStencil[D]: Create a unified lattice stencil with explicit parameters

Parameters: pattern: Stencil pattern (e.g., nearestNeighborStencil[4]()) localGeom: Local lattice dimensions (without ghosts) ghostWidth: Ghost width in each dimension simdGrid: SIMD grid for vectorization (use 1,1,1,1 for no SIMD)
proc newStencil[D: static int](name: string = ""): Stencil[D]
proc newStencilPattern[D: static int](name: string = ""): StencilPattern[D]
proc newStencilPoint[D: static int](offset: array[D, int]; id: int = 0): StencilPoint[ D]
proc newStencilPoint[D: static int](sd: SignedDirection): StencilPoint[D]
proc newStencilView[D: static int](stencil: StencilPattern[D]; localGeom: array[D, int]): StencilView[D]: Legacy constructor - creates a StencilView (use LatticeStencil instead)
proc nLanes[D: static int](s: LatticeStencil[D]): int {.inline.}
proc nOuter[D: static int](s: LatticeStencil[D]): int {.inline.}
proc nPoints[D: static int](s: StencilPattern[D]): int {.inline.}
proc nSites[D: static int](s: LatticeStencil[D]): int {.inline.}
proc offsetBufferSize[D: static int](s: LatticeStencil[D]): int {.inline.}
proc paddedToLocal[D: static int](paddedSite: int; localGeom: array[D, int]; ghostWidth: array[D, int]): int
proc paddedToLocal[D: static int](s: LatticeStencil[D]; paddedSite: int): int: Convert padded site index to local site index Returns -1 if site is in ghost region
proc pathToStencil[D: static int](path: openArray[PathStep]): StencilPattern[D]
proc plaquettePath(mu, nu: int): seq[PathStep] {....raises: [], tags: [], forbids: [].}
proc rectanglePath(mu, nu: int; a, b: int): seq[PathStep] {....raises: [], tags: [], forbids: [].}
proc shift[D: static int](s: LatticeStencil[D]; site: int; sd: SignedDirection): StencilShift[ D] {.inline.}: Get a stencil shift using SignedDirection
proc shift[D: static int](s: LatticeStencil[D]; site: int; sign: int; dir: int): StencilShift[ D] {.inline.}: Get a stencil shift for neighbor access

Usage: viewstencil.shift(n, +1, 0) # Forward neighbor in x viewstencil.shift(n, -1, 3) # Backward neighbor in t

The point index is computed from (sign, dir) using the pattern: Forward stencil points are at even indices (0, 2, 4, ...) Backward stencil points are at odd indices (1, 3, 5, ...)
proc step(dir: int; forward: bool = true): PathStep {....raises: [], tags: [], forbids: [].}

Iterators

iterator allDirections(nDim: int): SignedDirection {....raises: [], tags: [], forbids: [].}: Iterate over all directions (forward and backward)
iterator directions(nDim: int): Direction {....raises: [], tags: [], forbids: [].}
iterator neighbors[D: static int](s: LatticeStencil[D]; site: int): int {.inline.}
iterator points[D: static int](s: LatticeStencil[D]): int {.inline.}
iterator sites[D: static int](s: LatticeStencil[D]): int {.inline.}

Templates

template forEachNeighbor[D: static int](s: LatticeStencil[D]; site: int; pointVar, nbrVar: untyped; body: untyped)
template T(): Direction
template X(): Direction
template Y(): Direction
template Z(): Direction