Heat equation

A simple reference program that integrates the two-dimensional heat equation with an explicit nine-point finite-difference stencil in space and a forward-Euler step in time. The problem is decomposed over an NX × NY grid of PEs, each holding an N × N subdomain padded with a 2P wide ghost region on the bottom and on the left.

The code implements the Domain Translation algorithm presented at HPC Asia 2026 available at https://dl.acm.org/doi/10.1145/3773656.3773683 and https://arxiv.org/abs/2511.11542.

Despite its simplicity, this code performs at almost 90% of peak fp32 performance on the CS-3 or CS-4.

The animation shows a simulation run on 8 × 8 PEs with 4 × 4 grid points each and using D = 1 and dt = 1/8. Colors are on a log scale. Each frame has been shifted to compensate for the grid reassignment performed by the domain transation.

---

Algorithm

PDE and discretisation

The solver integrates the heat equation

$$\frac{du}{dt} = D \cdot \nabla^2 u$$

on a doubly periodic square grid (2D torus) of spacing h with diffusion coefficient D. The calculations are done in single precision. The computational stencil is:

         j-1      j      j+1
  i-1  [ 0.5·α     α     0.5·α ]
  i    [   α    1-6·α      α   ]     α = D·dt/(2·h²)
  i+1  [ 0.5·α     α     0.5·α ]

Total sweeps = niter · ninner. make-inputbin.py fixes niter=10, ninner=3 (30 sweeps, 5 timer samples with nsample=2) and dt=1/64, with h=1 and D=1.

Initial condition and boundary conditions

The device sets its own initial data. Every PE zeroes its x array and writes x[w][w] = params.origin. make-inputbin.py gives origin = 1 only to the (ix=0, iy=0) instance of each of the four programs, so the field starts as four unit deltas in the lowest 2×2 block of PEs, at global cells (0,0), (N,0), (0,N), (N,N). The x0.bin.* files are not loaded onto the device; they are only the oracle's copy of the same initial condition.

Boundaries are fully periodic (toroidal) in both directions, implemented in the fabric by the moats rather than in the kernel.

Domain decomposition

The problem is decomposed into a grid of NX × NY compute PEs, each owning an N × N interior grid of DOFs. Local arrays are (N+w) × (N+w+pad) with w = 2P and pad chosen in heateq.paint so the row stride N+w+pad is a multiple of 4 (to reduce memory bank conflicts). Interior cells live at local indices [w, N+w) in both axes; indices [0, w) are the ghost region. Global domain is N·NX × N·NY.

Note the ghost region is low-side only (left and below). compute() reads at most index i+P+1 ≤ N+w−1, so no high-side halo is needed or allocated. The exchange is consequently a 3-neighbour, not 4-neighbour, pattern: each PE sends right and up, and receives from left, below, and (via a relay) below-left.

Halo exchange

sendrecv() runs two parallel blocks, which also act as the ordering point between the two phases on the shared up/down sockets:

  phase 1   right[] <- x[w:N+w][N:N+w]      send right edge columns
            up[]    <- x[N:N+w][w:N+w]      send top edge rows
            x[w:N+w][0:w] <- left[]         receive into left ghost
            x[0:w][w:N+w] <- down[]         receive into bottom ghost

  phase 2   up[]    <- x[N:N+w][0:w]        forward the received left ghost
            x[0:w][0:w] <- down[]           receive into bottom-left ghost

Phase 2 is the standard two-step corner fill: the down-left diagonal neighbour's data arrives by travelling right one hop in phase 1 and up one hop in phase 2. This requires N ≥ 2P so that the rows [N, N+w) forwarded in phase 2 lie inside the rows [w, N+w) filled in phase 1 (with N=2, P=1 the two ranges coincide exactly). There is no explicit barrier; flow control is socket back-pressure.

Fabric topology, colors and routing

heateq.paint allocates six colors: c0, c1 (horizontal), c2, c3 (vertical), and the swap partners c0x, c3x. A 2×2 tile block (x00, x01, x10, x11) is replicated NX/2 × NY/2 times via hstackrep / vstackrep and placed at (1,1), so compute PEs occupy x ∈ [1,NX], y ∈ [1,NY]. The four programs differ only in which color each socket binds to, so that a sender's outbound color is the receiver's inbound color:

PE parity (global x, y)	program	right →	left ←	up →	down ←
odd x, odd y	x00	c1	c0	c2	c3
odd x, even y	x01	c1	c0	c3	c2
even x, odd y	x10	c0	c1	c2	c3
even x, even y	x11	c0	c1	c3	c2

