FEEC data structures: StencilVector and StencilMatrix#

This tutorial is the third FEEC tutorial in the developer guide. It focuses on the core linear-algebra data structures used for FEEC coefficients and operators in Struphy.

If you are new to Struphy FEEC, start with:

  1. FEEC basics

  2. FEEC boundary conditions

Here we drill down into:

  1. StencilVector: storage of scalar FE coefficients

  2. StencilMatrix: storage of stencil-like sparse operators

Particle data structures are intentionally not covered in this notebook.

1. Setup with the current Derham API#

We use the current API based on TensorProductGrid and DerhamOptions, as introduced in the FEEC basics tutorial.

[1]:

import numpy as np
from feectools.linalg.stencil import StencilMatrix, StencilVector

from struphy.feec.psydac_derham import Derham
from struphy.io.options import DerhamOptions
from struphy.topology.grids import TensorProductGrid

grid = TensorProductGrid(num_elements=(8, 6, 4))
derham = Derham(grid, DerhamOptions())

print('Discrete spaces available in coeff_spaces:')
for key in ['H1', 'Hcurl', 'Hdiv', 'L2']:
    print(f'  {key:5s} -> {type(derham.coeff_spaces[key]).__name__}')

Discrete spaces available in coeff_spaces:
  H1    -> StencilVectorSpace
  Hcurl -> BlockVectorSpace
  Hdiv  -> BlockVectorSpace
  L2    -> StencilVectorSpace

2. StencilVector in practice#

A StencilVector stores scalar FE coefficients in a distributed layout with ghost regions.

In FEEC, 0-form (H1) and 3-form (L2) coefficients are StencilVector objects.

[2]:

# Option A: build a zero vector from a coefficient space
x = derham.coeff_spaces['H1'].zeros()
print('x type from zeros():', type(x))

# Option B: project a function with the commuting projector P0
f = lambda e1, e2, e3: np.cos(2 * np.pi * e1) + 0.25 * np.sin(2 * np.pi * e2)
xh = derham.P0(f)
print('xh type from P0:', type(xh))

assert isinstance(xh, StencilVector)

x type from zeros(): <class 'feectools.linalg.stencil.StencilVector'>
xh type from P0: <class 'feectools.linalg.stencil.StencilVector'>

[3]:

print('Starts (global index owned by this rank):', xh.starts)
print('Ends   (global index owned by this rank):', xh.ends)
print('Pads   (ghost widths per direction):', xh.pads)
print('shape  (logical tensor shape):', xh.shape)
print('_data shape (includes ghost layers):', xh._data.shape)

Starts (global index owned by this rank): (np.int64(0), np.int64(0), np.int64(0))
Ends   (global index owned by this rank): (np.int64(7), np.int64(5), np.int64(3))
Pads   (ghost widths per direction): (1, 1, 1)
shape  (logical tensor shape): (np.int64(192),)
_data shape (includes ghost layers): (10, 8, 6)

2.1 Global indexing vs local storage#

Using xh[i, j, k] addresses coefficients by global FE indices.

Using xh._data[...] addresses the local storage array directly, where ghost layers are explicit.

[4]:

# Fill owned coefficients with a recognizable pattern through global indexing
s = xh.starts
e = xh.ends

for i in range(s[0], e[0] + 1):
    for j in range(s[1], e[1] + 1):
        for k in range(s[2], e[2] + 1):
            xh[i, j, k] = i + 10 * j + 100 * k

print('Sample values with global indexing:')
print('xh[s0, s1, s2] =', xh[s[0], s[1], s[2]])
print('xh[e0, e1, e2] =', xh[e[0], e[1], e[2]])

Sample values with global indexing:
xh[s0, s1, s2] = 0.0
xh[e0, e1, e2] = 357.0

[5]:

# Compare compact coefficient view with local storage
flat = xh.toarray()
print('xh.toarray() shape:', flat.shape)
print('First 8 compact coefficients:', flat[:8])

print('Direct local array shape (with ghosts):', xh._data.shape)

xh.toarray() shape: (192,)
First 8 compact coefficients: [  0. 100. 200. 300.  10. 110. 210. 310.]
Direct local array shape (with ghosts): (10, 8, 6)

