Added: Initial world generation tool

Tools/Terrain/internal/fluvial/bucketpq.go

@@ -0,0 +1,79 @@
+package fluvial
+
+// A monotone bucket priority queue, which is what priority-flood actually needs.
+//
+// The flood pops cells in non-decreasing elevation and never pushes anything below the cell it just popped:
+// a neighbour lower than the current front is raised to it and goes to the FIFO instead. That "monotone"
+// property is exactly the condition under which a bucket queue beats a binary heap, because the read cursor
+// only ever moves forward and both operations become an append and a scan. The heap was costing about
+// log2(3.2M) = 22 comparisons and as many cache misses per operation, on two thirds of the solve's runtime.
+//
+// Elevations are quantised into fixed-width buckets. Cells inside one bucket pop in an arbitrary but
+// deterministic order (last in, first out), so a spill point can be wrong by at most one bucket width. At a
+// centimetre against a 2 km elevation range that is far below the millimetre-per-cell epsilon the flood adds
+// anyway, and it is the same approximation an integer-elevation priority-flood makes by construction.
+type bucketPQ struct {
+	lo      float64
+	width   float64
+	buckets [][]int32
+	cur     int
+	count   int
+}
+
+const bucketWidthM = 0.01
+
+func newBucketPQ(loM, hiM float64) *bucketPQ {
+	if hiM <= loM {
+		hiM = loM + 1
+	}
+	// Headroom above the top: the flood raises cells by epsilon as it fills, so the highest key pushed can
+	// sit slightly above the terrain's own maximum.
+	n := int((hiM-loM)/bucketWidthM) + 64
+	return &bucketPQ{lo: loM, width: bucketWidthM, buckets: make([][]int32, n)}
+}
+
+func (q *bucketPQ) reset() {
+	for i := range q.buckets {
+		q.buckets[i] = q.buckets[i][:0]
+	}
+	q.cur = 0
+	q.count = 0
+}
+
+func (q *bucketPQ) len() int { return q.count }
+
+func (q *bucketPQ) push(elev float32, idx int32) {
+	b := int((float64(elev) - q.lo) / q.width)
+	if b < q.cur {
+		b = q.cur // monotone: never behind the cursor, whatever rounding says
+	}
+	if b >= len(q.buckets) {
+		b = len(q.buckets) - 1
+	}
+	q.buckets[b] = append(q.buckets[b], idx)
+	q.count++
+}
+
+// pop returns the lowest cell. The cursor only moves forward, so the total scan cost over a whole flood is
+// the number of buckets, not the number of pops.
+func (q *bucketPQ) pop() int32 {
+	for q.cur < len(q.buckets) && len(q.buckets[q.cur]) == 0 {
+		q.cur++
+	}
+	if q.cur >= len(q.buckets) {
+		return -1
+	}
+	b := q.buckets[q.cur]
+	v := b[len(b)-1]
+	q.buckets[q.cur] = b[:len(b)-1]
+	q.count--
+	return v
+}
+
+// frontElev is the elevation the cursor is at, which the FIFO compares itself against.
+func (q *bucketPQ) frontElev() float32 {
+	for q.cur < len(q.buckets) && len(q.buckets[q.cur]) == 0 {
+		q.cur++
+	}
+	return float32(q.lo + float64(q.cur)*q.width)
+}

Tools/Terrain/internal/fluvial/fluvial.go

