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