Simple cellular-automata-based liquid flow simulation, implemented in SQL

README.md

Simple cellular-automata-based liquid flow simulation, implemented in SQL

Loosely based on https://www.reddit.com/r/cellular_automata/comments/6jhdfw/i_used_1dimensional_cellular_automata_to_make_a/. The SQL code in fluid.sql should work with any contemporary version of PostgreSQL.

1. In fluid.sql, uncomment one of the few initial fluid container scenarios (they start with \COPY input(cells) FROM stdin). A sample container has already been uncommented for you. Feel free to create your own.

2. Run the SQL code, e.g., via

   $ psql -f fluid.sql
   

   to obtain a CSV file fluid.csv which represents the state of the container and its fluid as the latter changes over time.

3. Visualize the liquid flow in your terminal by feeding the CSV file into Python script fluid.py:

   $ python3 fluid.py fluid.csv
   

fluid.py

#! /usr/bin/env python3
# -*- coding: utf8 -*-

# Render the 1D liquid flow simulation generated by fluid-1D.sql
# (output in CSV format: ‹iter›|‹x›|‹ground›|‹water›, sorted by ‹iter›, ‹x›)
#
# Usage: ./fluid.py ‹CSV file›

import sys
import csv
import math
import time
from contextlib import contextmanager

# frames per second
FPS = 30

# a column of height 8 is rendered with h/8×█ (base) and a residual character ∊ ▁▂▃▄▅▆▇ (top)
def column(h):
  return (int(h / 8), [None, '▁','▂','▃','▄','▅','▆','▇'][int(h) % 8])

def stone(cell):
  return cell == '░'

# build ground[iter][x] and water[iter][x] from CSV input
with open(sys.argv[1]) as f:
  lines = [line.rstrip() for line in f]

ground = [ int(row[2])   for row in csv.reader(lines) ]
water  = [ float(row[3]) for row in csv.reader(lines) ]

# width/height of scenario (in characters)
width  = max([ int(row[1]) for row in csv.reader(lines) ])
height = int(max([ gw[0] + gw[1] for gw in zip(ground,water) ]) / 8) + 1

ground = [ ground[i:i+width] for i in range(0, len(ground), width) ]
water  = [ water[i:i+width]  for i in range(0, len(water) , width) ]

# number of iterations
iters  = len(ground)

# ANSI terminal escape codes
esc = chr(27)
cls = esc + "[2J" + esc + "[0;0H"
alt = esc + "[?1049h"
org = esc + "[?1049l"
# liquid/solid colors
liquid = esc + "[34;1m"
solid  = esc + "[0m"

# switch to alternate terminal buffer while rendering
@contextmanager
def alternate():
  try:
    print(alt, end='')
    yield
  finally:
    print(org, end='')

# render the grids for all iterations
# NB: each grids[‹iter›] is a list of column strings (flip that by 90° for screen rendering)
grids = [[]] * iters

for iter, gs in enumerate(ground):
  grids[iter] = [[]] * width
  # render ground
  for x, g in enumerate(gs):
    (base, _) = column(g)
    grids[iter][x] = '░' * base

for iter, ws in enumerate(water):
  # render water on top of ground
  for x, w in enumerate(ws):
    (base, top) = column(w)
    grids[iter][x] += '█' * base
    if top:
      grids[iter][x] += top
    grids[iter][x] = grids[iter][x].ljust(height, ' ')

with alternate():
  # screen rendering
  for iter, grid in enumerate(grids):
    print(cls + 'iteration #' + str(iter))
    for y in range(height, 0, -1):
      for x in range(0, width):
        cell = grids[iter][x][y-1]
        print((solid if stone(cell) else liquid) + cell, end='')
      print('\n', end='')
    print('\n')
    if iter == 0:
      input()
    else:
      time.sleep(1 / FPS)
  input()

fluid.sql

-- Liquid flow simulation
-- (example of a cell automaton in which a cell is influenced by as well
--  as influences its neighbours)
--
-- Loosely based on
-- https://www.reddit.com/r/cellular_automata/comments/6jhdfw/i_used_1dimensional_cellular_automata_to_make_a/
--
-- Visualize using ./fluid.py ‹CSV output file›