@@ -0,0 +1,512 @@
+// Package fluvial is the stream-power erosion solve: dh/dt = U - K * A^m * S^n, integrated implicitly up the
+// drainage stack by the method of Braun & Willett (2013).
+//
+// This is the reason the generator exists (D-47). Particle erosion carves the path each droplet happens to
+// take: it makes wear, but never a network. Stream power solves for drainage area first and erodes in
+// proportion to it, which is what produces a branching hierarchy, valleys whose size matches the area they
+// drain, and divides that sit where the basins either side of them put them.
+//
+// Four things happen per step, in this order:
+//
+//  1. Depressions are filled (priority-flood), because a D8 receiver graph containing a pit has no path to
+//     base level and the implicit solve has nothing to descend to. This is the only part of a step that is
+//     not O(n), so it runs every FillEvery steps, not every step.
+//  2. Receivers and the stack are computed: steepest descent to one of eight neighbours, then a depth-first
+//     ordering in which every node appears after its receiver.
+//  3. Drainage area is accumulated down the stack in reverse.
+//  4. The implicit update runs up the stack, so each node's receiver already holds its new height. This is
+//     what makes the scheme unconditionally stable in dt, and it is why a naive explicit solver is not an
+//     acceptable substitute at dt = 1500 yr.
+package fluvial
+
+import (
+	"math"
+
+	"salty/terrain/internal/field"
+	"salty/terrain/internal/thermal"
+)
+
+// Params are the stream-power constants. K is per year with A in m².
+type Params struct {
+	K         float64
+	M         float64
+	N         float64
+	DtYr      float64
+	Steps     int
+	Diffusion float64 // hillslope diffusivity, m²/yr
+	FillEvery int
+
+	// Landsliding. Stream power bounds nothing at small drainage area, so without this the hillslopes grow
+	// as steep as the uplift rate asks them to and the map becomes needles. TalusSlope is rise over run; 0
+	// disables. It runs inside the loop rather than after it, because a cap applied once at the end just
+	// shaves the tops off, while a cap applied throughout changes where the sediment goes.
+	TalusSlope    float64
+	ThermalEvery  int
+	ThermalPasses int
+
+	// CriticalAreaM2 is where channels begin. Below it the cell is a hillslope: it still rises, and
+	// diffusion and landsliding still shape it, but stream power does not incise it.
+	//
+	// This is not a tuning knob, it is a correctness fix. Stream power is a law about channels, and applying
+	// it at A = one cell says a cell that drains only itself should stand at U/(K*cellM^(2m)) — 63 degrees at
+	// the rates here. That is why the map came out as needles at a 39 degree median. It is also why relief
+	// was resolution-dependent (1020 m at 512², 2605 m at 1786² on one seed): halving the cell size halves
+	// the smallest A and steepens every divide, for ever. A critical area is a physical length, so the same
+	// landscape comes out at any resolution, which is the property the full-resolution run depends on.
+	CriticalAreaM2 float64
+
+	// ChannelTaper is the exponent on (A/Ac) below the channel head. 0 is no taper, 2 is strong.
+	ChannelTaper float64
+
+	// CriticalSlope is Sc in the nonlinear hillslope law, rise over run. Above zero it selects
+	// DiffuseNonlinear over plain linear diffusion and takes the repose clamp out of the step loop; see
+	// hillslope.go and Run. Zero keeps the old pairing of linear diffusion and an in-loop clamp.
+	CriticalSlope   float64
+	SlopeCap        float64 // where the flux stops stiffening, as a fraction of Sc
+	MaxHillslopeSub int     // the sub-step budget that bound buys
+}
+
+// Grid holds the flow topology and the scratch it is built from. Allocated once and reused across every
+// step: at 3.2 M cells the allocations would otherwise dominate the solve.
+type Grid struct {
+	W, H  int
+	CellM float64
+
+	// Base marks cells fixed at base level: the ocean. The map border is an outlet too, so what actually
+	// counts as "fixed" is Base plus the border, and that union is `fixed`. Keeping only Base in mind here is
+	// how the outlets themselves ended up being uplifted: on a map with no ocean, every border cell has
+	// Receiver == itself and Base == false, so base level rose two metres a step and the whole solve chased
+	// it. The steady-state test is what caught it.
+	Base  []bool
+	fixed []bool
+
+	Receiver []int32   // index of the cell this one drains to; itself for a base cell
+	Length   []float32 // distance to that receiver, metres
+	Stack    []int32   // every node after its receiver
+	Area     []float32 // drainage area, m²
+
+	seed      uint64 // the jitter's seed; see jitter.go and SetSeed
+	donorOff  []int32
+	donorList []int32
+	cursor    []int32
+	closed    []bool
+	pq        *bucketPQ
+	fifo      []int32
+	scratch   []float32
+}
+
+// SetElevationRange sizes the flood's bucket queue. Called once, with the manifest's elevation range plus a
+// margin, before the first step.
+func (g *Grid) SetElevationRange(loM, hiM float64) {
+	g.pq = newBucketPQ(loM, hiM)
+}
+
+var (
+	// D8, in the order (-1,-1) .. (1,1) skipping the centre. The order is fixed so that a run is reproducible,
+	// but it is no longer what decides a tie between two equally steep neighbours: a fixed order resolves
+	// every tie the same way and prints its preferred axis across any near-flat ground. See jitter.go.
+	dx8 = [8]int{-1, 0, 1, -1, 1, -1, 0, 1}
+	dy8 = [8]int{-1, -1, -1, 0, 0, 1, 1, 1}
+)
+
+func NewGrid(w, h int, cellM float64, base []bool) *Grid {
+	n := w * h
+	g := &Grid{
+		W: w, H: h, CellM: cellM, Base: base,
+		Receiver: make([]int32, n), Length: make([]float32, n), Stack: make([]int32, 0, n),
+		Area: make([]float32, n), donorOff: make([]int32, n+1), donorList: make([]int32, n),
+		closed: make([]bool, n), fifo: make([]int32, 0, n), scratch: make([]float32, n),
+	}
+	g.fixed = make([]bool, n)
+	for i := range g.fixed {
+		g.fixed[i] = g.isOutlet(i)
+	}
+	g.SetElevationRange(-2000, 4000)
+	return g
+}
+
+// FillDepressions raises closed depressions to their spill point, in place, using Barnes' improved
+// priority-flood with a plain FIFO beside the heap. The FIFO is the optimisation that matters: on real
+// terrain most cells are reached while descending into an already-flooded pit, and those never touch the
+// heap, which turns the cost from "half an hour over a run" into something affordable.
+//
+// The epsilon variant adds a millimetre of fall per cell across a flat, so filled lakes still route rather
+// than becoming a plateau the flow accumulator cannot leave. That millimetre is scattered per cell by a hash
+// of the index rather than applied uniformly: a uniform epsilon means the only gradient on a flat is the
+// flood's own traversal order, and the router then draws that order as rivers. See jitter.go.
+func (g *Grid) FillDepressions(h []float32, epsilon float32) {
+	n := g.W * g.H
+	for i := range g.closed {
+		g.closed[i] = false
+	}
+	g.pq.reset()
+	g.fifo = g.fifo[:0]
+
+	for i := 0; i < n; i++ {
+		if g.isOutlet(i) {
+			g.closed[i] = true
+			g.pq.push(h[i], int32(i))
+		}
+	}
+
+	head := 0
+	for g.pq.len() > 0 || head < len(g.fifo) {
+		var c int32
+		var celev float32
+		// Drain the FIFO while it cannot violate the queue's ordering.
+		if head < len(g.fifo) && (g.pq.len() == 0 || h[g.fifo[head]] <= g.pq.frontElev()) {
+			c = g.fifo[head]
+			head++
+			celev = h[c]
+		} else {
+			c = g.pq.pop()
+			if c < 0 {
+				break
+			}
+			celev = h[c]
+		}
+		cx := int(c) % g.W
+		cy := int(c) / g.W
+		for k := 0; k < 8; k++ {
+			nx, ny := cx+dx8[k], cy+dy8[k]
+			if nx < 0 || ny < 0 || nx >= g.W || ny >= g.H {
+				continue
+			}
+			ni := int32(ny*g.W + nx)
+			if g.closed[ni] {
+				continue
+			}
+			g.closed[ni] = true
+			if h[ni] <= celev {
+				h[ni] = celev + epsilon*(0.5+hash01(g.seed, ni))
+				g.fifo = append(g.fifo, ni)
+			} else {
+				g.pq.push(h[ni], ni)
+			}
+		}
+		// Compact the FIFO occasionally so it does not grow without bound over a whole flood.
+		if head > n/2 {
+			g.fifo = append(g.fifo[:0], g.fifo[head:]...)
+			head = 0
+		}
+	}
+}
+
+func (g *Grid) isOutlet(i int) bool {
+	if g.Base != nil && g.Base[i] {
+		return true
+	}
+	x, y := i%g.W, i/g.W
+	return x == 0 || y == 0 || x == g.W-1 || y == g.H-1
+}
+
+// ComputeReceivers picks the steepest downhill neighbour of each cell. A base cell, and any cell with no
+// lower neighbour, receives itself, which makes it a root of the stack.
+//
+// Two neighbours of equal steepness are separated by a hash of the cell and the direction rather than by the
+// fixed order of dx8/dy8. A fixed order always resolves a tie the same way, which on any near-flat surface
+// puts a systematic preference on one grid axis and shows up as rivers that run straight along it. The
+// perturbation is a tenth of a percent, so it decides near-ties and nothing else: a neighbour that is
+// genuinely steeper than another by more than that is still chosen.
+func (g *Grid) ComputeReceivers(h []float32) {
+	diag := float32(g.CellM * math.Sqrt2)
+	card := float32(g.CellM)
+	field.Rows(g.H, func(y0, y1 int) {
+		for y := y0; y < y1; y++ {
+			for x := 0; x < g.W; x++ {
+				i := int32(y*g.W + x)
+				if g.isOutlet(int(i)) {
+					g.Receiver[i] = i
+					g.Length[i] = card
+					continue
+				}
+				best := int32(-1)
+				bestJitter := float32(0)
+				bestLen := card
+				for k := 0; k < 8; k++ {
+					nx, ny := x+dx8[k], y+dy8[k]
+					if nx < 0 || ny < 0 || nx >= g.W || ny >= g.H {
+						continue
+					}
+					ni := int32(ny*g.W + nx)
+					l := card
+					if dx8[k] != 0 && dy8[k] != 0 {
+						l = diag
+					}
+					s := (h[i] - h[ni]) / l
+					if s <= 0 {
+						continue
+					}
+					// The tie-break, not a change of gradient: the comparison is jittered, the slope that
+					// is kept is not, so Length and the stream-power update see the true geometry.
+					sj := s * (1 + 1e-3*(hash01(g.seed, i*8+int32(k))-0.5))
+					if sj > bestJitter {
+						bestJitter, best, bestLen = sj, ni, l
+					}
+				}
+				if best < 0 {
+					g.Receiver[i] = i
+					g.Length[i] = card
+				} else {
+					g.Receiver[i] = best
+					g.Length[i] = bestLen
+				}
+			}
+		}
+	})
+}
+
+// BuildStack orders every node after its receiver, by counting donors into a CSR list and then walking it
+// depth-first from the roots. O(n), no recursion, and the order is fully determined by the receiver array,
+// so it does not vary between runs.
+func (g *Grid) BuildStack() {
+	n := g.W * g.H
+	for i := 0; i <= n; i++ {
+		g.donorOff[i] = 0
+	}
+	for i := 0; i < n; i++ {
+		r := g.Receiver[i]
+		if int(r) != i {
+			g.donorOff[r+1]++
+		}
+	}
+	for i := 0; i < n; i++ {
+		g.donorOff[i+1] += g.donorOff[i]
+	}
+	cursor := g.scratchInt32()
+	copy(cursor, g.donorOff[:n])
+	for i := 0; i < n; i++ {
+		r := g.Receiver[i]
+		if int(r) != i {
+			g.donorList[cursor[r]] = int32(i)
+			cursor[r]++
+		}
+	}
+
+	g.Stack = g.Stack[:0]
+	for i := 0; i < n; i++ {
+		if int(g.Receiver[i]) == i {
+			g.Stack = append(g.Stack, int32(i))
+		}
+	}
+	// Depth-first: everything already in the stack expands its donors, which land after it.
+	for read := 0; read < len(g.Stack); read++ {
+		c := g.Stack[read]
+		for d := g.donorOff[c]; d < g.donorOff[c+1]; d++ {
+			g.Stack = append(g.Stack, g.donorList[d])
+		}
+	}
+}
+
+// scratchInt32 reuses the float32 scratch as int32 storage; same width, and it saves a 12 MB allocation per
+// step at the geology grid.
+func (g *Grid) scratchInt32() []int32 {
+	if cap(g.cursor) < g.W*g.H {
+		g.cursor = make([]int32, g.W*g.H)
+	}
+	return g.cursor[:g.W*g.H]
+}
+
+// Accumulate sums drainage area down the stack in reverse, so every node has collected its whole upstream
+// catchment before its receiver is reached.
+func (g *Grid) Accumulate() {
+	cell := float32(g.CellM * g.CellM)
+	for i := range g.Area {
+		g.Area[i] = cell
+	}
+	for k := len(g.Stack) - 1; k >= 0; k-- {
+		i := g.Stack[k]
+		r := g.Receiver[i]
+		if r != i {
+			g.Area[r] += g.Area[i]
+		}
+	}
+}
+
+// StreamPower is the implicit update, walked up the stack. uplift is in metres per year and k is the local
+// erodibility; either may be nil for a uniform value.
+func (g *Grid) StreamPower(h []float32, uplift, k []float32, p Params) {
+	dt := p.DtYr
+	linear := math.Abs(p.N-1) < 1e-9
+	for _, i := range g.Stack {
+		if g.fixed[i] {
+			continue // base level is fixed: no uplift, no erosion
+		}
+		u := 0.0
+		if uplift != nil {
+			u = float64(uplift[i])
+		}
+		r := g.Receiver[i]
+		if r == i {
+			// A local minimum that the last flood has not reached yet. It still rises: skipping uplift here
+			// freezes exactly the cells that differential uplift is busy pushing up, which quietly removes
+			// the basins from the landscape between floods.
+			h[i] += float32(dt * u)
+			continue
+		}
+		kk := p.K
+		if k != nil {
+			kk *= float64(k[i])
+		}
+		hr := float64(h[r])
+		hi := float64(h[i]) + dt*u
+		area := float64(g.Area[i])
+		if p.CriticalAreaM2 > 0 && area < p.CriticalAreaM2 {
+			// Below the channel head, incision is tapered rather than switched off. Switching it off
+			// entirely is what broke: the material had nowhere to go, because the only remaining transport
+			// was landsliding, which caps slope but not height, and the map grew until 22 % of it clipped
+			// the elevation range. A taper suppresses the fine dissection that makes lowlands look like
+			// small mountains while still letting the hillslope shed its uplift into the network.
+			kk *= math.Pow(area/p.CriticalAreaM2, p.ChannelTaper)
+		}
+		a := math.Pow(area, p.M)
+		l := float64(g.Length[i])
+
+		var next float64
+		if linear {
+			f := kk * dt * a / l
+			next = (hi + f*hr) / (1 + f)
+		} else {
+			next = newtonStreamPower(hi, hr, kk*dt*a, l, p.N)
+		}
+		// A node may never fall below what it drains into; the implicit form only guarantees that while
+		// uplift has not raised the receiver past it.
+		if next < hr {
+			next = hr
+		}
+		h[i] = float32(next)
+	}
+}
+
+// newtonStreamPower solves h - hi + c*((h-hr)/l)^n = 0 for n != 1. Five iterations from the linear answer is
+// comfortably enough at the exponents anyone actually uses; it is here so the manifest's n is not a lie.
+func newtonStreamPower(hi, hr, c, l, n float64) float64 {
+	f := c / l
+	h := (hi + f*hr) / (1 + f) // the n = 1 answer, as a starting point
+	for iter := 0; iter < 5; iter++ {
+		d := h - hr
+		if d < 0 {
+			d = 0
+		}
+		s := d / l
+		fx := h - hi + c*math.Pow(s, n)/1
+		dfx := 1 + c*n*math.Pow(s, n-1)/l
+		if dfx == 0 {
+			break
+		}
+		step := fx / dfx
+		h -= step
+		if math.Abs(step) < 1e-6 {
+			break
+		}
+	}
+	if h < hr {
+		h = hr
+	}
+	return h
+}
+
+// Diffuse is hillslope diffusion, which rounds the divides and stops every channel head from being a needle.
+//
+// Explicit five-point diffusion is stable only while D*dt/dx² <= 0.25, and the defaults sit above that
+// (0.02 m²/yr at dt 1500 on 8 m cells is 0.47), so it sub-steps rather than quietly going unstable. This is
+// the sort of thing that shows up as a checkerboard three thousand steps in.
+func (g *Grid) Diffuse(h []float32, d, dt float64) {
+	if d <= 0 || dt <= 0 {
+		return
+	}
+	dx2 := g.CellM * g.CellM
+	total := d * dt / dx2
+	sub := int(math.Ceil(total / 0.2))
+	if sub < 1 {
+		sub = 1
+	}
+	alpha := float32(total / float64(sub))
+	src := h
+	tmp := g.scratch[:len(h)]
+	for s := 0; s < sub; s++ {
+		field.Rows(g.H, func(y0, y1 int) {
+			for y := y0; y < y1; y++ {
+				for x := 0; x < g.W; x++ {
+					i := y*g.W + x
+					if g.fixed[i] {
+						tmp[i] = src[i]
+						continue
+					}
+					c := src[i]
+					lap := clampAt(src, g.W, g.H, x-1, y) + clampAt(src, g.W, g.H, x+1, y) +
+						clampAt(src, g.W, g.H, x, y-1) + clampAt(src, g.W, g.H, x, y+1) - 4*c
+					tmp[i] = c + alpha*lap
+				}
+			}
+		})
+		copy(src, tmp)
+	}
+}
+
+func clampAt(a []float32, w, h, x, y int) float32 {
+	if x < 0 {
+		x = 0
+	} else if x >= w {
+		x = w - 1
+	}
+	if y < 0 {
+		y = 0
+	} else if y >= h {
+		y = h - 1
+	}
+	return a[y*w+x]
+}
+
+// Run is the whole solve. Progress is reported through log, which is what a five-minute budget needs to be
+// steerable: a run that is going wrong should say so at step 500, not at the end.
+func (g *Grid) Run(h []float32, uplift, k []float32, p Params, log func(step int, total int, elapsedPct float64)) {
+	fill := p.FillEvery
+	if fill < 1 {
+		fill = 1
+	}
+	for step := 0; step < p.Steps; step++ {
+		if step%fill == 0 {
+			g.FillDepressions(h, 1e-3)
+		}
+		g.ComputeReceivers(h)
+		g.BuildStack()
+		g.Accumulate()
+		g.StreamPower(h, uplift, k, p)
+		if p.CriticalSlope > 0 {
+			// The clamp still runs, and it still has to: a belt rising at millimetres a year asks for slopes
+			// no bounded-flux transport law can hold, which is a fact about the forcing and not about the
+			// scheme. What changes is the order and who gets the last word. The clamp cuts along the eight
+			// D8 directions and leaves grid-aligned pyramid faces; nonlinear diffusion then runs over the
+			// result with a symmetric five-point stencil and rounds them off before the next step sees them.
+			//
+			// Running the clamp only once at the end was tried and is worse: a thousand steps of unclamped
+			// growth arrive at it together, so it cuts deeply, and nothing runs afterwards to soften what it
+			// cut. The facets came back in the summits. Little and often, with diffusion last, is what keeps
+			// the constraint without printing the stencil.
+			if p.TalusSlope > 0 && p.ThermalEvery > 0 && step%p.ThermalEvery == 0 {
+				g.ClampToRepose(h, p.TalusSlope)
+			}
+			g.DiffuseNonlinear(h, p.Diffusion, p.CriticalSlope, p.SlopeCap, p.DtYr, p.MaxHillslopeSub)
+		} else {
+			g.Diffuse(h, p.Diffusion, p.DtYr)
+			if p.TalusSlope > 0 && p.ThermalEvery > 0 && step%p.ThermalEvery == 0 {
+				// The constraint first, which actually binds, then the transport, which puts scree at the
+				// foot of what the constraint cut.
+				g.ClampToRepose(h, p.TalusSlope)
+				if p.ThermalPasses > 0 {
+					thermal.Apply(h, g.W, g.H, g.CellM, p.TalusSlope, p.ThermalPasses, g.fixed, g.scratch)
+				}
+			}
+		}
+		if log != nil && p.Steps >= 10 && step%(p.Steps/10) == 0 {
+			log(step, p.Steps, float64(step)/float64(p.Steps)*100)
+		}
+	}
+	// One last fill so the result has no closed pits to hand to the detail passes, and one last routing so
+	// Area and Receiver describe the surface that is actually returned.
+	g.FillDepressions(h, 1e-3)
+	g.ComputeReceivers(h)
+	g.BuildStack()
+	g.Accumulate()
+}

