Skip to content

The pressure of a column, live

Hydrostatic and lithostatic pressure, the weight of what lies above: dp/dz = ρ g down a column, so the pressure at a cell's centre is the top's plus g dz times the density of the cells above it and half its own, p = p₀ + g dz (Σρ above + ρ / 2). Below, a pond ten metres deep, five metres of water over five of wet sand of twice its density, is laid out as three columns of 200 cells 5 cm high under a standard atmosphere, and the library's law (lusine.laws.pressure_column) gives the pressure in every cell on your GPU, every step. Beside it the same program computes the closed form, the weight of each layer above a cell's centre, p₀ + g ρ₁ z in the water and p₀ + g (ρ₁ H + ρ₂ (z − H)) in the sand, and checks every cell against it. A probe reads the middle column down, a cell a step, fifty steps a second.

This page runs the program on your GPU through WebGPU, so it needs a browser with WebGPU (a recent Chrome, Edge, Safari or Firefox). Without it the program does not run and nothing is drawn below.

What it shows

  • The strip above the controls is green while every check holds, and turns red for good once one fails.
  • The pond: the chart draws the pressure the probe reads against the cell it is at, counted from the top, and the closed form there, a point a step over one whole pass down the column, the two one over the other: a line that steepens where the water gives way to the sand, at cell 100. Below it, the depth of the probed cell's centre and the two pressures there.
  • The checks, each over every cell and held at every step: every cell, the pressure within four binary32 roundings of the closed form; never falls, the pressure going down every column never falling; no gravity, the same pond under no gravity at the top's pressure all the way down, exactly.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The cells' height, gravity, the top's pressure and the two layers are the program's settings, which the host may write; this page offers no control to change them yet.

The program

The pond, the probe and the checks, the panels the page draws, and the two procedures: start, run once before the first step, lays the two layers down; main, each step, takes the pressure and the closed form in every cell, the checks, and the probe's reading.

// The pressure column's live page: a pond ten metres deep, five of water over five of wet sand, three columns of 200
// cells each, the library's pressure in every cell beside the weight of what lies above it in closed form, read off
// down the middle column a cell a step, with the page's checks computed here on the GPU and shown by the panels the
// program declares.

namespace pressure_column;

lusine.version version = { { number = 1 } };

// The reader's settings: the cells' height, gravity, the pressure at the top (a standard atmosphere's), and the two
// layers' densities, the upper's the first hundred rows'.
host struct {
  float<m> height = 0.05 m; float<m/s2> gravity = 9.80665 m/s2; float<Pa> top = 101325.0 Pa;
  float<kg/m3> upper = 1000.0 kg/m3; float<kg/m3> lower = 2000.0 kg/m3; integer rows = 100;
} settings;

// The pond: x across and y down, row 0 the top. Each cell its density, its pressure, the closed form's there, the
// pressure with no gravity, the pressure where the probe is (nothing elsewhere), and whether it holds the checks.
struct cell {
  float<kg/m3> density; float<Pa> pressure; float<Pa> closed_form; float<Pa> weightless; float<Pa> probed;
  float<Pa> probed_closed; bool matches; bool monotone;
};
lattice<cell, 2> board(spacing = 1 m, boundary = wall);

// The probe: the row it reads, from the top down and round again (as a number, the cells above it), the depth of
// that row's centre, and the pressure there, the library's and the closed form's.
struct { integer row = -1; float cell; float<m> depth; float<Pa> pressure; float<Pa> closed_form; } probe;

// The checks: every cell's pressure within four roundings of the closed form; down every column the pressure never
// falling; and with no gravity, every cell at the top's pressure, exactly.
struct { bool matches = true; bool monotone = true; bool weightless = true; } checks;

// A strip of lamps the page draws, lit green while every check holds and red once one fails.
struct lamp { vector<float, 4> colour; };
lattice<lamp, 2> lamps(spacing = 1 m, boundary = wall);

vector<float, 4> lamp_colour(bool holds)
{
  return holds ? vector<float, 4>(0.15, 0.6, 0.3, 1.0) : vector<float, 4>(0.75, 0.15, 0.15, 1.0);
}

