nx_water_erosion.nx source
↩ module page · 781 lines · 32053 B
1// nx_water_erosion.nx -- full Mei 2007 hydraulic erosion simulation.
2//
3// Mei, Decaudin, Hu 2007: "Fast Hydraulic Erosion Simulation and
4// Visualization on GPU". Full implementation per cardinal
5// `feedback-maximum-capability-no-simplification`: this primitive
6// ships the entire shallow-water + sediment-transport algorithm,
7// not a slope-magnitude approximation.
8//
9// Algorithm (Mei 2007 sections 3.1 - 3.5):
10//
11// STAGE 1 -- WATER INCREMENT (rain)
12// d(t+dt) += dt * r_t where r_t = rain rate per cell
13//
14// STAGE 2 -- FLOW SIMULATION (pipe model)
15// For each cell + each cardinal direction X:
16// f_X(t+dt) = max(0, f_X(t) + dt * A * g * dh_X / L)
17// where dh_X = h(this) + d(this) - h(neighbour_X) - d(neighbour_X),
18// A = pipe cross-section, L = pipe length, g = gravity.
19//
20// Outflow scaling: if total_outflow * dt > available_water * cell_area,
21// scale all 4 outflows by min(1, available_water * cell_area /
22// (total_outflow * dt)).
23//
24// STAGE 3 -- WATER UPDATE (mass conservation)
25// d(t+dt) = d(t) + dt * (total_inflow - total_outflow) / cell_area
26//
27// STAGE 4 -- VELOCITY FIELD
28// delta_W_x = (f_R(left_neighbour) - f_L(this) + f_R(this) - f_L(right_neighbour)) / 2
29// u = delta_W_x / (cell_size * d_avg) (similar for v, y axis)
30//
31// STAGE 5 -- EROSION + DEPOSITION
32// C = K_c * |sin(local_tilt)| * |velocity| (sediment capacity)
33// if C > s: erode -> h -= K_s * (C - s); s += K_s * (C - s)
34// if C < s: deposit -> h += K_d * (s - C); s -= K_d * (s - C)
35//
36// STAGE 6 -- SEDIMENT TRANSPORT (semi-Lagrangian)
37// s_new(x, y) = s_old(x - u*dt, y - v*dt) -- bilinear interp
38//
39// STAGE 7 -- EVAPORATION
40// d(t+dt) = d(t) * (1 - K_e * dt)
41//
42// Per-cell state stored in a flat i64 buffer with 9 fields:
43// h, d, s, f_L, f_R, f_T, f_B, u, v
44//
45// Caller allocates the buffer (cell_count * 9 * 8 bytes), initialises
46// h from a heightmap source, and ticks the simulation N times.
47//
48// Q14 fixed-point throughout. Performance: tight per-cell loops with
49// neighbour reads via cell-index arithmetic; no virtual dispatch; no
50// per-tick allocation. 64x64 grid = 9 * 4096 = 36864 i64 = ~290 KB
51// per buffer.
52//
53// Convenience: the v1 stateless slope-magnitude API
54// (nx_water_erosion_delta + nx_water_slope_at + nx_water_erosion_at)
55// is preserved -- callers wanting a stateless quick-check don't have
56// to drive the full simulation. The 4 mode-routed strength tables
57// (GENTLE_WEATHERING / RIVERINE / FLASH_FLOOD / SUSTAINED) remain.
58//
59// genealogy_id: mei_2007_fast_hydraulic_erosion_gpu +
60// benes_2002_visual_simulation_of_hydraulic_erosion
61// lineage_id: nx_water_erosion_mei2007_full_q14_v2
62//
63// nx_safety_envelope:
64// intended_use: "Hydraulic water-erosion sim (Mei 2007 full
65// cellular automaton) -- deep-time terrain
66// sculpting"
67// sil_target: SIL1
68// asil_target: QM
69// dal_target: NONE
70// evidence: [Mei_2007_canonical_basis, Q14_fixed_point,
71// cellular_automaton_classical]
72// hazard_register: [bug-tape-flow-direction-mass-imbalance,
73// bug-tape-sediment-suspension-overflow]
74// residual_risk: "Procgen quality only."
75// verdict: NOT_YET_EVALUATED
76
77import "nx_syscalls.nx"
78import "nx_tier.nx"
79const NX_MAGIC_2000: i64 = 2000
80
81// ===== Q14 ==========================================================
82const NX_WATER_Q: nx_int = 16384
83
84// ===== Per-cell state layout ========================================
85// 9 i64 per cell. Index pattern: cell(x, y) base offset = (y * w + x) * 9.
86const NX_WATER_STRIDE: nx_int = 9
87const NX_WATER_OFF_H: nx_int = 0 // terrain height (signed metres, Q14)
88const NX_WATER_OFF_D: nx_int = 1 // water depth (metres, Q14)
89const NX_WATER_OFF_S: nx_int = 2 // suspended sediment depth (metres, Q14)
90const NX_WATER_OFF_FL: nx_int = 3 // outflow flux LEFT (m^3/s, Q14)
91const NX_WATER_OFF_FR: nx_int = 4 // outflow flux RIGHT
92const NX_WATER_OFF_FT: nx_int = 5 // outflow flux TOP (lower y)
93const NX_WATER_OFF_FB: nx_int = 6 // outflow flux BOTTOM (higher y)
94const NX_WATER_OFF_U: nx_int = 7 // x-velocity (m/s, Q14)
95const NX_WATER_OFF_V: nx_int = 8 // y-velocity (m/s, Q14)
96
97// ===== Default Mei 2007 parameters =================================
98// All in Q14. Caller can override per nx_water_erosion_params_default.
99// These are the constants Mei 2007 recommends for the visualisation;
100// real-world erosion uses much smaller K_s / K_d.
101const NX_WATER_DEFAULT_GRAVITY: nx_int = 160760 // 9.81 m/s^2
102const NX_WATER_DEFAULT_PIPE_AREA: nx_int = 16384 // 1.0 m^2
103const NX_WATER_DEFAULT_PIPE_LENGTH: nx_int = 16384 // 1.0 m
104const NX_WATER_DEFAULT_DT: nx_int = 1638 // 0.1 seconds per tick
105const NX_WATER_DEFAULT_K_CAPACITY: nx_int = 1638 // 0.1 (Mei's K_c)
106const NX_WATER_DEFAULT_K_SOIL: nx_int = 4915 // 0.3
107const NX_WATER_DEFAULT_K_DEPOSIT: nx_int = 4915 // 0.3
108const NX_WATER_DEFAULT_K_EVAP: nx_int = 164 // 0.01
109const NX_WATER_DEFAULT_CELL_SIZE: nx_int = 16384 // 1.0 m grid cell
110
111// Parameter-pack stride (caller passes 9 i64).
112const NX_WATER_PARAM_COUNT: nx_int = 9
113const NX_WATER_PARAM_GRAVITY: nx_int = 0
114const NX_WATER_PARAM_PIPE_AREA: nx_int = 1
115const NX_WATER_PARAM_PIPE_LENGTH: nx_int = 2
116const NX_WATER_PARAM_DT: nx_int = 3
117const NX_WATER_PARAM_K_CAPACITY: nx_int = 4
118const NX_WATER_PARAM_K_SOIL: nx_int = 5
119const NX_WATER_PARAM_K_DEPOSIT: nx_int = 6
120const NX_WATER_PARAM_K_EVAP: nx_int = 7
121const NX_WATER_PARAM_CELL_SIZE: nx_int = 8
122
123// ===== V1 mode sealed enum (preserved for stateless API) ===========
124const NX_WATER_MODE_GENTLE_WEATHERING: nx_int = 0
125const NX_WATER_MODE_RIVERINE: nx_int = 1
126const NX_WATER_MODE_FLASH_FLOOD: nx_int = 2
127const NX_WATER_MODE_SUSTAINED: nx_int = 3
128
129const NX_WATER_MODE_COUNT: nx_int = 4
130
131func nx_water_mode_is_valid(m: nx_int) -> nx_int {
132 if m == NX_WATER_MODE_GENTLE_WEATHERING { return 1 }
133 if m == NX_WATER_MODE_RIVERINE { return 1 }
134 if m == NX_WATER_MODE_FLASH_FLOOD { return 1 }
135 if m == NX_WATER_MODE_SUSTAINED { return 1 }
136 return 0
137}
138
139// ===== V1 mode params (preserved) ==================================
140func nx_water_mode_params(mode: nx_int, out: *i64) {
141 out[0] = NX_WATER_Q / 2
142 out[1] = NX_WATER_Q / 8
143 out[2] = 500
144 if mode == NX_WATER_MODE_GENTLE_WEATHERING {
145 out[0] = NX_WATER_Q / 16
146 out[1] = NX_WATER_Q / 64
147 out[2] = 50
148 }
149 if mode == NX_WATER_MODE_FLASH_FLOOD {
150 out[0] = 2 * NX_WATER_Q
151 out[1] = NX_WATER_Q / 2
152 out[2] = NX_MAGIC_2000
153 }
154 if mode == NX_WATER_MODE_SUSTAINED {
155 out[0] = NX_WATER_Q
156 out[1] = NX_WATER_Q / 32
157 out[2] = 1000
158 }
159}
160
161// ===== Default parameter pack ======================================
162// Writes the recommended Mei 2007 defaults into params (9 i64).
163func nx_water_erosion_params_default(params: *i64) {
164 params[NX_WATER_PARAM_GRAVITY] = NX_WATER_DEFAULT_GRAVITY
165 params[NX_WATER_PARAM_PIPE_AREA] = NX_WATER_DEFAULT_PIPE_AREA
166 params[NX_WATER_PARAM_PIPE_LENGTH] = NX_WATER_DEFAULT_PIPE_LENGTH
167 params[NX_WATER_PARAM_DT] = NX_WATER_DEFAULT_DT
168 params[NX_WATER_PARAM_K_CAPACITY] = NX_WATER_DEFAULT_K_CAPACITY
169 params[NX_WATER_PARAM_K_SOIL] = NX_WATER_DEFAULT_K_SOIL
170 params[NX_WATER_PARAM_K_DEPOSIT] = NX_WATER_DEFAULT_K_DEPOSIT
171 params[NX_WATER_PARAM_K_EVAP] = NX_WATER_DEFAULT_K_EVAP
172 params[NX_WATER_PARAM_CELL_SIZE] = NX_WATER_DEFAULT_CELL_SIZE
173}
174
175// ===== Cell-index helpers ===========================================
176func _w_cell_base(x: nx_int, y: nx_int, width: nx_int) -> nx_int {
177 return (y * width + x) * NX_WATER_STRIDE
178}
179
180// ===== Init: zero water + sediment + flow ==========================
181// Caller has already filled in heightmap (NX_WATER_OFF_H) per cell.
182// This zeros all other fields.
183func nx_water_erosion_init(state: *i64, width: nx_int, gridh: nx_int) {
184 var i: nx_int = 0
185 let total: nx_int = width * gridh
186 while i < total {
187 let base: nx_int = i * NX_WATER_STRIDE
188 state[base + NX_WATER_OFF_D] = 0
189 state[base + NX_WATER_OFF_S] = 0
190 state[base + NX_WATER_OFF_FL] = 0
191 state[base + NX_WATER_OFF_FR] = 0
192 state[base + NX_WATER_OFF_FT] = 0
193 state[base + NX_WATER_OFF_FB] = 0
194 state[base + NX_WATER_OFF_U] = 0
195 state[base + NX_WATER_OFF_V] = 0
196 i = i + 1
197 }
198}
199
200// ===== Stage 1: rain ===============================================
201// Adds `rain_amount_q14` to every cell's water depth.
202func nx_water_erosion_add_rain(
203 state: *i64, width: nx_int, gridh: nx_int, rain_amount_q14: nx_int
204) {
205 var i: nx_int = 0
206 let total: nx_int = width * gridh
207 while i < total {
208 let base: nx_int = i * NX_WATER_STRIDE
209 state[base + NX_WATER_OFF_D] = state[base + NX_WATER_OFF_D] + rain_amount_q14
210 i = i + 1
211 }
212}
213
214// ===== Stage 2: flow simulation ====================================
215// Pipe model. For each cell, compute outflow in each cardinal
216// direction based on height + water-level difference to that
217// neighbour. Then scale outflows so total doesn't exceed available
218// water * cell_area.
219func nx_water_erosion_compute_flow(
220 state: *i64,
221 width: nx_int,
222 gridh: nx_int,
223 params: *i64
224) {
225 let q: nx_int = NX_WATER_Q
226 let g: nx_int = params[NX_WATER_PARAM_GRAVITY]
227 let a: nx_int = params[NX_WATER_PARAM_PIPE_AREA]
228 let l: nx_int = params[NX_WATER_PARAM_PIPE_LENGTH]
229 let dt: nx_int = params[NX_WATER_PARAM_DT]
230 let cs: nx_int = params[NX_WATER_PARAM_CELL_SIZE]
231 let cs_sq: nx_int = cs * cs / q
232
233 // a * g * dt / l = pipe-model flow rate scaler. In Q14:
234 // scaler = a * g * dt / l (still Q14)
235 let scaler_partial: nx_int = a * g / q * dt / q
236 var scaler: nx_int = scaler_partial
237 if l != 0 { scaler = scaler_partial * q / l }
238
239 var y: nx_int = 0
240 while y < gridh {
241 var x: nx_int = 0
242 while x < width {
243 let base: nx_int = _w_cell_base(x, y, width)
244 let h_here: nx_int = state[base + NX_WATER_OFF_H]
245 let d_here: nx_int = state[base + NX_WATER_OFF_D]
246 let level_here: nx_int = h_here + d_here
247
248 // For each cardinal direction, compute pressure-driven
249 // outflow. Bordering edge -> no outflow that way.
250
251 // LEFT (-x)
252 var f_l_new: nx_int = state[base + NX_WATER_OFF_FL]
253 if x > 0 {
254 let nbase: nx_int = _w_cell_base(x - 1, y, width)
255 let lvl_n: nx_int = state[nbase + NX_WATER_OFF_H] + state[nbase + NX_WATER_OFF_D]
256 let dh: nx_int = level_here - lvl_n
257 f_l_new = f_l_new + scaler * dh / q
258 if f_l_new < 0 { f_l_new = 0 }
259 } else { f_l_new = 0 }
260
261 // RIGHT (+x)
262 var f_r_new: nx_int = state[base + NX_WATER_OFF_FR]
263 if x + 1 < width {
264 let nbase: nx_int = _w_cell_base(x + 1, y, width)
265 let lvl_n: nx_int = state[nbase + NX_WATER_OFF_H] + state[nbase + NX_WATER_OFF_D]
266 let dh: nx_int = level_here - lvl_n
267 f_r_new = f_r_new + scaler * dh / q
268 if f_r_new < 0 { f_r_new = 0 }
269 } else { f_r_new = 0 }
270
271 // TOP (-y)
272 var f_t_new: nx_int = state[base + NX_WATER_OFF_FT]
273 if y > 0 {
274 let nbase: nx_int = _w_cell_base(x, y - 1, width)
275 let lvl_n: nx_int = state[nbase + NX_WATER_OFF_H] + state[nbase + NX_WATER_OFF_D]
276 let dh: nx_int = level_here - lvl_n
277 f_t_new = f_t_new + scaler * dh / q
278 if f_t_new < 0 { f_t_new = 0 }
279 } else { f_t_new = 0 }
280
281 // BOTTOM (+y)
282 var f_b_new: nx_int = state[base + NX_WATER_OFF_FB]
283 if y + 1 < gridh {
284 let nbase: nx_int = _w_cell_base(x, y + 1, width)
285 let lvl_n: nx_int = state[nbase + NX_WATER_OFF_H] + state[nbase + NX_WATER_OFF_D]
286 let dh: nx_int = level_here - lvl_n
287 f_b_new = f_b_new + scaler * dh / q
288 if f_b_new < 0 { f_b_new = 0 }
289 } else { f_b_new = 0 }
290
291 // Outflow scaling. Available water = d_here * cell_size^2.
292 // Total outflow over dt = (f_L + f_R + f_T + f_B) * dt.
293 let total_out: nx_int = f_l_new + f_r_new + f_t_new + f_b_new
294 if total_out > 0 {
295 let available: nx_int = d_here * cs_sq / q
296 let demanded: nx_int = total_out * dt / q
297 if demanded > available {
298 // Scale factor = available / demanded; in Q14.
299 let k: nx_int = available * q / demanded
300 f_l_new = f_l_new * k / q
301 f_r_new = f_r_new * k / q
302 f_t_new = f_t_new * k / q
303 f_b_new = f_b_new * k / q
304 }
305 }
306
307 state[base + NX_WATER_OFF_FL] = f_l_new
308 state[base + NX_WATER_OFF_FR] = f_r_new
309 state[base + NX_WATER_OFF_FT] = f_t_new
310 state[base + NX_WATER_OFF_FB] = f_b_new
311
312 x = x + 1
313 }
314 y = y + 1
315 }
316}
317
318// ===== Stage 3 + 4: water update + velocity ========================
319// Mass-conservation update of d using neighbour-flow. Then compute
320// per-cell velocity (u, v) from the surrounding flow imbalance.
321func nx_water_erosion_update_water_and_velocity(
322 state: *i64,
323 width: nx_int,
324 gridh: nx_int,
325 params: *i64
326) {
327 let q: nx_int = NX_WATER_Q
328 let dt: nx_int = params[NX_WATER_PARAM_DT]
329 let cs: nx_int = params[NX_WATER_PARAM_CELL_SIZE]
330 let cs_sq: nx_int = cs * cs / q
331
332 var y: nx_int = 0
333 while y < gridh {
334 var x: nx_int = 0
335 while x < width {
336 let base: nx_int = _w_cell_base(x, y, width)
337 let f_l_here: nx_int = state[base + NX_WATER_OFF_FL]
338 let f_r_here: nx_int = state[base + NX_WATER_OFF_FR]
339 let f_t_here: nx_int = state[base + NX_WATER_OFF_FT]
340 let f_b_here: nx_int = state[base + NX_WATER_OFF_FB]
341 let total_out: nx_int = f_l_here + f_r_here + f_t_here + f_b_here
342
343 // Inflow from each neighbour's flux toward us.
344 var inflow: nx_int = 0
345 if x > 0 {
346 let nbase: nx_int = _w_cell_base(x - 1, y, width)
347 inflow = inflow + state[nbase + NX_WATER_OFF_FR]
348 }
349 if x + 1 < width {
350 let nbase: nx_int = _w_cell_base(x + 1, y, width)
351 inflow = inflow + state[nbase + NX_WATER_OFF_FL]
352 }
353 if y > 0 {
354 let nbase: nx_int = _w_cell_base(x, y - 1, width)
355 inflow = inflow + state[nbase + NX_WATER_OFF_FB]
356 }
357 if y + 1 < gridh {
358 let nbase: nx_int = _w_cell_base(x, y + 1, width)
359 inflow = inflow + state[nbase + NX_WATER_OFF_FT]
360 }
361
362 // d(t+dt) = d(t) + dt * (inflow - outflow) / cell_area
363 let d_old: nx_int = state[base + NX_WATER_OFF_D]
364 let net: nx_int = (inflow - total_out) * dt / q
365 var d_new: nx_int = d_old + (net * q / cs_sq)
366 if d_new < 0 { d_new = 0 } // physical floor; shouldn't drop below 0
367
368 // Avg depth for velocity normalisation (Mei eq.).
369 let d_avg: nx_int = (d_old + d_new) / 2
370
371 // Velocity: dW_x = (f_R(left_nbr) - f_L(this) + f_R(this) -
372 // f_L(right_nbr)) / 2.
373 var u: nx_int = 0
374 if d_avg > 0 {
375 var f_r_left: nx_int = 0
376 var f_l_right: nx_int = 0
377 if x > 0 {
378 let nbase: nx_int = _w_cell_base(x - 1, y, width)
379 f_r_left = state[nbase + NX_WATER_OFF_FR]
380 }
381 if x + 1 < width {
382 let nbase: nx_int = _w_cell_base(x + 1, y, width)
383 f_l_right = state[nbase + NX_WATER_OFF_FL]
384 }
385 let delta_w_x: nx_int = (f_r_left - f_l_here + f_r_here - f_l_right) / 2
386 u = delta_w_x * q / (cs * d_avg / q)
387 }
388 var v: nx_int = 0
389 if d_avg > 0 {
390 var f_b_top: nx_int = 0
391 var f_t_bottom: nx_int = 0
392 if y > 0 {
393 let nbase: nx_int = _w_cell_base(x, y - 1, width)
394 f_b_top = state[nbase + NX_WATER_OFF_FB]
395 }
396 if y + 1 < gridh {
397 let nbase: nx_int = _w_cell_base(x, y + 1, width)
398 f_t_bottom = state[nbase + NX_WATER_OFF_FT]
399 }
400 let delta_w_y: nx_int = (f_b_top - f_t_here + f_b_here - f_t_bottom) / 2
401 v = delta_w_y * q / (cs * d_avg / q)
402 }
403
404 state[base + NX_WATER_OFF_D] = d_new
405 state[base + NX_WATER_OFF_U] = u
406 state[base + NX_WATER_OFF_V] = v
407
408 x = x + 1
409 }
410 y = y + 1
411 }
412}
413
414// ===== Stage 5: erosion + deposition ==============================
415// Per-cell sediment capacity from velocity + local slope; erode or
416// deposit accordingly.
417func nx_water_erosion_terrain_update(
418 state: *i64,
419 width: nx_int,
420 gridh: nx_int,
421 params: *i64
422) {
423 let q: nx_int = NX_WATER_Q
424 let kc: nx_int = params[NX_WATER_PARAM_K_CAPACITY]
425 let ks: nx_int = params[NX_WATER_PARAM_K_SOIL]
426 let kd: nx_int = params[NX_WATER_PARAM_K_DEPOSIT]
427 let cs: nx_int = params[NX_WATER_PARAM_CELL_SIZE]
428
429 var y: nx_int = 0
430 while y < gridh {
431 var x: nx_int = 0
432 while x < width {
433 let base: nx_int = _w_cell_base(x, y, width)
434 let u: nx_int = state[base + NX_WATER_OFF_U]
435 let v: nx_int = state[base + NX_WATER_OFF_V]
436 // |velocity| approx = max(|u|, |v|) + 0.5 * min (Manhattan-
437 // octagonal cheap-norm; close enough for capacity sizing).
438 var u_abs: nx_int = u
439 if u_abs < 0 { u_abs = 0 - u_abs }
440 var v_abs: nx_int = v
441 if v_abs < 0 { v_abs = 0 - v_abs }
442 var vel_mag: nx_int = u_abs
443 if v_abs > vel_mag {
444 vel_mag = v_abs + u_abs / 2
445 } else {
446 vel_mag = u_abs + v_abs / 2
447 }
448
449 // Local slope sin (cheap approx): max neighbour height
450 // delta / cell size.
451 var max_dh: nx_int = 0
452 let h_here: nx_int = state[base + NX_WATER_OFF_H]
453 if x > 0 {
454 let nbase: nx_int = _w_cell_base(x - 1, y, width)
455 var dh: nx_int = h_here - state[nbase + NX_WATER_OFF_H]
456 if dh < 0 { dh = 0 - dh }
457 if dh > max_dh { max_dh = dh }
458 }
459 if x + 1 < width {
460 let nbase: nx_int = _w_cell_base(x + 1, y, width)
461 var dh: nx_int = h_here - state[nbase + NX_WATER_OFF_H]
462 if dh < 0 { dh = 0 - dh }
463 if dh > max_dh { max_dh = dh }
464 }
465 if y > 0 {
466 let nbase: nx_int = _w_cell_base(x, y - 1, width)
467 var dh: nx_int = h_here - state[nbase + NX_WATER_OFF_H]
468 if dh < 0 { dh = 0 - dh }
469 if dh > max_dh { max_dh = dh }
470 }
471 if y + 1 < gridh {
472 let nbase: nx_int = _w_cell_base(x, y + 1, width)
473 var dh: nx_int = h_here - state[nbase + NX_WATER_OFF_H]
474 if dh < 0 { dh = 0 - dh }
475 if dh > max_dh { max_dh = dh }
476 }
477 // sin(slope) ~ dh / cell_size (small-angle). Floor at
478 // ~0.05 so even flat ground has SOME capacity (Mei's lmin).
479 var slope_sin: nx_int = max_dh * q / cs
480 if slope_sin < q / 20 { slope_sin = q / 20 }
481
482 // Capacity = K_c * slope_sin * |velocity|.
483 let capacity: nx_int = kc * slope_sin / q * vel_mag / q
484
485 let s_here: nx_int = state[base + NX_WATER_OFF_S]
486 if capacity > s_here {
487 // Erode.
488 let delta: nx_int = (capacity - s_here) * ks / q
489 state[base + NX_WATER_OFF_H] = h_here - delta
490 state[base + NX_WATER_OFF_S] = s_here + delta
491 } else {
492 // Deposit.
493 let delta: nx_int = (s_here - capacity) * kd / q
494 state[base + NX_WATER_OFF_H] = h_here + delta
495 state[base + NX_WATER_OFF_S] = s_here - delta
496 }
497
498 x = x + 1
499 }
500 y = y + 1
501 }
502}
503
504// ===== Stage 6: sediment transport (semi-Lagrangian) ===============
505// Backtrack from (x, y) along (-u, -v) and sample sediment via
506// bilinear interpolation. Writes new sediment into a CALLER-PROVIDED
507// scratch buffer (not in-place, to avoid order-dependence) then
508// copies back.
509func nx_water_erosion_transport(
510 state: *i64,
511 scratch_s: *i64,
512 width: nx_int,
513 gridh: nx_int,
514 params: *i64
515) {
516 let q: nx_int = NX_WATER_Q
517 let dt: nx_int = params[NX_WATER_PARAM_DT]
518 let cs: nx_int = params[NX_WATER_PARAM_CELL_SIZE]
519
520 var y: nx_int = 0
521 while y < gridh {
522 var x: nx_int = 0
523 while x < width {
524 let base: nx_int = _w_cell_base(x, y, width)
525 let u: nx_int = state[base + NX_WATER_OFF_U]
526 let v: nx_int = state[base + NX_WATER_OFF_V]
527 // Backtrack position in cell coords: px = x - u*dt/cs.
528 let px_q: nx_int = x * q - u * dt / cs
529 let py_q: nx_int = y * q - v * dt / cs
530 // Integer + fractional decomposition (Q14).
531 var px_i: nx_int = px_q / q
532 var py_i: nx_int = py_q / q
533 if px_i < 0 { px_i = 0 }
534 if py_i < 0 { py_i = 0 }
535 if px_i >= width - 1 { px_i = width - 2 }
536 if py_i >= gridh - 1 { py_i = gridh - 2 }
537 if px_i < 0 { px_i = 0 }
538 if py_i < 0 { py_i = 0 }
539 let fx: nx_int = px_q - px_i * q
540 let fy: nx_int = py_q - py_i * q
541 var fx_c: nx_int = fx
542 if fx_c < 0 { fx_c = 0 }
543 if fx_c > q { fx_c = q }
544 var fy_c: nx_int = fy
545 if fy_c < 0 { fy_c = 0 }
546 if fy_c > q { fy_c = q }
547
548 // Bilinear sample sediment.
549 let s00: nx_int = state[_w_cell_base(px_i, py_i, width) + NX_WATER_OFF_S]
550 let s10: nx_int = state[_w_cell_base(px_i + 1, py_i, width) + NX_WATER_OFF_S]
551 let s01: nx_int = state[_w_cell_base(px_i, py_i + 1, width) + NX_WATER_OFF_S]
552 let s11: nx_int = state[_w_cell_base(px_i + 1, py_i + 1, width) + NX_WATER_OFF_S]
553 let s0: nx_int = s00 + (s10 - s00) * fx_c / q
554 let s1: nx_int = s01 + (s11 - s01) * fx_c / q
555 let sample: nx_int = s0 + (s1 - s0) * fy_c / q
556 scratch_s[y * width + x] = sample
557
558 x = x + 1
559 }
560 y = y + 1
561 }
562 // Copy scratch back into state.
563 var i: nx_int = 0
564 let total: nx_int = width * gridh
565 while i < total {
566 let base: nx_int = i * NX_WATER_STRIDE
567 state[base + NX_WATER_OFF_S] = scratch_s[i]
568 i = i + 1
569 }
570}
571
572// ===== Stage 7: evaporation =========================================
573func nx_water_erosion_evaporate(
574 state: *i64, width: nx_int, gridh: nx_int, params: *i64
575) {
576 let q: nx_int = NX_WATER_Q
577 let ke: nx_int = params[NX_WATER_PARAM_K_EVAP]
578 let dt: nx_int = params[NX_WATER_PARAM_DT]
579 // factor = 1 - K_e * dt. In Q14.
580 let factor: nx_int = q - ke * dt / q
581 var i: nx_int = 0
582 let total: nx_int = width * gridh
583 while i < total {
584 let base: nx_int = i * NX_WATER_STRIDE
585 state[base + NX_WATER_OFF_D] = state[base + NX_WATER_OFF_D] * factor / q
586 i = i + 1
587 }
588}
589
590// ===== Composed single tick ========================================
591// Runs all 7 stages in sequence per Mei 2007. Caller still drives
592// rain (call nx_water_erosion_add_rain before this, or just embed it
593// in the loop with this).
594func nx_water_erosion_tick(
595 state: *i64,
596 scratch_s: *i64,
597 width: nx_int,
598 gridh: nx_int,
599 rain_amount_q14: nx_int,
600 params: *i64
601) {
602 nx_water_erosion_add_rain(state, width, gridh, rain_amount_q14)
603 nx_water_erosion_compute_flow(state, width, gridh, params)
604 nx_water_erosion_update_water_and_velocity(state, width, gridh, params)
605 nx_water_erosion_terrain_update(state, width, gridh, params)
606 nx_water_erosion_transport(state, scratch_s, width, gridh, params)
607 nx_water_erosion_evaporate(state, width, gridh, params)
608}
609
610// ===== V1 stateless API (preserved) ================================
611// Quick-check primitive: given a single slope value, return per-mode
612// erosion delta. Use when full Mei 2007 simulation is overkill (e.g.
613// heightmap touch-up, isolated grading).
614func nx_water_erosion_delta(
615 slope_magnitude_q14: nx_int,
616 mode: nx_int
617) -> nx_int {
618 if slope_magnitude_q14 <= 0 { return 0 }
619 if nx_water_mode_is_valid(mode) == 0 { return 0 }
620 let params: *i64 = (sys_mmap(3 * NX_SIZEOF_NX_INT)) as *i64
621 nx_water_mode_params(mode, params)
622 let strength: nx_int = params[0]
623 let threshold: nx_int = params[1]
624 let max_delta: nx_int = params[2]
625 if slope_magnitude_q14 <= threshold { return 0 }
626 let excess: nx_int = slope_magnitude_q14 - threshold
627 var delta: nx_int = 0 - strength * excess / NX_WATER_Q
628 if delta < 0 - max_delta { delta = 0 - max_delta }
629 return delta
630}
631
632func nx_water_slope_at(
633 heightmap: *i64,
634 width: nx_int,
635 grid_height: nx_int,
636 x: nx_int,
637 y: nx_int
638) -> nx_int {
639 if x < 0 { return 0 }
640 if y < 0 { return 0 }
641 if x >= width { return 0 }
642 if y >= grid_height { return 0 }
643 let h_here: nx_int = heightmap[y * width + x]
644 var h_min: nx_int = h_here
645 if x > 0 {
646 let h_w: nx_int = heightmap[y * width + (x - 1)]
647 if h_w < h_min { h_min = h_w }
648 }
649 if x + 1 < width {
650 let h_e: nx_int = heightmap[y * width + (x + 1)]
651 if h_e < h_min { h_min = h_e }
652 }
653 if y > 0 {
654 let h_n: nx_int = heightmap[(y - 1) * width + x]
655 if h_n < h_min { h_min = h_n }
656 }
657 if y + 1 < grid_height {
658 let h_s: nx_int = heightmap[(y + 1) * width + x]
659 if h_s < h_min { h_min = h_s }
660 }
661 return h_here - h_min
662}
663
664func nx_water_erosion_at(
665 heightmap: *i64,
666 width: nx_int,
667 grid_height: nx_int,
668 x: nx_int,
669 y: nx_int,
670 mode: nx_int
671) -> nx_int {
672 let slope: nx_int = nx_water_slope_at(heightmap, width, grid_height, x, y)
673 return nx_water_erosion_delta(slope, mode)
674}
675
676// ===== Self-test ====================================================
677func main() -> i64 {
678 let q: nx_int = NX_WATER_Q
679
680 // T1: V1 mode validity preserved.
681 if nx_water_mode_is_valid(NX_WATER_MODE_GENTLE_WEATHERING) != 1 { return __syscall(93, 1, 0, 0, 0, 0, 0) }
682 if nx_water_mode_is_valid(NX_WATER_MODE_FLASH_FLOOD) != 1 { return __syscall(93, 2, 0, 0, 0, 0, 0) }
683 if nx_water_mode_is_valid(99) != 0 { return __syscall(93, 3, 0, 0, 0, 0, 0) }
684
685 // T2: V1 delta sanity. Flat -> 0. Above-threshold -> negative.
686 if nx_water_erosion_delta(0, NX_WATER_MODE_RIVERINE) != 0 { return __syscall(93, 4, 0, 0, 0, 0, 0) }
687 if nx_water_erosion_delta(q, NX_WATER_MODE_RIVERINE) >= 0 { return __syscall(93, 5, 0, 0, 0, 0, 0) }
688
689 // T3: Default-param pack.
690 let params: *i64 = (sys_mmap(NX_WATER_PARAM_COUNT * NX_SIZEOF_NX_INT)) as *i64
691 nx_water_erosion_params_default(params)
692 if params[NX_WATER_PARAM_GRAVITY] != NX_WATER_DEFAULT_GRAVITY { return __syscall(93, 10, 0, 0, 0, 0, 0) }
693 if params[NX_WATER_PARAM_K_EVAP] != NX_WATER_DEFAULT_K_EVAP { return __syscall(93, 11, 0, 0, 0, 0, 0) }
694
695 // T4: State init + zero fields. Use a small 4x4 grid.
696 let w: nx_int = 4
697 let h: nx_int = 4
698 let n: nx_int = w * h
699 let state: *i64 = (sys_mmap(n * NX_WATER_STRIDE * NX_SIZEOF_NX_INT)) as *i64
700 // Caller pre-fills heightmap. We'll write a slope from x=0 (high)
701 // to x=3 (low) -- water should flow rightward.
702 var i: nx_int = 0
703 while i < n {
704 let xi: nx_int = i % w
705 let base: nx_int = i * NX_WATER_STRIDE
706 state[base + NX_WATER_OFF_H] = (3 - xi) * 100 * q // slope down to +x
707 // Other fields random non-zero to verify init zeros them.
708 state[base + NX_WATER_OFF_D] = 99
709 state[base + NX_WATER_OFF_S] = 88
710 state[base + NX_WATER_OFF_FL] = 77
711 i = i + 1
712 }
713 nx_water_erosion_init(state, w, h)
714 // After init: D, S, F* all zero.
715 var k: nx_int = 0
716 while k < n {
717 let base: nx_int = k * NX_WATER_STRIDE
718 if state[base + NX_WATER_OFF_D] != 0 { return __syscall(93, 20, 0, 0, 0, 0, 0) }
719 if state[base + NX_WATER_OFF_S] != 0 { return __syscall(93, 21, 0, 0, 0, 0, 0) }
720 if state[base + NX_WATER_OFF_FL] != 0 { return __syscall(93, 22, 0, 0, 0, 0, 0) }
721 if state[base + NX_WATER_OFF_FR] != 0 { return __syscall(93, 23, 0, 0, 0, 0, 0) }
722 // Heightmap preserved.
723 let xi: nx_int = k % w
724 if state[base + NX_WATER_OFF_H] != (3 - xi) * 100 * q { return __syscall(93, 24, 0, 0, 0, 0, 0) }
725 k = k + 1
726 }
727
728 // T5: Rain adds to every cell uniformly.
729 nx_water_erosion_add_rain(state, w, h, 5 * q)
730 var k2: nx_int = 0
731 while k2 < n {
732 let base: nx_int = k2 * NX_WATER_STRIDE
733 if state[base + NX_WATER_OFF_D] != 5 * q { return __syscall(93, 30, 0, 0, 0, 0, 0) }
734 k2 = k2 + 1
735 }
736
737 // T6: Compute-flow produces rightward outflow on a leftward-high
738 // slope. After flow + water update, cell (0, 0) should have less
739 // water than (3, 0); right-flowing flux out of (0, 0) should be > 0.
740 nx_water_erosion_compute_flow(state, w, h, params)
741 let base_00: nx_int = _w_cell_base(0, 0, w)
742 let base_30: nx_int = _w_cell_base(3, 0, w)
743 if state[base_00 + NX_WATER_OFF_FR] <= 0 { return __syscall(93, 40, 0, 0, 0, 0, 0) }
744 // Cell (3, 0) is the lowest; should have NO right-outflow (edge).
745 if state[base_30 + NX_WATER_OFF_FR] != 0 { return __syscall(93, 41, 0, 0, 0, 0, 0) }
746
747 // T7: Full tick runs end-to-end without crashing.
748 let scratch: *i64 = (sys_mmap(n * NX_SIZEOF_NX_INT)) as *i64
749 nx_water_erosion_tick(state, scratch, w, h, q, params)
750 // After 1 tick on a slope, the highest cell should have lower
751 // water than the lowest (water flows downhill).
752 if state[base_00 + NX_WATER_OFF_D] >= state[base_30 + NX_WATER_OFF_D] {
753 return __syscall(93, 50, 0, 0, 0, 0, 0)
754 }
755
756 // T8: Evaporation reduces water without going negative.
757 let pre_d: nx_int = state[base_00 + NX_WATER_OFF_D]
758 nx_water_erosion_evaporate(state, w, h, params)
759 let post_d: nx_int = state[base_00 + NX_WATER_OFF_D]
760 if post_d > pre_d { return __syscall(93, 60, 0, 0, 0, 0, 0) }
761 if post_d < 0 { return __syscall(93, 61, 0, 0, 0, 0, 0) }
762
763 // T9: Multi-tick stability -- 10 ticks of simulation, total water
764 // should be bounded (rain in + evaporate ~~ steady-state).
765 var t: nx_int = 0
766 while t < 10 {
767 nx_water_erosion_tick(state, scratch, w, h, q / 10, params)
768 t = t + 1
769 }
770 // Just check no NaN/overflow: every cell's water in reasonable
771 // range.
772 var k3: nx_int = 0
773 while k3 < n {
774 let base: nx_int = k3 * NX_WATER_STRIDE
775 if state[base + NX_WATER_OFF_D] < 0 { return __syscall(93, 70, 0, 0, 0, 0, 0) }
776 if state[base + NX_WATER_OFF_D] > 1000 * q { return __syscall(93, 71, 0, 0, 0, 0, 0) }
777 k3 = k3 + 1
778 }
779
780 return 0
781}