Tools/Terrain/internal/fluvial/hillslope.go

@@ -0,0 +1,187 @@
+package fluvial
+
+import (
+	"math"
+
+	"salty/terrain/internal/field"
+)
+
+// Nonlinear hillslope transport: q = D*S / (1 - (S/Sc)^2), the Roering form.
+//
+// It replaces the pair of patches that stood in for a hillslope law — linear diffusion, which does not care
+// how steep the ground is, and ClampToRepose, which cares about nothing else — with one process that does
+// both jobs. As S goes to zero it *is* linear diffusion, q -> D*S, so the divides in the lowlands round over
+// exactly as before. As S approaches Sc the flux diverges, so the slope approaches Sc and never reaches it.
+//
+// The difference that shows is not in the numbers, it is in the shape. ClampToRepose cuts each cell down to
+// talus*distance along one of eight neighbour directions and pops cells in grid order within an elevation
+// bucket, so what it leaves is pyramids with faces aligned to the grid — the blocky, ruler-cut facets that
+// are visible in any preview of a mountain belt here. Nothing about that is geology; it is the D8 stencil
+// printed onto the landscape. Nonlinear diffusion approaches the same limiting angle *asymptotically* and
+// through a symmetric five-point stencil, so there is no cut, no facet and no preferred direction.
+//
+// It is also mass-conserving, which the clamp is not: the flux out of one cell is the flux into its
+// neighbour by construction, so material shed from a divide arrives at the foot of the slope rather than
+// being deleted.
+//
+// # The cost, and the honest limit
+//
+// The catch is stiffness. The tangent of the flux law, which is what sets the explicit time-step limit, is
+//
+//	D_eff(S) = D * (1 + u^2) / (1 - u^2)^2,  u = S/Sc
+//
+// and that goes to infinity at u = 1. At the defaults the linear Courant number D*dt/dx^2 is already 0.29,
+// so u = 0.9 alone asks for about seventy sub-steps a step and u = 0.95 for nearly three hundred. That is
+// not affordable over a thousand steps, so the stiffening is bounded: u is capped at SlopeCap, and if even
+// that exceeds MaxSubSteps the cap is lowered further to whatever the budget affords. The sub-step count is
+// then derived from the cap that was actually used, so the scheme stays inside its stability limit whatever
+// happens — it degrades by transporting less on the steepest ground, never by going unstable.
+//
+// Ground steeper than the cap therefore relaxes at a finite rate instead of an unbounded one, and on a
+// mountain belt rising at 2 mm/yr that is not fast enough on its own. ClampToRepose stays for exactly that,
+// as a safety pass rather than as the process that shapes the land: see Run.
+func (g *Grid) DiffuseNonlinear(h []float32, d, sc, slopeCap, dt float64, maxSub int) {
+	if d <= 0 || dt <= 0 {
+		return
+	}
+	if sc <= 0 {
+		g.Diffuse(h, d, dt) // no critical slope configured: the linear law, unchanged
+		return
+	}
+	if slopeCap <= 0 || slopeCap >= 1 {
+		slopeCap = 0.9
+	}
+	if maxSub < 1 {
+		maxSub = 1
+	}
+	dx := g.CellM
+	dx2 := dx * dx
+
+	// The steepest ground on the grid bounds D_eff for the whole call. Uplift is not applied in here and
+	// diffusion only relaxes slopes, so nothing can get steeper part-way through and invalidate the bound.
+	u := math.Min(g.maxSlopeRatio(h, sc), slopeCap)
+	f := stiffness(u)
+	// What the sub-step budget can pay for. Lowering the cap rather than truncating the sub-step count is
+	// what keeps this stable: a truncated count leaves alpha above 0.25 and the surface checkerboards a few
+	// hundred steps later, which is precisely the sort of failure that does not announce itself.
+	if budget := float64(maxSub) * 0.2 * dx2 / (d * dt); f > budget {
+		f = budget
+		u = invStiffness(f)
+	}
+	if f < 1 {
+		// The budget cannot buy even the linear law. It is not optional: D*dt/dx^2 alone may need several
+		// sub-steps and going without them is an unstable scheme, so MaxSubSteps bounds the nonlinear
+		// *enhancement* and never the stability floor underneath it.
+		f = 1
+		u = 0
+	}
+	sub := int(math.Ceil(d * f * dt / dx2 / 0.2))
+	if sub < 1 {
+		sub = 1
+	}
+	dtSub := dt / float64(sub)
+	coeff := float32(d * dtSub / dx2)
+	uCap := float32(u)
+	// The height difference across one cell that *is* Sc. flux works in height differences rather than
+	// slopes, so the cell spacing has to be folded into the critical value here; leaving it out makes u a
+	// factor of dx too large, which pins every face against the cap and quietly turns the whole law into
+	// linear diffusion with a constant multiplier.
+	dhCrit := float32(sc * dx)
+
+	src := h
+	tmp := g.scratch[:len(h)]
+	for s := 0; s < sub; s++ {
+		field.Rows(g.H, func(y0, y1 int) {
+			for y := y0; y < y1; y++ {
+				for x := 0; x < g.W; x++ {
+					i := y*g.W + x
+					if g.fixed[i] {
+						tmp[i] = src[i] // base level: held, and whatever arrives here has left the system
+						continue
+					}
+					c := src[i]
+					// The net inflow over the four faces. Each face is evaluated from both of its cells,
+					// which costs twice and buys a gather: no two goroutines ever write the same cell.
+					net := flux(clampAt(src, g.W, g.H, x-1, y)-c, dhCrit, uCap) +
+						flux(clampAt(src, g.W, g.H, x+1, y)-c, dhCrit, uCap) +
+						flux(clampAt(src, g.W, g.H, x, y-1)-c, dhCrit, uCap) +
+						flux(clampAt(src, g.W, g.H, x, y+1)-c, dhCrit, uCap)
+					tmp[i] = c + coeff*net
+				}
+			}
+		})
+		copy(src, tmp)
+	}
+}
+
+// flux is q/D for one face, in height differences rather than slopes: one factor of the cell spacing cancels
+// against the divergence and is carried in coeff instead. dhCrit is the height difference that corresponds to
+// Sc across one cell, so dh/dhCrit is exactly S/Sc. u is capped so the denominator cannot reach zero.
+func flux(dh, dhCrit, uCap float32) float32 {
+	u := dh / dhCrit
+	if u < 0 {
+		u = -u
+	}
+	if u > uCap {
+		u = uCap
+	}
+	return dh / (1 - u*u)
+}
+
+// stiffness is D_eff/D at a given u = S/Sc: the factor by which the nonlinear law shortens the stable step.
+func stiffness(u float64) float64 {
+	q := 1 - u*u
+	return (1 + u*u) / (q * q)
+}
+
+// invStiffness inverts it. Bisection because stiffness is monotone on [0,1) and this runs once per call, so
+// there is nothing to gain from being cleverer and something to lose from being wrong.
+func invStiffness(f float64) float64 {
+	if f <= 1 {
+		return 0
+	}
+	lo, hi := 0.0, 0.999999
+	for i := 0; i < 60; i++ {
+		mid := (lo + hi) / 2
+		if stiffness(mid) < f {
+			lo = mid
+		} else {
+			hi = mid
+		}
+	}
+	return lo
+}
+
+// maxSlopeRatio is the steepest face on the grid as a fraction of Sc. Cardinal neighbours only, because those
+// are the faces the five-point stencil actually transports across.
+func maxSlopeRatio(h []float32, w, hgt int, sc, cellM float64) float64 {
+	var maxDiff float32
+	for y := 0; y < hgt; y++ {
+		for x := 0; x < w; x++ {
+			i := y*w + x
+			c := h[i]
+			if x+1 < w {
+				if dv := abs32(h[i+1] - c); dv > maxDiff {
+					maxDiff = dv
+				}
+			}
+			if y+1 < hgt {
+				if dv := abs32(h[i+w] - c); dv > maxDiff {
+					maxDiff = dv
+				}
+			}
+		}
+	}
+	return float64(maxDiff) / cellM / sc
+}
+
+func (g *Grid) maxSlopeRatio(h []float32, sc float64) float64 {
+	return maxSlopeRatio(h, g.W, g.H, sc, g.CellM)
+}
+
+func abs32(v float32) float32 {
+	if v < 0 {
+		return -v
+	}
+	return v
+}