// The layers: the upper density in the first rows, the lower below.
float<kg/m3> layered(float<kg/m3> density, float<kg/m3> upper, float<kg/m3> lower, integer rows)
{
  return count(along(-y)) < rows ? upper : lower;
}

// The weight of what lies above a cell's centre, in closed form: p0 + g (rho_1 z) in the upper layer, and
// p0 + g (rho_1 H + rho_2 (z - H)) in the lower, z the centre's depth and H the upper layer's.
float<Pa> hydrostatic(float<kg/m3> density, float<m> height, float<m/s2> gravity, float<Pa> top, float<kg/m3> upper,
                      float<kg/m3> lower, integer rows)
{
  let z = (float(count(along(-y))) + 0.5) * height;
  let layer = float(rows) * height;
  if (z < layer) return top + gravity * upper * z;
  return top + gravity * (upper * layer + lower * (z - layer));
}

bool close_to(float<Pa> pressure, float<Pa> closed_form)
{
  return abs(pressure - closed_form) <= 0.00000024 * closed_form;
}

bool never_falls(float<Pa> pressure)
{
  if (count(along(-y)) == 0) return true;
  return pressure >= max(along(-y).pressure);
}

// The cell's pressure where it is the middle column's at the probe's row; nothing elsewhere.
float<Pa> at_probe(float<Pa> pressure, integer row)
{
  return count(along(-y)) == row && count(along(-x)) == 1 ? pressure : 0.0 Pa;
}

float<Pa> picked(table<float<Pa>> probed) { return sum(probed); }
bool all_hold(table<bool> matches) { return all(matches); }
bool none_at_top(table<float<Pa>> weightless, float<Pa> top) { return all(weightless == top); }

output table<lusine.widget> panels = {
  lusine.series { { panel = "the pond", chart = "pressure", quantity = probe.pressure, steps = 200,
                    against = probe.cell, label = "the column's" } },
  lusine.series { { panel = "the pond", chart = "pressure", quantity = probe.closed_form, steps = 200,
                    against = probe.cell, label = "closed form" } },
  lusine.value { { panel = "the pond", quantity = probe.depth, label = "depth" } },
  lusine.value { { panel = "the pond", quantity = probe.pressure, label = "pressure" } },
  lusine.value { { panel = "the pond", quantity = probe.closed_form, label = "closed form" } },
  lusine.value { { panel = "the checks", quantity = checks.matches, label = "every cell" } },
  lusine.value { { panel = "the checks", quantity = checks.monotone, label = "never falls" } },
  lusine.value { { panel = "the checks", quantity = checks.weightless, label = "no gravity" } },
};

void start()
{
  board.density = layered(board.density, settings.upper, settings.lower, settings.rows);
}

void main()
{
  board.pressure = lusine.laws.pressure_column(board.density, settings.height, settings.gravity, settings.top);
  board.weightless = lusine.laws.pressure_column(board.density, settings.height, gravity = 0.0 m/s2, settings.top);
  board.closed_form = hydrostatic(board.density, settings.height, settings.gravity, settings.top, settings.upper,
                                  settings.lower, settings.rows);
  board.matches = close_to(board.pressure, board.closed_form);
  board.monotone = never_falls(board.pressure);
  checks.matches = checks.matches && all_hold(board.matches);
  checks.monotone = checks.monotone && all_hold(matches = board.monotone);
  checks.weightless = checks.weightless && none_at_top(board.weightless, settings.top);
  probe.row = probe.row >= 2 * settings.rows - 1 ? 0 : probe.row + 1;
  probe.cell = float(probe.row);
  probe.depth = (float(probe.row) + 0.5) * settings.height;
  board.probed = at_probe(board.pressure, probe.row);
  board.probed_closed = at_probe(pressure = board.closed_form, probe.row);
  probe.pressure = picked(board.probed);
  probe.closed_form = picked(probed = board.probed_closed);
  lamps.colour = lamp_colour(holds = checks.matches && checks.monotone && checks.weightless);
}

The law's other checks (bit for bit against a binary32 model and the float64 reference over rock of random densities, the weight all carried, the columns independent, and planted faults the checks each catch) are on its proof page.