Example
sph_double_dam_column_2d
Two-dam SPH setup with a central obstacle column.
// ============================================================================
// 2D SPH Double-Dam Simulation With Central Column
//
// Layout:
// - wide rectangular tank
// - left fluid column at the bottom-left
// - right fluid column at the bottom-right
// - rigid rectangular obstacle column near the center-bottom
//
// Uses the same compact cell-bin traversal as the optimized dam-break example.
// The obstacle is handled with simple reflective collision, not SPH wall
// particles, so this stays consistent with the existing reflective-wall model.
// ============================================================================
function clamp(real64 x, real64 lo, real64 hi) -> real64 {
if (x < lo) {
return lo
}
if (x > hi) {
return hi
}
return x
}
program main {
use array: full_int64
string out_dir = "sph_double_dam_column_2d_output"
ensure_dir(out_dir)
// Wide rectangular tank.
int64 n = 1800
real64 h = 5.0
real64 edge = 5.0
real64 width = 700.0
real64 height = 280.0
real64 m = 1.0
real64 rho0 = 1.0
real64 dt = 0.05
real64 k = 20.0
real64 g = 0.2
real64 visc = 0.5
real64 eps = 0.1
real64 bdry_reflec = 0.5
real64 bdry_proj_dist = 0.0
int64 max_steps = 12000
int64 output_every = 50
int64 progress_every = 100
int64 quiet_steps = 0
int64 quiet_steps_needed = 200
real64 steady_ke_per_particle_tol = 0.0002
real64 x_min = edge
real64 y_min = edge
real64 x_max = width - edge
real64 y_max = height - edge
// Two bottom corner dams.
int64 left_nx = 18
int64 left_ny = 50
int64 right_nx = 18
int64 right_ny = 50
int64 left_count = left_nx * left_ny
// Central rigid obstacle column.
real64 col_width = 10.0 * h
real64 col_height = 22.0 * h
real64 col_x0 = 0.5 * (x_min + x_max) - 0.5 * col_width
real64 col_x1 = col_x0 + col_width
real64 col_y0 = y_min
real64 col_y1 = y_min + col_height
int64 grid_nx = int64((x_max - x_min) / h) + 1
int64 grid_ny = int64((y_max - y_min) / h) + 1
int64 num_cells = grid_nx * grid_ny
// Fluid state.
real64 x[n], y[n], u[n], v[n]
real64 accel_x[n], accel_y[n], rho[n], p[n]
// Compact cell bins and reorder scratch.
int64 particle_cell[n], particle_order[n]
int64 cell_counts[num_cells], cell_starts[num_cells], cell_offsets[num_cells]
real64 x_tmp[n], y_tmp[n], u_tmp[n], v_tmp[n]
real64 accel_x_tmp[n], accel_y_tmp[n], rho_tmp[n], p_tmp[n]
int64 step = 0, i = 0, ix = 0, iy = 0
int64 cell = 0, cx = 0, cy = 0, ncx = 0, ncy = 0
int64 sx0 = 0, sx1 = 0, sy0 = 0, sy1 = 0
int64 start = 0, count = 0, pos = 0, j = 0, running = 0
int64 right_i = 0
real64 dx = 0.0, dy = 0.0, rij = 0.0, r2 = 0.0, t = 0.0
real64 grad_fac = 0.0, lapv = 0.0
real64 grad_px = 0.0, grad_py = 0.0
real64 visc_x = 0.0, visc_y = 0.0
real64 rho_i = 0.0, avg_rho = 0.0, max_speed = 0.0, speed = 0.0
real64 kinetic_energy = 0.0, ke_per_particle = 0.0
real64 bx0 = 0.0, bx1 = 0.0, by0 = 0.0, by1 = 0.0
real64 jitter_x = 0.0, jitter_y = 0.0
real64 h2 = 0.0, poly6_coeff = 0.0, spiky_grad_coeff = 0.0, visc_lap_coeff = 0.0
real64 left_pen = 0.0, right_pen = 0.0, top_pen = 0.0
string csv = "", frame_path = ""
h2 = h * h
poly6_coeff = 315.0 / (64.0 * 3.141592653589793 * h * h * h * h * h * h * h * h * h)
spiky_grad_coeff = -45.0 / (3.141592653589793 * h * h * h * h * h * h)
visc_lap_coeff = 45.0 / (3.141592653589793 * h * h * h * h * h * h)
// Left dam.
do iy = 0, left_ny - 1 {
do ix = 0, left_nx - 1 {
i = ix + left_nx * iy
jitter_x = rand()
jitter_y = rand()
x[i] = x_min + h * real64(ix) + h * eps * jitter_x
y[i] = y_min + h * real64(iy) + h * eps * jitter_y
u[i] = 0.0
v[i] = 0.0
accel_x[i] = 0.0
accel_y[i] = 0.0
rho[i] = rho0
p[i] = 0.0
}
}
// Right dam.
do iy = 0, right_ny - 1 {
do ix = 0, right_nx - 1 {
right_i = left_count + ix + right_nx * iy
if (right_i < n) {
jitter_x = rand()
jitter_y = rand()
x[right_i] = x_max - h * real64(right_nx - ix - 1) - h * eps * jitter_x
y[right_i] = y_min + h * real64(iy) + h * eps * jitter_y
u[right_i] = 0.0
v[right_i] = 0.0
accel_x[right_i] = 0.0
accel_y[right_i] = 0.0
rho[right_i] = rho0
p[right_i] = 0.0
}
}
}
string meta = "max_steps=" + to_string(max_steps) + "\n"
meta += "output_every=" + to_string(output_every) + "\n"
meta += "n=" + to_string(n) + "\n"
meta += "h=" + to_string(h) + "\n"
meta += "dt=" + to_string(dt) + "\n"
meta += "steady_ke_per_particle_tol=" + to_string(steady_ke_per_particle_tol) + "\n"
meta += "quiet_steps_needed=" + to_string(quiet_steps_needed) + "\n"
meta += "x_min=" + to_string(x_min) + "\n"
meta += "x_max=" + to_string(x_max) + "\n"
meta += "y_min=" + to_string(y_min) + "\n"
meta += "y_max=" + to_string(y_max) + "\n"
meta += "col_x0=" + to_string(col_x0) + "\n"
meta += "col_x1=" + to_string(col_x1) + "\n"
meta += "col_y0=" + to_string(col_y0) + "\n"
meta += "col_y1=" + to_string(col_y1) + "\n"
write_file(out_dir + "/meta.txt", meta)
print("SPH 2D double dam: n=", n, "h=", h, "dt=", dt, "steady_ke=", steady_ke_per_particle_tol)
step = 0
while (step <= max_steps && quiet_steps < quiet_steps_needed) {
cell_counts = full_int64(num_cells, 0)
do i = 0, n - 1 {
cx = int64(clamp((x[i] - x_min) / h, 0.0, real64(grid_nx - 1)))
cy = int64(clamp((y[i] - y_min) / h, 0.0, real64(grid_ny - 1)))
cell = cx + cy * grid_nx
particle_cell[i] = cell
cell_counts[cell] += 1
}
running = 0
do cell = 0, num_cells - 1 {
cell_starts[cell] = running
cell_offsets[cell] = running
running += cell_counts[cell]
}
do i = 0, n - 1 {
cell = particle_cell[i]
particle_order[cell_offsets[cell]] = i
cell_offsets[cell] += 1
}
// Reorder by cell to make neighbor spans contiguous.
do pos = 0, n - 1 {
i = particle_order[pos]
x_tmp[pos] = x[i]
y_tmp[pos] = y[i]
u_tmp[pos] = u[i]
v_tmp[pos] = v[i]
accel_x_tmp[pos] = accel_x[i]
accel_y_tmp[pos] = accel_y[i]
rho_tmp[pos] = rho[i]
p_tmp[pos] = p[i]
}
x = x_tmp
y = y_tmp
u = u_tmp
v = v_tmp
accel_x = accel_x_tmp
accel_y = accel_y_tmp
rho = rho_tmp
p = p_tmp
// Density
do i = 0, n - 1 {
rho_i = 0.0
cx = int64(clamp((x[i] - x_min) / h, 0.0, real64(grid_nx - 1)))
cy = int64(clamp((y[i] - y_min) / h, 0.0, real64(grid_ny - 1)))
sx0 = max(0, cx - 1)
sx1 = min(grid_nx - 1, cx + 1)
sy0 = max(0, cy - 1)
sy1 = min(grid_ny - 1, cy + 1)
do ncy = sy0, sy1 {
do ncx = sx0, sx1 {
cell = ncx + ncy * grid_nx
start = cell_starts[cell]
count = cell_counts[cell]
do pos = 0, count - 1 {
j = start + pos
dx = x[i] - x[j]
dy = y[i] - y[j]
r2 = dx * dx + dy * dy
if (r2 <= h2) {
t = h2 - r2
rho_i += m * poly6_coeff * t * t * t
}
}
}
}
rho[i] = rho_i
}
do i = 0, n - 1 {
p[i] = k * (rho[i] - rho0)
}
// Acceleration
do i = 0, n - 1 {
accel_x[i] = 0.0
accel_y[i] = -g
cx = int64(clamp((x[i] - x_min) / h, 0.0, real64(grid_nx - 1)))
cy = int64(clamp((y[i] - y_min) / h, 0.0, real64(grid_ny - 1)))
sx0 = max(0, cx - 1)
sx1 = min(grid_nx - 1, cx + 1)
sy0 = max(0, cy - 1)
sy1 = min(grid_ny - 1, cy + 1)
do ncy = sy0, sy1 {
do ncx = sx0, sx1 {
cell = ncx + ncy * grid_nx
start = cell_starts[cell]
count = cell_counts[cell]
do pos = 0, count - 1 {
j = start + pos
if (j != i) {
dx = x[i] - x[j]
dy = y[i] - y[j]
r2 = dx * dx + dy * dy
if (r2 > 0.0 && r2 <= h2) {
rij = sqrt(r2)
grad_fac = spiky_grad_coeff * (h - rij) * (h - rij) / rij
grad_px = m * (p[i] + p[j]) * grad_fac * dx / (2.0 * rho[i] * rho[j])
grad_py = m * (p[i] + p[j]) * grad_fac * dy / (2.0 * rho[i] * rho[j])
lapv = visc_lap_coeff * (h - rij)
visc_x = m * visc * (u[j] - u[i]) * lapv / (rho[i] * rho[j])
visc_y = m * visc * (v[j] - v[i]) * lapv / (rho[i] * rho[j])
accel_x[i] += grad_px + visc_x
accel_y[i] += grad_py + visc_y
}
}
}
}
}
}
// Explicit Euler
do i = 0, n - 1 {
u[i] += dt * accel_x[i]
v[i] += dt * accel_y[i]
x[i] += dt * u[i]
y[i] += dt * v[i]
}
// Tank walls.
do i = 0, n - 1 {
if (x[i] < x_min) {
u[i] *= -bdry_reflec
x[i] = x_min + bdry_proj_dist
} else if (x[i] > x_max) {
u[i] *= -bdry_reflec
x[i] = x_max - bdry_proj_dist
}
if (y[i] < y_min) {
v[i] *= -bdry_reflec
y[i] = y_min + bdry_proj_dist
} else if (y[i] > y_max) {
v[i] *= -bdry_reflec
y[i] = y_max - bdry_proj_dist
}
}
// Central obstacle column.
do i = 0, n - 1 {
if (x[i] >= col_x0 && x[i] <= col_x1 && y[i] >= col_y0 && y[i] <= col_y1) {
left_pen = abs(x[i] - col_x0)
right_pen = abs(col_x1 - x[i])
top_pen = abs(col_y1 - y[i])
if (top_pen <= left_pen && top_pen <= right_pen) {
v[i] *= -bdry_reflec
y[i] = col_y1 + bdry_proj_dist
} else if (left_pen <= right_pen) {
u[i] *= -bdry_reflec
x[i] = col_x0 - bdry_proj_dist
} else {
u[i] *= -bdry_reflec
x[i] = col_x1 + bdry_proj_dist
}
}
}
max_speed = 0.0
kinetic_energy = 0.0
do i = 0, n - 1 {
speed = sqrt(u[i] * u[i] + v[i] * v[i])
max_speed = max(max_speed, speed)
kinetic_energy += 0.5 * m * (u[i] * u[i] + v[i] * v[i])
}
ke_per_particle = kinetic_energy / real64(n)
if (ke_per_particle < steady_ke_per_particle_tol) {
quiet_steps += 1
} else {
quiet_steps = 0
}
if (step % progress_every == 0) {
avg_rho = 0.0
do i = 0, n - 1 {
avg_rho += rho[i]
}
avg_rho = avg_rho / real64(n)
print("step", step, "/", max_steps, "avg_rho", avg_rho, "max_speed", max_speed, "ke_per_particle", ke_per_particle, "quiet", quiet_steps)
}
if (step % output_every == 0) {
csv = "x,y,vx,vy,rho,p\n"
do i = 0, n - 1 {
csv += to_string(x[i]) + "," + to_string(y[i]) + ","
csv += to_string(u[i]) + "," + to_string(v[i]) + ","
csv += to_string(rho[i]) + "," + to_string(p[i]) + "\n"
}
frame_path = out_dir + "/frame_" + to_string(step) + ".csv"
write_file(frame_path, csv)
}
step += 1
}
avg_rho = 0.0
max_speed = 0.0
bx0 = x[0]
bx1 = x[0]
by0 = y[0]
by1 = y[0]
do i = 0, n - 1 {
avg_rho += rho[i]
speed = sqrt(u[i] * u[i] + v[i] * v[i])
max_speed = max(max_speed, speed)
bx0 = min(bx0, x[i])
bx1 = max(bx1, x[i])
by0 = min(by0, y[i])
by1 = max(by1, y[i])
}
avg_rho = avg_rho / real64(n)
string summary = "particles=" + to_string(n) + "\n"
summary += "executed_steps=" + to_string(step) + "\n"
summary += "quiet_steps=" + to_string(quiet_steps) + "\n"
summary += "avg_rho=" + to_string(avg_rho) + "\n"
summary += "ke_per_particle=" + to_string(ke_per_particle) + "\n"
summary += "max_speed=" + to_string(max_speed) + "\n"
summary += "bounds=" + to_string(bx0) + "," + to_string(bx1) + "," + to_string(by0) + "," + to_string(by1) + "\n"
write_file(out_dir + "/summary.txt", summary)
print("Done:", n, "particles, avg_rho=", avg_rho, "bounds=", bx0, bx1, by0, by1)
}