Tools/Terrain/internal/fluvial/hillslope_test.go

@@ -0,0 +1,239 @@
+package fluvial
+
+import (
+	"math"
+	"testing"
+)
+
+// The three properties the nonlinear law is being trusted for. Each one is a thing the pair it replaces got
+// wrong, so each is worth a test rather than an assurance.
+
+// TestDiffuseNonlinearConservesMass is the property ClampToRepose does not have: material shed from a divide
+// has to arrive somewhere, not be deleted.
+//
+// The border is always an outlet, so the sum over the whole grid cannot be conserved by construction — base
+// level is a sink and is meant to be. The check is therefore over an interior that the disturbance never
+// reaches: a bump in the middle of a grid big enough that nothing has diffused to the edge by the time the
+// run ends.
+func TestDiffuseNonlinearConservesMass(t *testing.T) {
+	const (
+		w, h  = 101, 101
+		cellM = 10.0
+		sc    = 0.7
+	)
+	base := make([]bool, w*h)
+	g := NewGrid(w, h, cellM, base)
+	g.SetElevationRange(-100, 2000)
+
+	field := make([]float32, w*h)
+	for y := 0; y < h; y++ {
+		for x := 0; x < w; x++ {
+			d := math.Hypot(float64(x-w/2), float64(y-h/2)) * cellM
+			field[y*w+x] = float32(math.Max(0, 300-1.5*d)) // a cone well past Sc
+		}
+	}
+
+	sum := func() float64 {
+		var s float64
+		for y := 2; y < h-2; y++ {
+			for x := 2; x < w-2; x++ {
+				s += float64(field[y*w+x])
+			}
+		}
+		return s
+	}
+	before := sum()
+	for i := 0; i < 40; i++ {
+		g.DiffuseNonlinear(field, 0.02, sc, 0.9, 1500, 24)
+	}
+	after := sum()
+
+	// The cone is 300 m tall and the interior holds millions of cubic metres; a tenth of a percent is a very
+	// tight bound on forty steps of an explicit scheme in float32.
+	rel := math.Abs(after-before) / before
+	t.Logf("interior mass %.1f -> %.1f, relative change %.2e", before, after, rel)
+	if rel > 1e-3 {
+		t.Errorf("interior mass changed by %.3f%%; the flux is not antisymmetric", rel*100)
+	}
+}
+
+// TestDiffuseNonlinearLimitsSlope is the self-limiting property, and the only honest way to test it is under
+// uplift. Without uplift every diffusion law flattens everything eventually, nonlinear included, so a
+// relaxing cone proves nothing. What distinguishes the two laws is where they come to rest against a forcing:
+// linear diffusion has no limiting angle at all and lets relief grow to U*L^2/(2D), which at these numbers is
+// a kilometre and slopes many times Sc, while the nonlinear law's flux diverges as the slope approaches Sc so
+// the landscape settles near it whatever U is. That is the entire reason for the change, so it is the test.
+func TestDiffuseNonlinearLimitsSlope(t *testing.T) {
+	// The forcing is chosen so the question is about the law and not about the sub-step budget. A hillslope of
+	// half-width L under uplift U comes to rest, under the linear law, at a maximum slope of U*L/D; here that is
+	// 0.8, twice Sc, so linear diffusion visibly fails to limit. The nonlinear law can hold Sc only while its
+	// flux at the cap, D*Sc/(1-uCap^2), still exceeds U*L, and at these numbers it does with room to spare — so
+	// a failure here is the law's, not the budget's. Push U much higher and no bounded-flux law holds Sc; that
+	// is the regime ClampToRepose exists for, and Run keeps it for exactly that reason.
+	const (
+		w, h    = 41, 41
+		cellM   = 10.0
+		sc      = 0.4
+		upliftM = 2e-4 // 0.2 mm/yr
+		dt      = 1000.0
+		steps   = 5000
+	)
+	grow := func(nonlinear bool) float64 {
+		base := make([]bool, w*h)
+		g := NewGrid(w, h, cellM, base)
+		g.SetElevationRange(-100, 8000)
+		f := make([]float32, w*h)
+		for i := 0; i < steps; i++ {
+			for j := range f {
+				if !g.fixed[j] {
+					f[j] += float32(upliftM * dt)
+				}
+			}
+			if nonlinear {
+				g.DiffuseNonlinear(f, 0.05, sc, 0.9, dt, 24)
+			} else {
+				g.Diffuse(f, 0.05, dt)
+			}
+		}
+		return maxCardinalSlope(f, w, h, cellM)
+	}
+	lin := grow(false)
+	non := grow(true)
+	t.Logf("after %.1f Myr at %.1f mm/yr: linear reaches slope %.3f, nonlinear %.3f (Sc %.3f)",
+		steps*dt/1e6, upliftM*1000, lin, non, sc)
+
+	if non > sc {
+		t.Errorf("nonlinear settled at %.3f, above Sc %.3f: the flux is not stiffening", non, sc)
+	}
+	if non < sc*0.4 {
+		t.Errorf("nonlinear settled at %.3f, far below Sc %.3f: it is over-transporting", non, sc)
+	}
+	// The discriminating statement: under one forcing, the linear law overshoots the critical slope and the
+	// nonlinear law does not. If linear stays under it too, the forcing was too gentle to test anything.
+	if lin <= sc {
+		t.Errorf("linear only reached %.3f against Sc %.3f; the forcing is too weak to tell the laws apart", lin, sc)
+	}
+}
+
+// TestDiffuseNonlinearIsStable catches the failure that does not announce itself. An explicit scheme run past
+// its stability limit does not blow up on the first step; it grows a checkerboard over hundreds of them, and
+// by then the run is finished and the artefact looks like texture. A checkerboard is the mode a five-point
+// stencil goes unstable in, so it is what the test starts from: a stable scheme damps it towards flat.
+//
+// The settings put the sub-step logic where it has to choose. The initial field is far past Sc, so the cap
+// engages; the budget is well below what u = 0.95 would want, so the cap has to be lowered rather than the
+// sub-step count truncated. Truncating is the tempting, wrong branch and is what this is here to catch.
+//
+// Only the interior is measured. The border is an outlet and is held fixed by design, so it keeps its initial
+// values for ever and reading it back tells you nothing about the scheme.
+func TestDiffuseNonlinearIsStable(t *testing.T) {
+	const (
+		w, h  = 64, 64
+		cellM = 8.0
+		sc    = 0.7
+	)
+	base := make([]bool, w*h)
+	g := NewGrid(w, h, cellM, base)
+	g.SetElevationRange(-1000, 4000)
+
+	// The checkerboard goes in the interior only. The border is an outlet and is held fixed, so a
+	// checkerboard written across it is a permanent forcing that keeps re-injecting the mode into the first
+	// interior ring — the scheme would then be blamed for a boundary condition.
+	field := make([]float32, w*h)
+	for y := 1; y < h-1; y++ {
+		for x := 1; x < w-1; x++ {
+			if (x+y)%2 == 0 {
+				field[y*w+x] = 200
+			}
+		}
+	}
+	for i := range field {
+		if g.fixed[i] {
+			field[i] = 100 // flat base level, the mean of the checkerboard
+		}
+	}
+	for i := 0; i < 500; i++ {
+		g.DiffuseNonlinear(field, 0.02, sc, 0.95, 1500, 24)
+	}
+	lo, hi := float32(math.Inf(1)), float32(math.Inf(-1))
+	for y := 1; y < h-1; y++ {
+		for x := 1; x < w-1; x++ {
+			v := field[y*w+x]
+			if math.IsNaN(float64(v)) || math.IsInf(float64(v), 0) {
+				t.Fatalf("hillslope diffusion produced %v at %d,%d", v, x, y)
+			}
+			if v < lo {
+				lo = v
+			}
+			if v > hi {
+				hi = v
+			}
+		}
+	}
+	t.Logf("after 500 steps the interior spans %.3f..%.3f m, from a 200 m checkerboard", lo, hi)
+	if hi-lo > 1 {
+		t.Errorf("the checkerboard is still %.1f m after 500 steps: it is not being damped", hi-lo)
+	}
+}
+
+// TestDiffuseNonlinearMatchesLinearWhenGentle pins the other end of the law. Well below Sc the two must agree
+// closely, because that is the claim that lets this replace linear diffusion outright rather than sit beside
+// it: the lowlands must not change when the switch is thrown.
+func TestDiffuseNonlinearMatchesLinearWhenGentle(t *testing.T) {
+	const (
+		w, h  = 64, 64
+		cellM = 10.0
+		sc    = 1.0 // Sc far above anything in the field, so u stays near zero
+	)
+	base := make([]bool, w*h)
+	a := make([]float32, w*h)
+	for y := 0; y < h; y++ {
+		for x := 0; x < w; x++ {
+			// A gentle bump: peak slope about 0.01, one percent of Sc.
+			d := math.Hypot(float64(x-w/2), float64(y-h/2)) * cellM
+			a[y*w+x] = float32(3 * math.Exp(-d*d/(2*100*100)))
+		}
+	}
+	b := make([]float32, w*h)
+	copy(b, a)
+
+	ga := NewGrid(w, h, cellM, base)
+	ga.SetElevationRange(-100, 100)
+	gb := NewGrid(w, h, cellM, base)
+	gb.SetElevationRange(-100, 100)
+	for i := 0; i < 20; i++ {
+		ga.DiffuseNonlinear(a, 0.02, sc, 0.9, 1500, 24)
+		gb.Diffuse(b, 0.02, 1500)
+	}
+
+	var worst float64
+	for i := range a {
+		if d := math.Abs(float64(a[i] - b[i])); d > worst {
+			worst = d
+		}
+	}
+	t.Logf("worst divergence from linear diffusion over 20 steps: %.4f m", worst)
+	if worst > 0.01 {
+		t.Errorf("nonlinear and linear diffusion differ by %.4f m at u ~ 0.01; they should agree", worst)
+	}
+}
+
+func maxCardinalSlope(h []float32, w, hgt int, cellM float64) float64 {
+	var worst float64
+	for y := 0; y < hgt; y++ {
+		for x := 0; x < w; x++ {
+			i := y*w + x
+			if x+1 < w {
+				if s := math.Abs(float64(h[i+1]-h[i])) / cellM; s > worst {
+					worst = s
+				}
+			}
+			if y+1 < hgt {
+				if s := math.Abs(float64(h[i+w]-h[i])) / cellM; s > worst {
+					worst = s
+				}
+			}
+		}
+	}
+	return worst
+}