2.2 Relation to vector-valued FE spaces#

For Hcurl and Hdiv, coefficient vectors are block-structured (one block per component). Each block is itself stencil-based.

For projector calls like P1, vector-valued arguments should be passed as a 3-tuple of scalar callables (f1, f2, f3), one callable per component.

[6]:

g1 = lambda e1, e2, e3: e1 + e2
g2 = lambda e1, e2, e3: e2 + e3
g3 = lambda e1, e2, e3: e3 + e1
u_hcurl = derham.P1((g1, g2, g3))

print('Type of Hcurl coefficients:', type(u_hcurl))
print('Number of blocks/components:', len(u_hcurl.blocks))
print('Type of first block:', type(u_hcurl.blocks[0]))

Type of Hcurl coefficients: <class 'feectools.linalg.block.BlockVector'>
Number of blocks/components: 3
Type of first block: <class 'feectools.linalg.stencil.StencilVector'>

3. StencilMatrix in practice#

A StencilMatrix represents a linear operator between stencil vector spaces. It stores row ownership plus a compact set of local diagonals (stencil offsets), not a dense matrix.

[7]:

V0 = derham.coeff_spaces['H1']
A = StencilMatrix(V0, V0)

print('A type:', type(A))
print('A codomain starts:', A.codomain.starts)
print('A codomain ends  :', A.codomain.ends)
print('A pads:', A.pads)
print('A._data shape:', A._data.shape)

A type: <class 'feectools.linalg.stencil.StencilMatrix'>
A codomain starts: (np.int64(0), np.int64(0), np.int64(0))
A codomain ends  : (np.int64(7), np.int64(5), np.int64(3))
A pads: (1, 1, 1)
A._data shape: (10, 8, 6, 3, 3, 3)

3.1 Understanding matrix indexing#

In 3D, indexing is A[i, j, k, s1, s2, s3]:

  • (i, j, k) are global row indices owned by the current rank

  • (s1, s2, s3) are stencil offsets (diagonal at (0, 0, 0))

[8]:

s = A.codomain.starts
e = A.codomain.ends

# Build a simple 1D second-derivative-like stencil in eta1 direction
for i in range(s[0], e[0] + 1):
    for j in range(s[1], e[1] + 1):
        for k in range(s[2], e[2] + 1):
            A[i, j, k, 0, 0, 0] = 2.0
            A[i, j, k, -1, 0, 0] = -1.0
            A[i, j, k, 1, 0, 0] = -1.0

print('Stencil entries at one sample row:')
print('diag      =', A[s[0], s[1], s[2], 0, 0, 0])
print('left diag =', A[s[0], s[1], s[2], -1, 0, 0])
print('rightdiag =', A[s[0], s[1], s[2], 1, 0, 0])

Stencil entries at one sample row:
diag      = 2.0
left diag = -1.0
rightdiag = -1.0

[9]:

# Apply matrix to a vector and inspect the result
y = A.dot(xh)
print('Type of y = A.dot(xh):', type(y))
print('y.toarray() shape:', y.toarray().shape)
print('First 8 entries of y.toarray():', y.toarray()[:8])

Type of y = A.dot(xh): <class 'feectools.linalg.stencil.StencilVector'>
y.toarray() shape: (192,)
First 8 entries of y.toarray(): [ -1.70710678  98.29289322 198.29289322 298.29289322   8.07638687
 108.07638687 208.07638687 308.07638687]

4. Practical developer notes#

  1. toarray() is ideal for debugging small problems, but avoid it in production-scale runs because it materializes owned coefficients as a dense 1D array.

  2. Prefer operations at the vector/matrix level (A.dot(x), projector calls, operator composition) over manual loops in performance-critical code.

  3. In parallel runs, ghost regions and ownership become essential for correctness. This notebook keeps the focus on conceptual usage, while MPI details can be explored in specialized FEEC tests/examples.

For context on FEEC spaces and projectors, revisit the basics tutorial. For boundary-condition effects on projected fields, see the boundary-conditions tutorial.