nx_fluid.nx source

1// nx_fluid.nx -- ★R6 of the Infinigen ladder: FLUID SIMULATION (their FLIP-fluids axis). A height-field WAVE 2// simulation (the shallow-water / wave-equation grid -- the standard real-time water technique; FLIP is the 3// particle variant). Real physics, all integer: a height grid h + velocity grid v integrated by the discrete 4// wave equation v += c^2*dt * laplacian(h) ; v *= damp ; h += v , with REFLECTING boundaries so VOLUME is 5// exactly conserved (sum(laplacian)=0 => sum(v) stays 0 => sum(h) constant). Disturbances propagate outward at 6// finite speed, reflect off shores, and damp -- not a flat plane. license_tier: ORIGINAL 7import "nx_syscalls.nx" 8import "nx_itrig.nx" 9import "nx_vecmath.nx" 10const FL_MAGIC_1024: i64 = 1024 11const FL_MAGIC_1048576: i64 = 1048576 12const FL_MAGIC_9973: i64 = 9973 13const FL_MAGIC_4096: i64 = 4096 14const FL_MAGIC_65536: i64 = 65536 15 16const FL_N: i64 = 168 // grid cells per axis 17const FL_K: i64 = 100 // c^2*dt (fx1024) -- CFL-stable (4*K/1024 = 0.39, well under the 2D limit) 18const FL_DAMP: i64 = 1013 // velocity damping per step (~0.989) -- removes energy, keeps sum(v)=0 19static FL_H: i64 // height grid ptr (FL_N*FL_N i64, fx) 20static FL_V: i64 // velocity grid ptr 21static FL_STEP: i64 // step counter (rotates the conservation-remainder offset) 22 23func fl_isqrt(x: i64) -> i64 { return vm_isqrt(x) } 24func fl_iabs(v: i64) -> i64 { if v<0 {return 0-v} return v } 25func fl_clampi(v: i64, lo: i64, hi: i64) -> i64 { if v<lo {return lo} if v>hi {return hi} return v } 26 27func fl_init() -> i64 { 28 if FL_H == 0 { FL_H = sys_mmap(FL_N*FL_N*8) as i64 } 29 if FL_V == 0 { FL_V = sys_mmap(FL_N*FL_N*8) as i64 } 30 let h: *i64 = FL_H as *i64 31 let v: *i64 = FL_V as *i64 32 var i: i64 = 0 33 while i < FL_N*FL_N { h[i]=0; v[i]=0; i=i+1 } 34 FL_STEP = 0 35 return 0 36} 37 38// place a smooth circular disturbance (a raindrop / splash) of amplitude amp at (cx,cy), radius r. 39func fl_drop(cx: i64, cy: i64, r: i64, amp: i64) -> i64 { 40 let h: *i64 = FL_H as *i64 41 var y: i64 = cy - r 42 while y <= cy + r { 43 var x: i64 = cx - r 44 while x <= cx + r { 45 if x >= 0 { if x < FL_N { if y >= 0 { if y < FL_N { 46 let d2: i64 = (x-cx)*(x-cx) + (y-cy)*(y-cy) 47 if d2 < r*r { 48 // cosine-ish bump: amp * (1 + cos(pi*d/r))/2 via (r^2-d^2)/r^2 squared, smooth 49 let f: i64 = (r*r - d2)*FL_MAGIC_1024/(r*r) 50 h[y*FL_N+x] = h[y*FL_N+x] + amp*f*f/FL_MAGIC_1048576 51 } 52 } } } } 53 x = x + 1 54 } 55 y = y + 1 56 } 57 return 0 58} 59 60// one wave-equation step. REFLECTING boundaries (out-of-grid neighbour = centre) => volume-conserving. 61func fl_step() -> i64 { 62 let h: *i64 = FL_H as *i64 63 let v: *i64 = FL_V as *i64 64 // pass 1: v += K * laplacian(h) ; v *= damp (reads h only) 65 var y: i64 = 0 66 while y < FL_N { 67 var x: i64 = 0 68 while x < FL_N { 69 let c: i64 = h[y*FL_N+x] 70 var hl: i64 = c 71 var hr: i64 = c 72 var hu: i64 = c 73 var hd: i64 = c 74 if x > 0 { hl = h[y*FL_N+x-1] } 75 if x < FL_N-1 { hr = h[y*FL_N+x+1] } 76 if y > 0 { hu = h[(y-1)*FL_N+x] } 77 if y < FL_N-1 { hd = h[(y+1)*FL_N+x] } 78 let lap: i64 = hl + hr + hu + hd - 4*c 79 var nv: i64 = v[y*FL_N+x] + lap*FL_K/FL_MAGIC_1024 80 nv = nv*FL_DAMP/FL_MAGIC_1024 81 v[y*FL_N+x] = nv 82 x = x + 1 83 } 84 y = y + 1 85 } 86 // ★conservation projection: integer truncation biases the per-cell v-update so sum(v) drifts off 0 (which 87 // would slowly change total volume). Subtract the mean velocity -> sum(v)=0 exactly (removes the spurious 88 // uniform-drift mode; real waves, whose mean is ~0, are untouched). Guarantees mass conservation. 89 var vs: i64 = 0 90 var k: i64 = 0 91 while k < FL_N*FL_N { vs = vs + v[k]; k = k + 1 } 92 let vm: i64 = vs / (FL_N*FL_N) 93 let rem: i64 = vs - vm*(FL_N*FL_N) // the truncation leftover (same sign as vs) 94 // remove the mean from every cell (sum(v) now == rem) 95 var i: i64 = 0 96 while i < FL_N*FL_N { v[i] = v[i] - vm; i = i + 1 } 97 // zero the leftover EXACTLY by distributing -+1 to |rem| cells -- but at a ROTATING offset each step so the 98 // corrections spread uniformly over time (a fixed offset would pump energy into one region). sum(v)=0 exactly. 99 FL_STEP = FL_STEP + 1 100 if rem != 0 { 101 let off: i64 = (FL_STEP*FL_MAGIC_9973) % (FL_N*FL_N) 102 var cnt: i64 = rem 103 var dir: i64 = 0 - 1 104 if rem < 0 { cnt = 0 - rem; dir = 1 } 105 var c2: i64 = 0 106 while c2 < cnt { 107 let idx: i64 = (off + c2) % (FL_N*FL_N) 108 v[idx] = v[idx] + dir 109 c2 = c2 + 1 110 } 111 } 112 // pass 2: h += v 113 i = 0 114 while i < FL_N*FL_N { h[i] = h[i] + v[i]; i = i + 1 } 115 return 0 116} 117 118func fl_h(x: i64, y: i64) -> i64 { let h: *i64 = FL_H as *i64; return h[y*FL_N+x] } 119func fl_n() -> i64 { return FL_N } 120// total signed volume (conserved) and total abs-energy (damps) 121func fl_volume() -> i64 { let h: *i64 = FL_H as *i64; var s: i64=0; var i: i64=0; while i<FL_N*FL_N { s=s+h[i]; i=i+1 } return s } 122func fl_energy() -> i64 { let h: *i64 = FL_H as *i64; var s: i64=0; var i: i64=0; while i<FL_N*FL_N { s=s+fl_iabs(h[i]); i=i+1 } return s } 123// peak amplitude max|h| -- the clean damping signal (decays from spreading + damping; unlike sum|h| which 124// grows as the wave develops troughs) 125func fl_maxabs() -> i64 { let h: *i64 = FL_H as *i64; var m: i64=0; var i: i64=0; while i<FL_N*FL_N { let a: i64=fl_iabs(h[i]); if a>m {m=a} i=i+1 } return m } 126// wavefront radius from (cx,cy): farthest disturbed cell (|h|>thr) 127func fl_front(cx: i64, cy: i64, thr: i64) -> i64 { 128 let h: *i64 = FL_H as *i64 129 var best: i64 = 0 130 var y: i64 = 0 131 while y < FL_N { 132 var x: i64 = 0 133 while x < FL_N { 134 if fl_iabs(h[y*FL_N+x]) > thr { 135 let d: i64 = fl_isqrt((x-cx)*(x-cx)+(y-cy)*(y-cy)) 136 if d > best { best = d } 137 } 138 x = x + 1 139 } 140 y = y + 1 141 } 142 return best 143} 144 145// render the height field as a shaded water surface into fb (W*H): normal from the height gradient, sun diffuse 146// + specular glint + deep-water base + sky tint. Ripples read as light/dark bands. sun = fx1024 unit dir. 147func fl_render(fb: *i64, W: i64, H: i64, sunx: i64, suny: i64, sunz: i64) -> i64 { 148 let h: *i64 = FL_H as *i64 149 let sl: i64 = fl_isqrt(sunx*sunx+suny*suny+sunz*sunz) 150 let lx: i64 = sunx*FL_MAGIC_1024/sl 151 let ly: i64 = suny*FL_MAGIC_1024/sl 152 let lz: i64 = sunz*FL_MAGIC_1024/sl 153 var py: i64 = 0 154 while py < H { 155 var px: i64 = 0 156 while px < W { 157 let cx: i64 = fl_clampi(px*FL_N/W, 1, FL_N-2) 158 let cy: i64 = fl_clampi(py*FL_N/H, 1, FL_N-2) 159 let dhx: i64 = h[cy*FL_N+cx+1] - h[cy*FL_N+cx-1] 160 let dhy: i64 = h[(cy+1)*FL_N+cx] - h[(cy-1)*FL_N+cx] 161 // surface normal (height scaled down so slopes are gentle): (-dhx, S, -dhy) 162 var nx: i64 = 0 - dhx*4 163 var ny: i64 = FL_MAGIC_4096 164 var nz: i64 = 0 - dhy*4 165 let nl: i64 = fl_isqrt(nx*nx+ny*ny+nz*nz) 166 if nl > 0 { nx=nx*FL_MAGIC_1024/nl; ny=ny*FL_MAGIC_1024/nl; nz=nz*FL_MAGIC_1024/nl } 167 var nd: i64 = (nx*lx+ny*ly+nz*lz)/FL_MAGIC_1024 168 if nd < 0 { nd = 0 } 169 // specular: reflect the view (straight down, 0,1,0) not used; use half-vector with view=up 170 let hx: i64 = lx 171 let hyv: i64 = ly + FL_MAGIC_1024 172 let hz: i64 = lz 173 let hl2: i64 = fl_isqrt(hx*hx+hyv*hyv+hz*hz) 174 var sp: i64 = (nx*hx+ny*hyv+nz*hz)/hl2 175 if sp < 0 { sp = 0 } 176 sp = sp*sp/FL_MAGIC_1024; sp = sp*sp/FL_MAGIC_1024; sp = sp*sp/FL_MAGIC_1024 // ^8 tight glint 177 // deep-water base + sky tint by upward-facing + diffuse 178 let lit: i64 = 40 + nd*70/FL_MAGIC_1024 179 var r: i64 = 26*lit/100 + ny*24/FL_MAGIC_1024 + sp*200/FL_MAGIC_1024 180 var g: i64 = 78*lit/100 + ny*30/FL_MAGIC_1024 + sp*205/FL_MAGIC_1024 181 var b: i64 = 128*lit/100 + ny*36/FL_MAGIC_1024 + sp*210/FL_MAGIC_1024 182 if r>255{r=255} if g>255{g=255} if b>255{b=255} 183 if r<0{r=0} if g<0{g=0} if b<0{b=0} 184 fb[py*W+px] = r + g*256 + b*FL_MAGIC_65536 185 px = px + 1 186 } 187 py = py + 1 188 } 189 return 0 190}