Tools/Terrain/internal/fluvial/jitter.go

@@ -0,0 +1,36 @@
+package fluvial
+
+// Deterministic per-cell jitter, and why a router needs one.
+//
+// D8 lets a cell drain to one of eight neighbours, so every channel is a chain of 0, 45 and 90 degree
+// segments. In mountains the slope hides it. On a plain it is the dominant artefact, and for a specific
+// reason: across a filled flat the only gradient present is the priority-flood's own epsilon, one millimetre
+// a cell, applied in the order the flood happened to reach the cells. The router then faithfully follows the
+// flood's traversal geometry and draws it as rivers — ruler-straight diagonals, the polygonal network that
+// killed the first attempt at flat plains.
+//
+// The fix is to stop the epsilon being uniform. A hash of the cell index scatters it by plus or minus half,
+// which is far below anything that matters to the solve (a millimetre against metre-scale relief) and far
+// above the difference the flood's ordering would otherwise leave, so the descent direction on a flat is
+// decided by the hash rather than by scan order. The same hash breaks near-ties between two equally steep
+// neighbours, which is the other place a fixed direction order leaks a grid axis into the result.
+//
+// It is a hash rather than a random source because cross-cutting rule 12 is determinism from a seed: the
+// value for a cell must not depend on how many cells were visited before it, on which goroutine ran, or on
+// how many steps the solve has taken.
+
+// hash01 is splitmix64 finalised to the unit interval. Cheap, no state, and well enough distributed that
+// neighbouring indices get unrelated values — which is the whole requirement here.
+func hash01(seed uint64, i int32) float32 {
+	x := seed ^ (uint64(uint32(i)) * 0x9e3779b97f4a7c15)
+	x ^= x >> 30
+	x *= 0xbf58476d1ce4e5b9
+	x ^= x >> 27
+	x *= 0x94d049bb133111eb
+	x ^= x >> 31
+	return float32(x>>11) / float32(1<<53)
+}
+
+// SetSeed ties the jitter to the run's seed, so two seeds do not share the same flat-routing geometry.
+// Zero is a perfectly good seed; it is the default and nothing depends on it being set.
+func (g *Grid) SetSeed(seed int64) { g.seed = uint64(seed)*0x9e3779b97f4a7c15 + 0x243f6a8885a308d3 }