-- Initial world
--
-- One cell of the autmaton = one column (ground + water) in the world
--
DROP TABLE IF EXISTS fluid;
CREATE TABLE fluid (
  x      int     NOT NULL,                         -- x coordinate
  ground int     NOT NULL CHECK (ground % 8 = 0),  -- elevation at point x (divisible by 8, because of character-based rendering)
  water  numeric NOT NULL,                         -- water volume at point x
  PRIMARY KEY (x)
);


-- A repertoire of initial worlds
--
DROP TABLE IF EXISTS input;
CREATE TEMPORARY TABLE input (
  y     serial,
  cells text NOT NULL
);



-- \COPY input(cells) FROM stdin
-- ░▄           ▄░
-- ░██         ██░
-- ░██         ██░
-- ░██         ██░
-- ░██         ██░
-- ░░░░░░░░░░░░░░░
-- \.


-- \COPY input(cells) FROM stdin
-- ░                  ███░
-- ░                █████░
-- ░                █████░
-- ░                █████░
-- ░                █████░
-- ░                █████░
-- ░                █████░
-- ░                █████░
-- ░                █████░
-- ░░░░░░░░░░░░░░░░░░░░░░░
-- \.


\COPY input(cells) FROM stdin
░         █████        ░
░         █████        ░
░         █████        ░
░         █████        ░
░         █████        ░
░         █████        ░
░         █████        ░
░         █████        ░
░         █████        ░
░░░░░░░░░░░░░░░░░░░░░░░░
\.


-- \COPY input(cells) FROM stdin
-- ░                      ░
-- ░                      ░
-- ░                      ░
-- ░████████      ████████░
-- ░████████      ████████░
-- ░████████      ████████░
-- ░████████      ████████░
-- ░████████      ████████░
-- ░████████      ████████░
-- ░░░░░░░░░░░░░░░░░░░░░░░░
-- \.


-- \COPY input(cells) FROM stdin
-- ░                        ▄░
-- ░                     ████░
-- ░                   ██████░
-- ░░    ░░░         ████████░
-- ░░░  ░░░░░░      █████████░
-- ░░░░░░░░░░░░     █████████░
-- ░░░░░░░░░░░░░░░░░░░░░░░░░░░
-- \.

TABLE input;

-- Turn input table format into valid initial world
--
INSERT INTO fluid(x, ground, water)
  SELECT col.x,
         COALESCE(SUM(8) FILTER (WHERE col.cell = '░')                                 , 0) AS ground,
         COALESCE(SUM(array_position(array['▁','▂','▃','▄','▅','▆','▇','█'], col.cell)), 0) AS water
  FROM   input AS row,
         LATERAL string_to_table(row.cells, NULL) WITH ORDINALITY AS col(cell,x)
  GROUP BY col.x;

TABLE fluid
ORDER BY x;

-- influence of neighboring cells on position x
DROP TYPE IF EXISTS influence CASCADE;
CREATE TYPE influence AS (
  x        int,
  Δwater   numeric,
  Δkinetic numeric
);