Router settings are [C>R] on the send color and [L>C] on the receive color horizontally, [C>U] / [D>C] vertically. Because colors alternate with parity, adjacent links never collide, and only c0 (rightward, out of column NX) and c3 (upward, out of row NY) ever leave the compute block.

Those two are wrapped by four one-PE-wide moat strips — mdown at y=0, mup at y=NY+1, mleft at x=0, mright at x=NX+1 — using the X swap modifier with colorpair(c0 c0x) / colorpair(c3 c3x):

mup    : paint(c3  [DXD])
mup    : paint(c3x [DXD])
mdown  : paint(c3  [UXU])
mdown  : paint(c3x [UXU])
mleft  : paint(c0  [RXR])
mleft  : paint(c0x [RXR])
mright : paint(c0  [LXL])
mright : paint(c0x [LXL])

The return legs cross the compute block on the anti-colors — corerect : paint(c0x [R>L]), corerect : paint(c3x [U>D]). So a wavelet sent right out of column NX on c0 is reflected leftward as c0x at mright, runs back across the whole row, is reflected rightward as c0 at mleft, and enters column 1 — an x00/x01 PE whose left socket is c0. The vertical wrap is identical with c3/c3x. Fabric footprint is (NX+2) × (NY+2) PEs (48 for t1; the four corners are unpainted).

Phases of one timestep

1. exchange phase 1 (all four streams concurrent in one parallel block);
2. exchange phase 2 (corner relay);
3. compute() — nine scaled accumulations into y, then x ← y;
4. every nsample outer iterations, recordtime() appends three hardware timer words to tstamp.

---

Source files

File	Role
bench.paint	the file paint is invoked on: a copy of heateq.paint
heateq.paint	floorplan: colors, the 2×2 program block, hstackrep/vstackrep tiling, routing, and the four color-swap moats
heateq9pt.w	the PE program: data layout, stencil weights, halo exchange, stencil sweep, hardware-timer sampling, main loop
timer.w	imported by heateq9pt.w: timerinit() starts the hardware timer, recordtime() appends a sample to tstamp
Makefile	PAINTVARS, PAINTFLAGS, and the miset/miget wiring for params, x, y and tstamp
test.rc	the configuration rows (see Parameters)
make-inputbin.py	writes params.bin and the reference x0.bin.* files
check-inoutbin.py	NumPy oracle; reassembles the global field from the four per-program dumps and compares
heat_animation.gif	the animation above

---

Parameters

Set per row of test.rc; PAINTVARS = N P NX NY passes the first four to paint as set=NAME=value, where they appear as bare globals in heateq.paint.

Variable	Meaning
NX	number of compute PEs in x. Must be even (the block is 2×2 and make-inputbin.py rejects odd values).
NY	number of compute PEs in y. Must be even.
N	interior subdomain edge length per PE. Global grid is N·NX × N·NY. Must satisfy N ≥ 2P.
P	half the ghost depth: w = 2P ghost cells on the low side of each axis, and also the per-sweep translation of the solution. Only P=1 is exercised.
ARCH	sdr (WSE-3).

Runtime parameters travel in params.bin as the struct params of each PE: the sp values origin, h, dt, D, then the counts niter, ninner, nsample, also stored as sp (check-inoutbin.py raises if one is not a whole number). make-inputbin.py fixes them all: h=1, D=1, dt=1/64, niter=10, ninner=3 and nsample=2.

---

Building and running

Example run commands:

make TEST=t1 # build, simulate and check test t1 in test.rc make TEST='' # every row in test.rc make TEST=t1 GOAL=debugx # also create trace files for the portrait debugger make clean # remove all test-/ directories

### Correctness check

Commands in the make file extracts relevant application data from the core file after program execution.

Then `check-inoutbin.py` is run to check correctness of output. `check-inoutbin.py`
1. reassembles the global `(N·NY) × (N·NX)` field by de-interleaving the four
   per-program dumps and stripping both the ghost region and the row padding,
2. assembles and prints timing data
3. reruns the identical nine-point scheme in NumPy for `niter·ninner` timesteps.
4. raises if any output value is non-finite, and raises
   `"In accurate solution"` if the relative 2-norm error exceeds `1e-6`.

## Run it

The experimental SDK installs this project beside the tutorial, as `projects/heateq9pt`. Run it like a tutorial: `tungsten-tutorials -f projects/heateq9pt`.