Tools/Terrain/internal/fluvial/repose.go

@@ -0,0 +1,74 @@
+package fluvial
+
+import "math"
+
+// ClampToRepose enforces a maximum slope everywhere: no cell may stand above a neighbour by more than
+// talus * distance. It returns the mean thickness removed, in metres.
+//
+// This replaces iterating thermal.Apply inside the solve, which could not do the job however many passes it
+// was given (measured: 3, 10 and 40 passes all left the steepest land slope at 64 degrees against a 22 degree
+// repose). The reason is structural rather than a bug. That routine moves half the excess downhill, so on a
+// *uniform* over-steep slope every cell sheds exactly as much as it receives, the net change is zero, and the
+// slope is a fixed point. It relaxes only where the downhill flux diverges — which is why it cuts a cone,
+// whose contours converge, and why it cannot touch a planar hillside.
+//
+// So the constraint is imposed directly instead. This is the priority-flood mirrored: pop cells in ascending
+// elevation, and lower any neighbour standing higher than the repose angle allows. Because a lowered cell is
+// set to h[c] + talus*d, which is at or above the elevation just popped, the queue stays monotone and the
+// bucket queue works unchanged. One pass, O(n) with the bucket queue, and the constraint holds globally when
+// it returns.
+//
+// It is not mass-conserving: the material is removed rather than piled at the foot of the slope. That is the
+// deliberate simplification, because in this landscape the foot of a hillslope is a channel and the channel
+// exports the sediment anyway. The mean thickness removed is returned so a run can report it, and a run that
+// removes a suspicious amount is saying its uplift and its repose angle disagree.
+func (g *Grid) ClampToRepose(h []float32, talus float64) float64 {
+	if talus <= 0 {
+		return 0
+	}
+	n := g.W * g.H
+	for i := range g.closed {
+		g.closed[i] = false
+	}
+	g.pq.reset()
+	for i := 0; i < n; i++ {
+		g.pq.push(h[i], int32(i))
+	}
+
+	card := talus * g.CellM
+	diag := talus * g.CellM * math.Sqrt2
+	var removed float64
+
+	for g.pq.len() > 0 {
+		c := g.pq.pop()
+		if c < 0 {
+			break
+		}
+		if g.closed[c] {
+			continue
+		}
+		g.closed[c] = true
+		cx, cy := int(c)%g.W, int(c)/g.W
+		for k := 0; k < 8; k++ {
+			nx, ny := cx+dx8[k], cy+dy8[k]
+			if nx < 0 || ny < 0 || nx >= g.W || ny >= g.H {
+				continue
+			}
+			ni := int32(ny*g.W + nx)
+			if g.closed[ni] || g.fixed[ni] {
+				continue
+			}
+			allow := card
+			if dx8[k] != 0 && dy8[k] != 0 {
+				allow = diag
+			}
+			limit := h[c] + float32(allow)
+			if h[ni] > limit {
+				removed += float64(h[ni] - limit)
+				h[ni] = limit
+				g.pq.push(limit, ni)
+			}
+		}
+	}
+	return removed / float64(n)
+}