-- Save flow simulation result so that we can copy it to a CSV file below
DROP TABLE IF EXISTS simulation;
CREATE TEMPORARY TABLE simulation AS
--
WITH RECURSIVE
sim(iter,x,ground,water,kinetic) AS (
  SELECT 0 AS iter, f.x, f.ground, f.water, 0.0 AS kinetic
  FROM   fluid AS f
    UNION ALL
  (--  prepare for multiple references
   --  to the recursive table sim
   -- ┌────────────────────────────────┐
   WITH sim(iter,x,ground,water,kinetic) AS (
     TABLE sim
   )
   SELECT s0.iter + 1 AS iter, s0.x, s0.ground,
          --  if cell s0 has not been updated, Δwater ≡ NULL (we then keep our current water volume)
          --           ┌───────────────────────┐
          s0.water   + COALESCE(agg.Δwater  , 0) AS water,   -- water volume at x
          s0.kinetic + COALESCE(agg.Δkinetic, 0) AS kinetic  -- kinetic energy at x (< 0: towards left neighbor, > 0: towards right neighbor)
   FROM   --  first reference to recursive table sim:
          --  iterate over all cells s0 to find the changes that apply to it
          --  ┌───┐
          sim AS s0
          -- aggregate all the influences to be applied to current cell s0 (if there are none these yield NULL)
          --                      ┌─────────────┐             ┌────────────────┐
            LEFT OUTER JOIN
          LATERAL (SELECT infs.x, SUM(infs.Δwater) AS Δwater, SUM(infs.Δkinetic) AS Δkinetic
                          -- encodes the rules of the fluid flow automaton:
                          -- the SELECT yields a single array of entries (x,Δwater,Δkinetic)
                          -- to indicate that the cell at x needs to change
                          -- its water volum by Δwater and its kinetic energy by Δkinetic
                   FROM   (SELECT (-- flow to the left
                                   --           potential energy
                                   --        ┌───────────────────┐
                                   CASE WHEN s1.ground + s1.water - s1.kinetic > LAG(s1.ground, 1) OVER horizontal + LAG(s1.water, 1) OVER horizontal + LAG(s1.kinetic, 1) OVER horizontal
                                   THEN array[ROW(s1.x-1,  LEAST(s1.water, s1.ground + s1.water - s1.kinetic - (LAG(s1.ground, 1) OVER horizontal + LAG(s1.water, 1) OVER horizontal + LAG(s1.kinetic, 1) OVER horizontal)) / 4, 0.0),
                                              ROW(s1.x  , -LEAST(s1.water, s1.ground + s1.water - s1.kinetic - (LAG(s1.ground, 1) OVER horizontal + LAG(s1.water, 1) OVER horizontal + LAG(s1.kinetic, 1) OVER horizontal)) / 4, 0.0),
                                              ROW(s1.x-1, 0.0, -LAG(s1.kinetic, 1) OVER horizontal / 2 - LEAST(s1.water, s1.ground + s1.water - s1.kinetic - (LAG(s1.ground, 1) OVER horizontal + LAG(s1.water, 1) OVER horizontal + LAG(s1.kinetic, 1) OVER horizontal)) / 4)
                                             ] :: influence[]
                                   END
                                   ||
                                   -- flow to the right
                                   CASE WHEN s1.ground + s1.water + s1.kinetic > LEAD(s1.ground, 1) OVER horizontal + LEAD(s1.water, 1) OVER horizontal - LEAD(s1.kinetic, 1) OVER horizontal
                                   THEN array[ROW(s1.x+1,  LEAST(s1.water, s1.ground + s1.water + s1.kinetic - (LEAD(s1.ground, 1) OVER horizontal + LEAD(s1.water, 1) OVER horizontal - LEAD(s1.kinetic, 1) OVER horizontal)) / 4, 0.0),
                                              ROW(s1.x  , -LEAST(s1.water, s1.ground + s1.water + s1.kinetic - (LEAD(s1.ground, 1) OVER horizontal + LEAD(s1.water, 1) OVER horizontal - LEAD(s1.kinetic, 1) OVER horizontal)) / 4, 0.0),
                                              ROW(s1.x+1, 0.0, -LEAD(s1.kinetic, 1) OVER horizontal / 2 + LEAST(s1.water, s1.ground + s1.water + s1.kinetic - (LEAD(s1.ground, 1) OVER horizontal + LEAD(s1.water, 1) OVER horizontal - LEAD(s1.kinetic, 1) OVER horizontal)) / 4)
                                             ] :: influence[]
                                   END
                                  ) AS influence
                           --     second reference to recursive table sim
                           --     ┌───────┐
                           FROM   sim AS s1
                           -- window that allows us to inspect cells in the horizontal neighborhood
                           --     ┌───────────────────────────┐
                           WINDOW horizontal AS (ORDER BY s1.x)
                          ) AS inf(influence),
                           --   turn array into table of (x,Δwater,Δkinetic) influence entries
                           --     ┌───────────────────────────┐
                          LATERAL unnest(inf.influence) AS infs
                   GROUP BY infs.x
                   ) AS agg(x, Δwater, Δkinetic)
                   -- find those influences that relate to current cell s0
                   -- ┌────────────┐
                   ON (s0.x = agg.x)
   WHERE  s0.iter < 300
  ) -- inner WITH (non-recursive, allow multiple references to table sim)
) -- top-level WITH RECURSIVE
SELECT s.iter, s.x, s.ground, s.water
FROM   sim AS s
ORDER BY s.iter, s.x;

TABLE simulation
ORDER BY iter, x
LIMIT 100;

-- Export table simulation in CSV format for rendering in the terminal
-- (see Python program fluid.py)
\COPY simulation TO 'fluid.csv' WITH (FORMAT csv);

Very good. Though, who needs Python to plot things 😅

