lib/accy/src/profiling/sph/kernel.zig
daab053ee43316e1809a84551d573ddd1e5bf3d2
1 const std = @import("std");
2 const accy = @import("accy");
3
4 const Allocator = std.mem.Allocator;
5 const Kernel = accy.kernel;
6
7 pub const particle_count: u32 = 26880;
8 pub const domain_min: f32 = -1.0;
9 pub const domain_max: f32 = 1.0;
10 pub const grid_dims: u32 = 64;
11 pub const cell_size: f32 = (domain_max - domain_min) / @as(f32, @floatFromInt(grid_dims));
12 pub const spacing: f32 = 2.0 / 256.0;
13 pub const smoothing: f32 = 2.0 * spacing;
14 pub const particle_mass: f32 = 1.0;
15 pub const stiffness: f32 = 350.0;
16 pub const gravity: f32 = -9.8;
17 pub const viscosity: f32 = 12.0;
18 pub const dt: f32 = 4.0e-5;
19 pub const wall_damping: f32 = 0.5;
20 pub const wall_margin: f32 = spacing * 0.5;
21 pub const min_dist2: f32 = 1.0e-12;
22
23 const tau: f64 = 2.0 * std.math.pi;
24
25 pub fn poly6Scale() f64 {
26 const h: f64 = smoothing;
27 return 4.0 / (0.5 * tau * std.math.pow(f64, h, 8.0));
28 }
29
30 pub fn spikyGradScale() f64 {
31 const h: f64 = smoothing;
32 return 30.0 / (0.5 * tau * std.math.pow(f64, h, 5.0));
33 }
34
35 pub fn viscLaplacianScale() f64 {
36 const h: f64 = smoothing;
37 return 40.0 / (0.5 * tau * std.math.pow(f64, h, 5.0));
38 }
39
40 pub fn poly6(dist2: f64) f64 {
41 const h2 = @as(f64, smoothing) * smoothing;
42 if (dist2 >= h2) return 0.0;
43 const gap = h2 - dist2;
44 return poly6Scale() * gap * gap * gap;
45 }
46
47 pub fn restDensity() f64 {
48 var total: f64 = 0.0;
49 var di: i64 = -2;
50 while (di <= 2) : (di += 1) {
51 var dj: i64 = -2;
52 while (dj <= 2) : (dj += 1) {
53 const fx = @as(f64, @floatFromInt(di)) * spacing;
54 const fy = @as(f64, @floatFromInt(dj)) * spacing;
55 total += @as(f64, particle_mass) * poly6(fx * fx + fy * fy);
56 }
57 }
58 return total;
59 }
60
61 pub fn pressureOf(density: f64) f64 {
62 return @max(@as(f64, stiffness) * (density - restDensity()), 0.0);
63 }
64
65 pub const threads_per_block: u32 = 256;
66 pub const blocks_per_launch: u32 = particle_count / threads_per_block;
67 pub const offsets_stride: u32 = blocks_per_launch;
68
69 pub const block_columns: u32 = 112;
70 pub const block_rows: u32 = 240;
71
72 pub fn seedDamBreak(xs: []f32, ys: []f32, vxs: []f32, vys: []f32) void {
73 std.debug.assert(xs.len == particle_count);
74 const corner_x: f32 = domain_min + 2.0 * spacing;
75 const corner_y: f32 = domain_min + 2.0 * spacing;
76 for (0..particle_count) |index| {
77 const column = index % block_columns;
78 const row = index / block_columns;
79 const jitter = 0.06 * spacing * (hashUnit(index) - 0.5);
80 xs[index] = corner_x + @as(f32, @floatFromInt(column)) * spacing + jitter;
81 ys[index] = corner_y + @as(f32, @floatFromInt(row)) * spacing;
82 vxs[index] = 0.0;
83 vys[index] = 0.0;
84 }
85 }
86
87 fn hashUnit(seed: usize) f32 {
88 var state: u64 = @as(u64, @intCast(seed)) +% 0x9e3779b97f4a7c15;
89 state = (state ^ (state >> 30)) *% 0xbf58476d1ce4e5b9;
90 state = (state ^ (state >> 27)) *% 0x94d049bb133111eb;
91 state ^= state >> 31;
92 return @as(f32, @floatFromInt(state & 0xffffff)) / 16777215.0;
93 }
94
95 const WalkContext = struct {
96 builder: *Kernel.Builder,
97 sorted: Kernel.Value,
98 offsets: Kernel.Value,
99 cx: Kernel.Value,
100 cy: Kernel.Value,
101 };
102
103 const RowBounds = struct {
104 start: Kernel.Value,
105 end: Kernel.Value,
106 };
107
108 fn emitCellAxis(b: *Kernel.Builder, coordinate: Kernel.Value) !Kernel.Value {
109 const origin = try b.constantFloat(.f32, domain_min);
110 const inv_cell = try b.constantFloat(.f32, 1.0 / cell_size);
111 const scaled = try b.mul(try b.sub(coordinate, origin), inv_cell);
112 const raw_index = try b.cast(scaled, .i32);
113 const zero = try b.constantInt(.i32, 0);
114 const last = try b.constantInt(.i32, @as(i64, grid_dims) - 1);
115 return b.min(try b.max(raw_index, zero), last);
116 }
117
118 fn emitRowBounds(ctx: WalkContext, comptime row_offset: i64) !RowBounds {
119 const b = ctx.builder;
120 const zero = try b.constantInt(.i32, 0);
121 const one = try b.constantInt(.i32, 1);
122 const dims = try b.constantInt(.i32, grid_dims);
123 const cells_total = try b.constantInt(.i32, @as(i64, grid_dims) * grid_dims);
124 const count_value = try b.constantInt(.i32, particle_count);
125
126 const col_lo = try b.max(try b.sub(ctx.cx, one), zero);
127 const col_hi = try b.min(try b.add(ctx.cx, one), try b.sub(dims, one));
128
129 const offset_const = try b.constantInt(.i32, row_offset);
130 const ny = try b.add(ctx.cy, offset_const);
131 const row_valid = try b.and_(
132 try b.compare(.ge, ny, zero),
133 try b.compare(.lt, ny, dims),
134 );
135 const row_base = try b.mul(ny, dims);
136 const first_cell = try b.add(row_base, col_lo);
137 const next_cell = try b.add(try b.add(row_base, col_hi), one);
138 const has_next = try b.compare(.lt, next_cell, cells_total);
139
140 const stride = try b.constantInt(.i32, offsets_stride);
141 const first_safe = try b.select(row_valid, first_cell, zero);
142 const start_loaded = try b.load(ctx.offsets, try b.castIndex(try b.mul(first_safe, stride)));
143 const start_value = try b.cast(start_loaded, .i32);
144
145 const next_guard = try b.and_(row_valid, has_next);
146 const next_safe = try b.select(next_guard, next_cell, zero);
147 const end_loaded = try b.load(ctx.offsets, try b.castIndex(try b.mul(next_safe, stride)));
148 const end_value = try b.select(has_next, try b.cast(end_loaded, .i32), count_value);
149
150 const start_position = try b.select(row_valid, start_value, zero);
151 const end_position = try b.select(row_valid, end_value, zero);
152 return .{
153 .start = try b.castIndex(start_position),
154 .end = try b.castIndex(end_position),
155 };
156 }
157
158 pub fn buildDensityGraph(allocator: Allocator) !Kernel.Graph {
159 var builder = try Kernel.Builder.init(allocator, Kernel.Builder.Limits.standard, "accy_sph_density", &.{
160 Kernel.dynamicBuffer(.f32),
161 Kernel.dynamicBuffer(.f32),
162 Kernel.dynamicBuffer(.f32),
163 Kernel.dynamicBuffer(.i32),
164 Kernel.dynamicBuffer(.f32),
165 });
166 errdefer builder.deinit();
167
168 const density = builder.argument(0);
169 const xs = builder.argument(1);
170 const ys = builder.argument(2);
171 const sorted = builder.argument(3);
172 const offsets = builder.argument(4);
173
174 const gid = try builder.globalId(.x);
175 const px = try builder.load(xs, gid);
176 const py = try builder.load(ys, gid);
177 const walk_ctx = WalkContext{
178 .builder = &builder,
179 .sorted = sorted,
180 .offsets = offsets,
181 .cx = try emitCellAxis(&builder, px),
182 .cy = try emitCellAxis(&builder, py),
183 };
184
185 const zero_f = try builder.constantFloat(.f32, 0.0);
186 const radius2 = try builder.constantFloat(.f32, smoothing * smoothing);
187 const step = try builder.constantIndex(1);
188
189 const VisitContext = struct {
190 xs: Kernel.Value,
191 ys: Kernel.Value,
192 sorted: Kernel.Value,
193 px: Kernel.Value,
194 py: Kernel.Value,
195 radius2: Kernel.Value,
196 zero_f: Kernel.Value,
197 };
198 const visit_ctx = VisitContext{
199 .xs = xs,
200 .ys = ys,
201 .sorted = sorted,
202 .px = px,
203 .py = py,
204 .radius2 = radius2,
205 .zero_f = zero_f,
206 };
207
208 var total = zero_f;
209 inline for ([_]i64{ -1, 0, 1 }) |row_offset| {
210 const bounds = try emitRowBounds(walk_ctx, row_offset);
211 total = try builder.fold(bounds.start, bounds.end, step, total, visit_ctx, struct {
212 fn visit(b: *Kernel.Builder, position: Kernel.Value, acc: Kernel.Value, ctx: VisitContext) !Kernel.Value {
213 const candidate = try b.load(ctx.sorted, position);
214 const j = try b.castIndex(candidate);
215 const dx = try b.sub(ctx.px, try b.load(ctx.xs, j));
216 const dy = try b.sub(ctx.py, try b.load(ctx.ys, j));
217 const dist2 = try b.add(try b.mul(dx, dx), try b.mul(dy, dy));
218 const gap = try b.max(try b.sub(ctx.radius2, dist2), ctx.zero_f);
219 const gap3 = try b.mul(gap, try b.mul(gap, gap));
220 return b.add(acc, gap3);
221 }
222 }.visit);
223 }
224
225 const scale = try builder.constantFloat(.f32, @floatCast(poly6Scale() * particle_mass));
226 try builder.store(try builder.mul(total, scale), density, gid);
227 try builder.return_();
228 return builder.finish();
229 }
230
231 pub fn buildForceGraph(allocator: Allocator) !Kernel.Graph {
232 var builder = try Kernel.Builder.init(allocator, Kernel.Builder.Limits.standard, "accy_sph_force", &.{
233 Kernel.dynamicBuffer(.f32),
234 Kernel.dynamicBuffer(.f32),
235 Kernel.dynamicBuffer(.f32),
236 Kernel.dynamicBuffer(.f32),
237 Kernel.dynamicBuffer(.f32),
238 Kernel.dynamicBuffer(.f32),
239 Kernel.dynamicBuffer(.f32),
240 Kernel.dynamicBuffer(.i32),
241 Kernel.dynamicBuffer(.f32),
242 });
243 errdefer builder.deinit();
244
245 const fx_out = builder.argument(0);
246 const fy_out = builder.argument(1);
247 const xs = builder.argument(2);
248 const ys = builder.argument(3);
249 const vxs = builder.argument(4);
250 const vys = builder.argument(5);
251 const density = builder.argument(6);
252 const sorted = builder.argument(7);
253 const offsets = builder.argument(8);
254
255 const gid = try builder.globalId(.x);
256 const px = try builder.load(xs, gid);
257 const py = try builder.load(ys, gid);
258 const vx_i = try builder.load(vxs, gid);
259 const vy_i = try builder.load(vys, gid);
260 const rho_i = try builder.load(density, gid);
261
262 const zero_f = try builder.constantFloat(.f32, 0.0);
263 const half = try builder.constantFloat(.f32, 0.5);
264 const rest = try builder.constantFloat(.f32, @floatCast(restDensity()));
265 const stiffness_value = try builder.constantFloat(.f32, stiffness);
266 const p_i = try builder.max(try builder.mul(stiffness_value, try builder.sub(rho_i, rest)), zero_f);
267
268 const walk_ctx = WalkContext{
269 .builder = &builder,
270 .sorted = sorted,
271 .offsets = offsets,
272 .cx = try emitCellAxis(&builder, px),
273 .cy = try emitCellAxis(&builder, py),
274 };
275
276 const Walk = struct {
277 position: Kernel.Value,
278 fx: Kernel.Value,
279 fy: Kernel.Value,
280 };
281 const VisitContext = struct {
282 xs: Kernel.Value,
283 ys: Kernel.Value,
284 vxs: Kernel.Value,
285 vys: Kernel.Value,
286 density: Kernel.Value,
287 sorted: Kernel.Value,
288 end: Kernel.Value,
289 px: Kernel.Value,
290 py: Kernel.Value,
291 vx_i: Kernel.Value,
292 vy_i: Kernel.Value,
293 p_i: Kernel.Value,
294 rest: Kernel.Value,
295 stiffness: Kernel.Value,
296 zero_f: Kernel.Value,
297 half: Kernel.Value,
298 one_index: Kernel.Value,
299 radius2: Kernel.Value,
300 min_dist2: Kernel.Value,
301 h_value: Kernel.Value,
302 spiky: Kernel.Value,
303 visc: Kernel.Value,
304 };
305
306 var fx_total = zero_f;
307 var fy_total = zero_f;
308 inline for ([_]i64{ -1, 0, 1 }) |row_offset| {
309 const bounds = try emitRowBounds(walk_ctx, row_offset);
310 const visit_ctx = VisitContext{
311 .xs = xs,
312 .ys = ys,
313 .vxs = vxs,
314 .vys = vys,
315 .density = density,
316 .sorted = sorted,
317 .end = bounds.end,
318 .px = px,
319 .py = py,
320 .vx_i = vx_i,
321 .vy_i = vy_i,
322 .p_i = p_i,
323 .rest = rest,
324 .stiffness = stiffness_value,
325 .zero_f = zero_f,
326 .half = half,
327 .one_index = try builder.constantIndex(1),
328 .radius2 = try builder.constantFloat(.f32, smoothing * smoothing),
329 .min_dist2 = try builder.constantFloat(.f32, min_dist2),
330 .h_value = try builder.constantFloat(.f32, smoothing),
331 .spiky = try builder.constantFloat(.f32, @floatCast(spikyGradScale() * particle_mass)),
332 .visc = try builder.constantFloat(.f32, @floatCast(viscLaplacianScale() * particle_mass * viscosity)),
333 };
334 const walked = try builder.whileLoop(Walk{
335 .position = bounds.start,
336 .fx = fx_total,
337 .fy = fy_total,
338 }, visit_ctx, struct {
339 fn keepGoing(b: *Kernel.Builder, walk: Walk, ctx: VisitContext) !Kernel.Value {
340 return b.compare(.lt, walk.position, ctx.end);
341 }
342 }.keepGoing, struct {
343 fn step(b: *Kernel.Builder, walk: Walk, ctx: VisitContext) !Walk {
344 const candidate = try b.load(ctx.sorted, walk.position);
345 const j = try b.castIndex(candidate);
346 const dx = try b.sub(ctx.px, try b.load(ctx.xs, j));
347 const dy = try b.sub(ctx.py, try b.load(ctx.ys, j));
348 const dist2 = try b.add(try b.mul(dx, dx), try b.mul(dy, dy));
349 const within = try b.and_(
350 try b.compare(.lt, dist2, ctx.radius2),
351 try b.compare(.gt, dist2, ctx.min_dist2),
352 );
353 const r = try b.sqrt(try b.max(dist2, ctx.min_dist2));
354 const gap = try b.sub(ctx.h_value, r);
355 const inv_r = try b.div(try b.constantFloat(.f32, 1.0), r);
356 const dir_x = try b.mul(dx, inv_r);
357 const dir_y = try b.mul(dy, inv_r);
358
359 const rho_j = try b.load(ctx.density, j);
360 const inv_rho_j = try b.div(try b.constantFloat(.f32, 1.0), rho_j);
361 const p_j = try b.max(try b.mul(ctx.stiffness, try b.sub(rho_j, ctx.rest)), ctx.zero_f);
362 const press = try b.mul(
363 try b.mul(try b.mul(ctx.spiky, try b.mul(gap, gap)), inv_rho_j),
364 try b.mul(try b.add(ctx.p_i, p_j), ctx.half),
365 );
366 const visc_common = try b.mul(try b.mul(ctx.visc, gap), inv_rho_j);
367 const dfx = try b.add(
368 try b.mul(press, dir_x),
369 try b.mul(visc_common, try b.sub(try b.load(ctx.vxs, j), ctx.vx_i)),
370 );
371 const dfy = try b.add(
372 try b.mul(press, dir_y),
373 try b.mul(visc_common, try b.sub(try b.load(ctx.vys, j), ctx.vy_i)),
374 );
375 return .{
376 .position = try b.add(walk.position, ctx.one_index),
377 .fx = try b.add(walk.fx, try b.select(within, dfx, ctx.zero_f)),
378 .fy = try b.add(walk.fy, try b.select(within, dfy, ctx.zero_f)),
379 };
380 }
381 }.step);
382 fx_total = walked.fx;
383 fy_total = walked.fy;
384 }
385
386 try builder.store(fx_total, fx_out, gid);
387 try builder.store(fy_total, fy_out, gid);
388 try builder.return_();
389 return builder.finish();
390 }
391
392 pub fn buildIntegrateGraph(allocator: Allocator) !Kernel.Graph {
393 var builder = try Kernel.Builder.init(allocator, Kernel.Builder.Limits.standard, "accy_sph_integrate", &.{
394 Kernel.dynamicBuffer(.f32),
395 Kernel.dynamicBuffer(.f32),
396 Kernel.dynamicBuffer(.f32),
397 Kernel.dynamicBuffer(.f32),
398 Kernel.dynamicBuffer(.f32),
399 Kernel.dynamicBuffer(.f32),
400 Kernel.dynamicBuffer(.f32),
401 Kernel.dynamicBuffer(.f32),
402 Kernel.dynamicBuffer(.f32),
403 Kernel.dynamicBuffer(.f32),
404 Kernel.dynamicBuffer(.f32),
405 });
406 errdefer builder.deinit();
407
408 const x_out = builder.argument(0);
409 const y_out = builder.argument(1);
410 const vx_out = builder.argument(2);
411 const vy_out = builder.argument(3);
412 const xs = builder.argument(4);
413 const ys = builder.argument(5);
414 const vxs = builder.argument(6);
415 const vys = builder.argument(7);
416 const fxs = builder.argument(8);
417 const fys = builder.argument(9);
418 const density = builder.argument(10);
419
420 const gid = try builder.globalId(.x);
421 const dt_value = try builder.constantFloat(.f32, dt);
422 const gravity_value = try builder.constantFloat(.f32, gravity);
423 const damping = try builder.constantFloat(.f32, -wall_damping);
424 const wall_lo = try builder.constantFloat(.f32, domain_min + wall_margin);
425 const wall_hi = try builder.constantFloat(.f32, domain_max - wall_margin);
426
427 const inv_rho = try builder.div(try builder.constantFloat(.f32, 1.0), try builder.load(density, gid));
428 const ax = try builder.mul(try builder.load(fxs, gid), inv_rho);
429 const ay = try builder.add(try builder.mul(try builder.load(fys, gid), inv_rho), gravity_value);
430
431 const vx_next = try builder.add(try builder.load(vxs, gid), try builder.mul(ax, dt_value));
432 const vy_next = try builder.add(try builder.load(vys, gid), try builder.mul(ay, dt_value));
433 const x_next = try builder.add(try builder.load(xs, gid), try builder.mul(vx_next, dt_value));
434 const y_next = try builder.add(try builder.load(ys, gid), try builder.mul(vy_next, dt_value));
435
436 const x_inside = try builder.and_(
437 try builder.compare(.ge, x_next, wall_lo),
438 try builder.compare(.le, x_next, wall_hi),
439 );
440 const y_inside = try builder.and_(
441 try builder.compare(.ge, y_next, wall_lo),
442 try builder.compare(.le, y_next, wall_hi),
443 );
444 const vx_final = try builder.select(x_inside, vx_next, try builder.mul(vx_next, damping));
445 const vy_final = try builder.select(y_inside, vy_next, try builder.mul(vy_next, damping));
446 const x_final = try builder.min(try builder.max(x_next, wall_lo), wall_hi);
447 const y_final = try builder.min(try builder.max(y_next, wall_lo), wall_hi);
448
449 try builder.store(x_final, x_out, gid);
450 try builder.store(y_final, y_out, gid);
451 try builder.store(vx_final, vx_out, gid);
452 try builder.store(vy_final, vy_out, gid);
453 try builder.return_();
454 return builder.finish();
455 }
456
457 pub fn referenceDensity(xs: []const f32, ys: []const f32, target: usize) f64 {
458 var total: f64 = 0.0;
459 for (0..xs.len) |other| {
460 const dx = @as(f64, xs[target]) - @as(f64, xs[other]);
461 const dy = @as(f64, ys[target]) - @as(f64, ys[other]);
462 total += @as(f64, particle_mass) * poly6(dx * dx + dy * dy);
463 }
464 return total;
465 }
466
467 pub const Force = struct {
468 fx: f64,
469 fy: f64,
470 magnitude_x: f64,
471 magnitude_y: f64,
472 };
473
474 pub fn referenceForce(
475 xs: []const f32,
476 ys: []const f32,
477 vxs: []const f32,
478 vys: []const f32,
479 densities: []const f32,
480 target: usize,
481 ) Force {
482 const h: f64 = smoothing;
483 const radius2 = h * h;
484 const p_i = pressureOfDevice(densities[target]);
485 var force = Force{ .fx = 0.0, .fy = 0.0, .magnitude_x = 0.0, .magnitude_y = 0.0 };
486 for (0..xs.len) |other| {
487 const dx = @as(f64, xs[target]) - @as(f64, xs[other]);
488 const dy = @as(f64, ys[target]) - @as(f64, ys[other]);
489 const dist2 = dx * dx + dy * dy;
490 if (dist2 >= radius2 or dist2 <= min_dist2) continue;
491 const r = @sqrt(dist2);
492 const gap = h - r;
493 const dir_x = dx / r;
494 const dir_y = dy / r;
495 const rho_j = @as(f64, densities[other]);
496 const p_j = pressureOfDevice(densities[other]);
497 const press = spikyGradScale() * particle_mass * gap * gap / rho_j * (p_i + p_j) * 0.5;
498 const visc_common = viscLaplacianScale() * particle_mass * viscosity * gap / rho_j;
499 const term_x = press * dir_x + visc_common * (@as(f64, vxs[other]) - @as(f64, vxs[target]));
500 const term_y = press * dir_y + visc_common * (@as(f64, vys[other]) - @as(f64, vys[target]));
501 force.fx += term_x;
502 force.fy += term_y;
503 force.magnitude_x += @abs(term_x);
504 force.magnitude_y += @abs(term_y);
505 }
506 return force;
507 }
508
509 fn pressureOfDevice(density: f32) f64 {
510 return @max(@as(f64, stiffness) * (@as(f64, density) - restDensity()), 0.0);
511 }
512
513 pub fn referenceIntegrate(
514 x: f32,
515 y: f32,
516 vx: f32,
517 vy: f32,
518 fx: f32,
519 fy: f32,
520 density: f32,
521 ) [4]f32 {
522 const inv_rho = 1.0 / density;
523 const vx_next = vx + fx * inv_rho * dt;
524 const vy_next = vy + (fy * inv_rho + gravity) * dt;
525 const x_next = x + vx_next * dt;
526 const y_next = y + vy_next * dt;
527 const wall_lo = domain_min + wall_margin;
528 const wall_hi = domain_max - wall_margin;
529 const x_in = x_next >= wall_lo and x_next <= wall_hi;
530 const y_in = y_next >= wall_lo and y_next <= wall_hi;
531 return .{
532 std.math.clamp(x_next, wall_lo, wall_hi),
533 std.math.clamp(y_next, wall_lo, wall_hi),
534 if (x_in) vx_next else vx_next * -wall_damping,
535 if (y_in) vy_next else vy_next * -wall_damping,
536 };
537 }
538
539 test "graphs verify" {
540 inline for (.{ buildDensityGraph, buildForceGraph, buildIntegrateGraph }) |build| {
541 var graph = try build(std.testing.allocator);
542 defer graph.deinit();
543 try graph.verify();
544 }
545 }
546
547 test "rest density matches the seeded lattice interior" {
548 const allocator = std.testing.allocator;
549 const xs = try allocator.alloc(f32, particle_count);
550 defer allocator.free(xs);
551 const ys = try allocator.alloc(f32, particle_count);
552 defer allocator.free(ys);
553 const vxs = try allocator.alloc(f32, particle_count);
554 defer allocator.free(vxs);
555 const vys = try allocator.alloc(f32, particle_count);
556 defer allocator.free(vys);
557 seedDamBreak(xs, ys, vxs, vys);
558
559 const interior = (block_rows / 2) * block_columns + block_columns / 2;
560 const density = referenceDensity(xs, ys, interior);
561 const rest = restDensity();
562 try std.testing.expect(rest > 0.0);
563 try std.testing.expect(@abs(density - rest) / rest < 0.05);
564 }
565
566 test "reference force pushes compressed neighbors apart" {
567 const xs = [_]f32{ 0.0, spacing * 0.5 };
568 const ys = [_]f32{ 0.0, 0.0 };
569 const vxs = [_]f32{ 0.0, 0.0 };
570 const vys = [_]f32{ 0.0, 0.0 };
571 const dense: f32 = @floatCast(restDensity() * 1.5);
572 const densities = [_]f32{ dense, dense };
573 const force = referenceForce(xs[0..], ys[0..], vxs[0..], vys[0..], densities[0..], 0);
574 try std.testing.expect(force.fx < 0.0);
575 try std.testing.expect(force.magnitude_x >= @abs(force.fx));
576 const mirrored = referenceForce(xs[0..], ys[0..], vxs[0..], vys[0..], densities[0..], 1);
577 try std.testing.expect(mirrored.fx > 0.0);
578 }
579
580 test "reference integrate reflects at walls with damping" {
581 const rest: f32 = @floatCast(restDensity());
582 const escaping = referenceIntegrate(domain_max - wall_margin, 0.0, 10.0, 0.0, 0.0, 0.0, rest);
583 try std.testing.expectEqual(domain_max - wall_margin, escaping[0]);
584 try std.testing.expect(escaping[2] < 0.0);
585 const falling = referenceIntegrate(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, rest);
586 try std.testing.expect(falling[3] < 0.0);
587 }