87 lines
3.7 KiB
Go
87 lines
3.7 KiB
Go
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
|
|
}
|
|
// Pushed with a jittered bucket, not a plain one. The constraint this pass imposes is isotropic; the
|
|
// order it imposed it in was not. Every cell went in in flat-index order and the queue pops last-in
|
|
// first-out within a bucket, so on ground flat to within a centimetre - which is most of a hillside -
|
|
// cells popped bottom-right to top-left, and whichever popped first decided which of its neighbours got
|
|
// cut. That is where the grid-aligned pyramid faces came from, and it is one hash away from not being
|
|
// there. See bucketpq.go.
|
|
g.pq.reset()
|
|
for i := 0; i < n; i++ {
|
|
x, y := i%g.W, i/g.W
|
|
g.pq.pushJittered(h[i], int32(i), (hashXY(g.seed, g.worldX(x), g.worldY(y), jitterReposeOrder)-0.5)*2*reposeOrderBuckets)
|
|
}
|
|
|
|
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
|
|
}
|
|
// The same tie-break ComputeReceivers uses and for the same reason: a fixed allowance resolves
|
|
// every near-tie the same way and prints its preferred axis. A tenth of a percent, keyed on the
|
|
// cell being cut, so what a cell is allowed in a direction does not depend on which neighbour
|
|
// reached it first.
|
|
allow *= 1 + 1e-3*(float64(hashXY(g.seed, g.worldX(nx), g.worldY(ny), int32(k)+jitterReposeAllow))-0.5)
|
|
limit := h[c] + float32(allow)
|
|
if h[ni] > limit {
|
|
removed += float64(h[ni] - limit)
|
|
h[ni] = limit
|
|
g.pq.pushJittered(limit, ni, (hashXY(g.seed, g.worldX(nx), g.worldY(ny), jitterReposeOrder)-0.5)*2*reposeOrderBuckets)
|
|
}
|
|
}
|
|
}
|
|
return removed / float64(n)
|
|
}
|