select i, chart
from 
  (select min(iter) as lo, max(iter) as hi from simulation) as lohi,
  lateral generate_series(lo, hi) as s(i),
  lateral (
    with
     
      -- This part is what you can modify to adapt to your own needs
      --------------------------------------------------------------
     
      -- Your data producing query here 
      source (key, value) as (
        select x, water + 1
        from simulation
        where iter = i
      ),
       
      -- Some configuration items:
      constants as (
        select
         
          -- the height of the y axis
          15 as height, 
     
          -- the width of the x axis, if normalise_x, otherwise, ignored
          20 as width, 
     
          -- the bar characters
          '##' as characters,
     
          -- the characters between bars
          ' ' as separator,
           
          -- the padding of the labels on the y axis
          10 as label_pad, 
           
          -- whether to normalise the data on the x axis by
          -- - filling gaps (if int, bigint, numeric, timestamp, 
          --   timestamptz)
          -- - scaling the x axis to "width"
          false as normalise_x
      ),
       
      -- The rest doesn't need to be touched
      --------------------------------------
       
      -- Pre-calculated dimensions of the source data
      source_dimensions (kmin, kmax, kstep, vmin, vmax) as (
        select
          min(key), max(key), 
          (max(key) - min(key)) / max(width), 
          min(water), max(water)
        from source, constants, simulation
      ),
       
      -- Normalised data, which fills the gaps in case the key data
      -- type can be generated with generate_series (int, bigint, 
      -- numeric, timestamp, timestamptz)
      source_normalised (key, value) as (
        select k, coalesce(sum(source.value), 0)
        from source_dimensions
          cross join constants
          cross join lateral 
            generate_series(kmin, kmax, kstep) as t (k)
          left join source 
            on source.key >= t.k and source.key < t.k + kstep
        group by k
      ),
     
      -- Replace source_normalised by source if you don't like the 
      -- normalised version
      actual_source (i, key, value) as (
        select row_number() over (order by key), key, value 
        from source_normalised, constants
        where normalise_x
        union all
        select row_number() over (order by key), key, value
        from source, constants
        where not normalise_x
      ),
         
      -- Pre-calculated dimensions of the actual data
      actual_dimensions (
        kmin, kmax, kstep, vmin, vmax, width_or_count
      ) as (
        select
          min(key), max(key), 
          (max(key) - min(key)) / max(width), 
          min(vmin), max(vmax), 
          case
            when every(normalise_x) then least(max(width), count(*)::int) 
            else count(*)::int
          end
        from actual_source, constants, source_dimensions
      ),
       
      -- Additional convenience
      dims_and_consts as (
        with
          temp as (
            select *, 
            (length(characters) + length(separator)) 
              * width_or_count as bar_width
          from actual_dimensions, constants
        )
        select *,
          (bar_width - length(kmin::text) - length(kmax::text)) 
            as x_label_pad
        from temp
      ),
       
      -- A cartesian product for all (x, y) data points
      x (x) as (
        select generate_series(1, width_or_count) from dims_and_consts
      ),
      y (y) as (
        select generate_series(1, height) from dims_and_consts
      ),
     
      -- Rendering the ASCII chart
      chart (rn, chart) as (
        select
          y,
          lpad(y * (vmax - vmin) / height || '', label_pad) 
            || ' | '
            || string_agg(
                 case
                   when height * actual_source.value / (vmax - vmin) 
                     >= y then characters 
                   else repeat(' ', length(characters)) 
                 end, separator 
                 order by x
               )
        from
          x left join actual_source on actual_source.i = x, 
          y, dims_and_consts
        group by y, vmin, vmax, height, label_pad
        union all
        select
          0, 
          repeat('-', label_pad) 
            || '-+-'
            || repeat('-', bar_width)
        from dims_and_consts
        union all
        select
          -1, 
          repeat(' ', label_pad) 
            || ' | '
            || case
                 when x_label_pad < 1 then ''
                 else kmin || repeat(' ', x_label_pad) || kmax 
               end
        from dims_and_consts
      )
    select chart
    from chart
    order by rn desc
  ) as c
;

Hilarious, Lukas! That leaves rendering the ground column only. A proper choice of characters ▁▂▃▄▅▆▇█' and this is ready to go SQL-only.

[ The original simulation was part of a university course on SQL. There's a limit to how much SQL a student can stomach at a time. 🙃 ]

That leaves rendering the ground column only

Yeah, I did have a stacked bar chart rendering somewhere, but I couldn't find it, so I figured water only will be good enough for now.

Agreed on the stomaching part 😅