The pressure of a column¶
Hydrostatic and lithostatic pressure, the weight of what lies above.
dp/dz = rho g down the column: 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 = p0 + g dz (sum of rho above + rho / 2). The column runs along the lattice's y axis, y down (row 0 the top, as an image's): the cells above an item are along(-y), from its neighbour to the lattice's edge, nearest first; on a periodic y axis, once round. The cell's height dz is the lattice's spacing along y, the gravity the world's (lusine.sky.surface_gravity, or lusine.conventions.g_n), the top's pressure the atmosphere's or zero in a vacuum, all passed in. A float sum over along is added nearest first, one at a time. Source: Turcotte & Schubert, Geodynamics, 3rd ed., 2.2 (lithostatic pressure); Kundu & Cohen, Fluid Mechanics, 6th ed., 1.8 (hydrostatics). Checked against a reference in numpy, which has its own tests.
- Source
- Turcotte & Schubert, Geodynamics, 3rd ed., 2.2 (lithostatic pressure); Kundu & Cohen, Fluid Mechanics, 6th ed., 1.8 (hydrostatics).
- Proof
- curve, invariant, limit.
pressure_column.trait
namespace lusine.laws;
float<Pa> pressure_column(float<kg/m3> density, float<m> height, float<m/s2> gravity, float<Pa> top)
{
return top + gravity * height * (sum(along(-y).density) + 0.5 * density);
}
Hydrostatic and lithostatic pressure, the weight of what lies above, p = p0 + g dz (the density of the cells above + half the cell's own). The cells above an item are along(-y), from its neighbour to the lattice's edge, nearest first: the column runs down the lattice's y axis, row 0 the top.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.kernels._si import G0
F32 = np.float32
U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24
The kernel¶
The law as the library states it, and as this page runs it.
print(trait_run.source('laws/pressure_column.trait'))
// lusine.laws.pressure_column: hydrostatic and lithostatic pressure, the weight of what lies above. dp/dz = rho g down
// the column: 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 = p0 + g dz (sum of rho above + rho / 2). The column runs along the lattice's y axis, y down (row 0 the
// top, as an image's): the cells above an item are along(-y), from its neighbour to the lattice's edge, nearest
// first; on a periodic y axis, once round. The cell's height dz is the lattice's spacing along y, the
// gravity the world's (lusine.sky.surface_gravity, or lusine.conventions.g_n), the top's pressure the atmosphere's or
// zero in a vacuum, all passed in. A float sum over along is added nearest first, one at a time.
// Source: Turcotte & Schubert, Geodynamics, 3rd ed., 2.2 (lithostatic pressure); Kundu & Cohen, Fluid Mechanics, 6th
// ed., 1.8 (hydrostatics).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.laws;
float<Pa> pressure_column(float<kg/m3> density, float<m> height, float<m/s2> gravity, float<Pa> top)
{
return top + gravity * height * (sum(along(-y).density) + 0.5 * density);
}
The program¶
A plane of cells, x across and y down, each holding its density, the host's; the settings hold the cells' height, the gravity and the top's pressure, and main() binds the kernel over the cells. column(density, ...) runs it once for a density laid out [x, y]. The planted faults below run the same program with the page's own wrong kernel in the library's place.
PROGRAM = '''namespace col;
lusine.version version = { { number = 1 } };
struct cell { host float<kg/m3> density; float<Pa> pressure; };
lattice<cell, 2> board(spacing = 1 m, boundary = wall);
host struct { float<m> height; float<m/s2> gravity; float<Pa> top; } settings;
%s
void main() { board.pressure = %s(board.density, settings.height, settings.gravity, settings.top); }
'''
LAW = PROGRAM % ('', 'lusine.laws.pressure_column')
print(LAW)
def column(density, height=1.0, gravity=G0, top=0.0, program=LAW):
density = np.asarray(density, dtype=float)
run = trait_run.Program(['laws/pressure_column.trait'], program, {'col.board': list(density.shape)})
run.store('col.board', 'density', density)
for name, value in (('height', height), ('gravity', gravity), ('top', top)):
run.store('col.settings', name, value)
return run.steps('col.main', 1, [('col.board', 'pressure')])[0][0]
namespace col;
lusine.version version = { { number = 1 } };
struct cell { host float<kg/m3> density; float<Pa> pressure; };
lattice<cell, 2> board(spacing = 1 m, boundary = wall);
host struct { float<m> height; float<m/s2> gravity; float<Pa> top; } settings;
void main() { board.pressure = lusine.laws.pressure_column(board.density, settings.height, settings.gravity, settings.top); }
Curve: p0 + rho g z¶
Three columns of water 200 cells of 5 cm deep under a standard atmosphere: at every depth the pressure is the top's plus ρ g z, z the depth of the cell's centre, to binary32's rounding of the few operations (the density sums are exact here: whole multiples of 1000).
def water_holds(p, n=200, dz=0.05):
z = (np.arange(n) + 0.5) * dz
want = 101325.0 + 1000.0 * G0 * z
return bool(np.all(np.abs(p - want[None, :]) <= 4 * U * want[None, :]))
p = column(np.full((3, 200), 1000.0), 0.05, G0, 101325.0)
assert water_holds(p)
z = (np.arange(200) + 0.5) * 0.05
print(f'at {z[-1]:.3f} m: {p[1, -1]:.2f} Pa, rho g z + p0 = {101325 + 1000 * G0 * z[-1]:.2f}')
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.plot(p[1] / 1000, -z)
ax.set_xlabel('pressure (kPa)')
ax.set_ylabel('height (m)')
plt.show()
at 9.975 m: 199146.34 Pa, rho g z + p0 = 199146.33
Reference: bit for bit, the column in binary32¶
Rock of random densities, 7 columns of 60 cells: the kernel's pressure is a binary32 model's bit for bit, the sum over the cells above added nearest first, one addition at a time, as Trait sums along an axis, and then g dz times (the sum + half the cell's own), plus the top, each operation rounded. The numpy reference in float64 (over rays.parallel's rays down the column) agrees within the bound a sum of n positive terms keeps whatever its order, (n + 4) × 2⁻²⁴ relative.
sys.path.insert(0, os.path.join(root, 'kernels'))
from lusine.kernels.laws.pressure_column import pressure_column as reference
def model(density, height, gravity, top):
rho = F32(density)
out = np.empty(rho.shape, dtype=F32)
for x in range(rho.shape[0]):
for y in range(rho.shape[1]):
above = F32(0.0)
for j in range(y - 1, -1, -1): # along(-y): the neighbour first, the edge last
above = F32(above + rho[x, j])
out[x, y] = F32(top) + F32(F32(F32(gravity) * F32(height)) * F32(above + F32(F32(0.5) * rho[x, y])))
return out
rng = np.random.default_rng(93)
rock = F32(rng.uniform(1500.0, 3300.0, (7, 60))).astype(float)
dz, top = 2.5, 101325.0
p = column(rock, dz, G0, top)
assert np.array_equal(F32(p), model(rock, dz, G0, top))
want = reference(rock, F32(dz), F32(G0), F32(top), down=(0, 1))
rel = np.max(np.abs(p - want) / want)
assert rel <= (60 + 4) * U, rel / U
print(f'{p.size} cells bit for bit with the binary32 model; the float64 reference within {rel / U:.1f} x 2^-24')
420 cells bit for bit with the binary32 model; the float64 reference within 4.5 x 2^-24
Invariant: down a column it never falls, and the bottom carries everything¶
Every addition is of a density at least zero, so going down a column the pressure never falls, in binary32 too; the bottom cell's pressure less half its own weight is the weight of the whole column, g dz Σρ, within the sum's rounding. Each column is its own: changing one leaves every other bit for bit.
q = column(rock, 1.0, G0)
assert np.all(np.diff(q, axis=1) >= 0)
carried = q[:, -1] + G0 * rock[:, -1] / 2
assert np.all(np.abs(carried - G0 * rock.sum(axis=1)) <= 64 * U * G0 * rock.sum(axis=1))
changed = rock.copy()
changed[3] *= 0.5
r = column(changed, 1.0, G0)
assert np.array_equal(np.delete(r, 3, axis=0), np.delete(q, 3, axis=0)) and not np.array_equal(r[3], q[3])
print('monotone down every column; the weight all carried; the columns independent')
monotone down every column; the weight all carried; the columns independent
Limit: nothing above, or no gravity¶
An empty column, or one under no gravity, is at its top's pressure all the way down, exactly; the top row has nothing above it, so its pressure is the top's plus half its own weight.
assert np.all(column(np.zeros((4, 10)), 1.0, G0, 5.0) == 5.0)
assert np.all(column(np.full((4, 10), 2500.0), 1.0, 0.0, 5.0) == 5.0)
top_row = column(rock, 1.0, G0)[:, 0]
assert np.array_equal(F32(top_row), F32(F32(G0) * F32(F32(0.5) * F32(rock[:, 0]))))
print('the top pressure throughout; the top row half its own weight')
the top pressure throughout; the top row half its own weight
Planted faults¶
The checks above catch a law gone wrong. Two wrong kernels, written here in Trait and run through the same program in the library's place: the column summed from below (the weight of what lies under a cell, along(+y)), and the cell's whole weight counted rather than half. Each fails the water curve and the binary32 model.
BELOW = PROGRAM % ('''
float<Pa> beneath(float<kg/m3> density, float<m> height, float<m/s2> gravity, float<Pa> top)
{
return top + gravity * height * (sum(along(+y).density) + 0.5 * density);
}
''', 'beneath')
WHOLE = PROGRAM % ('''
float<Pa> whole(float<kg/m3> density, float<m> height, float<m/s2> gravity, float<Pa> top)
{
return top + gravity * height * (sum(along(-y).density) + density);
}
''', 'whole')
for name, wrong in (('summed from below', BELOW), ('the whole cell', WHOLE)):
caught = [not water_holds(column(np.full((3, 200), 1000.0), 0.05, G0, 101325.0, wrong)),
not np.array_equal(F32(column(rock, dz, G0, top, wrong)), model(rock, dz, G0, top))]
assert all(caught), name
print(f'PLANTED, {name}: the curve and the model each catch it')
PLANTED, summed from below: the curve and the model each catch it PLANTED, the whole cell: the curve and the model each catch it