Tools/Terrain/internal/fluvial/repose_test.go

@@ -0,0 +1,50 @@
+package fluvial
+
+import (
+	"math"
+	"testing"
+)
+
+// TestClampToReposeCutsACone is the smallest possible check on the constraint: a cone far steeper than the
+// repose angle must come back at or under it.
+func TestClampToReposeCutsACone(t *testing.T) {
+	const (
+		w, h  = 81, 81
+		cellM = 10.0
+		talus = 0.4 // about 22 degrees
+	)
+	base := make([]bool, w*h)
+	field := make([]float32, w*h)
+	for y := 0; y < h; y++ {
+		for x := 0; x < w; x++ {
+			d := math.Hypot(float64(x-w/2), float64(y-h/2)) * cellM
+			field[y*w+x] = float32(math.Max(0, 800-2.0*d)) // slope 2.0, five times repose
+		}
+	}
+	g := NewGrid(w, h, cellM, base)
+	g.SetElevationRange(-100, 2000)
+	removed := g.ClampToRepose(field, talus)
+
+	worst := 0.0
+	for y := 0; y < h; y++ {
+		for x := 0; x < w; x++ {
+			for k := 0; k < 8; k++ {
+				nx, ny := x+dx8[k], y+dy8[k]
+				if nx < 0 || ny < 0 || nx >= w || ny >= h {
+					continue
+				}
+				d := cellM
+				if dx8[k] != 0 && dy8[k] != 0 {
+					d = cellM * math.Sqrt2
+				}
+				if s := float64(field[y*w+x]-field[ny*w+nx]) / d; s > worst {
+					worst = s
+				}
+			}
+		}
+	}
+	t.Logf("removed %.2f m mean; steepest slope now %.3f (repose %.3f)", removed, worst, talus)
+	if worst > talus*1.02 {
+		t.Errorf("steepest slope %.3f exceeds repose %.3f", worst, talus)
+	}
+}