fuzzyworld

bubbles / files / src / physics / Water.ts

Water.ts

7.5 KB · history · blame

1/**
2 * Water surface simulation — 1D Shallow Water Equations (SWE), staggered grid.
3 *
4 * Grid layout:
5 * h[i] column height at cell centre i = 0 … N−1
6 * measured from the "gravity wall" (the wall gravity pulls toward).
7 * u[j] horizontal velocity at face j = 0 … N
8 * u[0] = u[N] = 0 (solid walls, sealed container)
9 *
10 * Governing equations:
11 * ∂h/∂t + H̄ · ∂u/∂x = 0 (mass continuity, nonlinear flux)
12 * ∂u/∂t = −g_sw · ∂h/∂x + aₓ (depth-averaged momentum)
13 *
14 * Full 360° support: when the device is flipped upside-down (verticalG < 0),
15 * h[i] is reinterpreted as measuring from the ceiling rather than the floor.
16 * getSurfaceY() accounts for this so callers need no special-casing.
17 *
18 * The container is sealed — no spilling. Volume is conserved permanently.
19 */
20export class Water {
21 /** Column heights (cell centres). */
22 private readonly h: Float32Array;
23 /** Face velocities. Length = columns + 1; u[0] = u[N] = 0 always. */
24 private readonly u: Float32Array;
25 /** Scratch buffer for continuity update — avoids in-place aliasing errors. */
26 private readonly hNext: Float32Array;
27
28 /** Conserved water volume [px²]. Never changes — sealed container. */
29 private volume: number;
30 /** Whether gravity currently points toward the ceiling (phone upside-down). */
31 private _isFlipped = false;
32 private _lastVerticalSign = 1;
33
34 private readonly gSW: number;
35 private readonly kDamp: number;
36 private static readonly CFL = 0.85;
37
38 constructor(
39 public readonly columns: number = 80,
40 public glassW: number = 300,
41 public glassH: number = 500,
42 fillFrac: number = 0.55,
43 /** Override wave gravity (default 980). Lower = slower/more viscous feel. */
44 gSW: number = 980,
45 /** Override damping coefficient (default 1.5). Higher = more viscous. */
46 kDamp: number = 1.5,
47 ) {
48 this.gSW = gSW;
49 this.kDamp = kDamp;
50 const h0 = glassH * fillFrac;
51 this.h = new Float32Array(columns).fill(h0);
52 this.hNext = new Float32Array(columns);
53 this.u = new Float32Array(columns + 1).fill(0);
54 this.volume = h0 * glassW;
55 }
56
57 get columnWidth(): number { return this.glassW / this.columns; }
58 get averageHeight(): number { return this.volume / this.glassW; }
59 /** True when gravity points toward the ceiling (phone held upside-down). */
60 get isFlipped(): boolean { return this._isFlipped; }
61
62 /**
63 * Height of the water column at index i, measured from the gravity wall
64 * (floor when upright, ceiling when flipped).
65 */
66 getColumnHeight(i: number): number {
67 return this.h[Math.max(0, Math.min(this.columns - 1, i))]!;
68 }
69
70 step(dt: number, lateralG: number, verticalG: number): void {
71 // ── Detect gravity direction flip (phone turned upside-down) ────────
72 const vSign = verticalG >= 0 ? 1 : -1;
73 if (vSign !== this._lastVerticalSign) {
74 this._lastVerticalSign = vSign;
75 this._isFlipped = !this._isFlipped;
76 const h0 = this.volume / this.glassW;
77 this.h.fill(h0);
78 this.u.fill(0);
79 }
80
81 const N = this.columns;
82 const dx = this.glassW / N;
83
84 // ── Adaptive sub-steps: keep CFL = dt_sub/dx * c_max < CFL_target ────
85 let maxH = 1;
86 for (let i = 0; i < N; i++) if (this.h[i]! > maxH) maxH = this.h[i]!;
87 const cMax = Math.sqrt(this.gSW * maxH);
88 const subSteps = Math.max(1, Math.ceil(dt * cMax / (dx * Water.CFL)));
89 const dtS = dt / subSteps;
90 const dampFactor = Math.exp(-this.kDamp * dtS);
91
92 for (let s = 0; s < subSteps; s++) {
93 // ── 1. Momentum: pressure gradient + tilt body-force + damping ────
94 for (let j = 1; j < N; j++) {
95 this.u[j] = (
96 this.u[j]!
97 - this.gSW * (dtS / dx) * (this.h[j]! - this.h[j - 1]!)
98 + lateralG * dtS
99 ) * dampFactor;
100 }
101 // Wall BCs — no flow through side walls
102 this.u[0] = 0;
103 this.u[N] = 0;
104
105 // ── 2. Continuity: upwind nonlinear flux into scratch buffer ──────
106 for (let i = 0; i < N; i++) {
107 const uR = this.u[i + 1]!;
108 const hR = uR >= 0 ? this.h[i]! : (i + 1 < N ? this.h[i + 1]! : this.h[i]!);
109 const uL = this.u[i]!;
110 const hL = uL <= 0 ? this.h[i]! : (i > 0 ? this.h[i - 1]! : this.h[i]!);
111 this.hNext[i] = this.h[i]! - (dtS / dx) * (uR * hR - uL * hL);
112 }
113 this.h.set(this.hNext);
114
115 // ── 3. Clamp: gravity wall and opposite wall (sealed) ─────────────
116 for (let i = 0; i < N; i++) {
117 if (this.h[i]! < 0) this.h[i] = 0;
118 else if (this.h[i]! > this.glassH) this.h[i] = this.glassH;
119 }
120 }
121
122 // ── 4. Volume correction (sealed — no spill, only clamp rounding) ────
123 let sum = 0;
124 for (let i = 0; i < N; i++) sum += this.h[i]!;
125 const actualVol = sum * dx;
126 if (actualVol > 1 && Math.abs(actualVol - this.volume) > 0.5) {
127 const delta = (this.volume - actualVol) / this.glassW;
128 for (let i = 0; i < N; i++) {
129 this.h[i] = Math.max(0, Math.min(this.glassH, this.h[i]! + delta));
130 }
131 }
132 }
133
134 /**
135 * Y coordinate (in glass space: 0 = top, glassH = bottom) of the
136 * water–air interface at horizontal position x.
137 *
138 * When not flipped: surface is above the water (water fills below).
139 * When flipped: surface is below the water (water fills above).
140 * Callers can use `water.isFlipped` to know which side the water is on.
141 */
142 getSurfaceY(x: number): number {
143 const col = Math.max(0, Math.min(this.columns - 1, Math.floor(x / this.columnWidth)));
144 // isFlipped: h[col] measures from ceiling → surface is that far from top
145 // normal: h[col] measures from floor → surface is (glassH − h) from top
146 return this._isFlipped ? this.h[col]! : this.glassH - this.h[col]!;
147 }
148
149 /**
150 * Seed horizontal face velocities near an impact point, creating an
151 * outward-spreading ripple. velocityY > 0 means downward impact.
152 */
153 disturb(localX: number, radius: number, velocityY: number): void {
154 const dx = this.columnWidth;
155 const N = this.columns;
156 const centerFace = localX / dx;
157 const spreadFaces = (radius / dx) * 1.5;
158 const strength = velocityY * 0.03;
159
160 for (let j = 1; j < N; j++) {
161 const dist = j - centerFace;
162 if (Math.abs(dist) > spreadFaces * 2) continue;
163 const falloff = Math.exp(-(dist * dist) / (spreadFaces * spreadFaces));
164 this.u[j] = this.u[j]! + strength * Math.sign(dist) * falloff;
165 }
166 }
167
168 resetFill(fillFrac: number = 0.55): void {
169 const h0 = this.glassH * fillFrac;
170 this.h.fill(h0);
171 this.u.fill(0);
172 this.volume = h0 * this.glassW;
173 }
174
175 resize(newGlassW: number, newGlassH: number): void {
176 const scaleH = newGlassH / this.glassH;
177 for (let i = 0; i < this.columns; i++) this.h[i] = this.h[i]! * scaleH;
178 this.volume = this.volume * scaleH * (newGlassW / this.glassW);
179 this.glassW = newGlassW;
180 this.glassH = newGlassH;
181 this.u.fill(0);
182 }
183}