Exchange over a level: diffusion, to the unit¶
Diffusion over a level, conserved to the unit.
Fick's first law in the explicit finite-volume form: the flux through each face between an item and its neighbour is r = D dt / dx^2 times the difference of their amounts, truncated to a whole unit toward zero (diffusion_flux, the face's law), and an item's inflow is the sum over its faces (adjacent: two a line, four a plane, six a space; none past a wall), bound board.amount += lusine.laws.diffusion(...). Each face's flux is a whole unit before the sum (a kernel applied to each member), and the face's law is antisymmetric, diffusion_flux(a, b) = -diffusion_flux(b, a), so what one item gains through a face its neighbour loses, and the level's total is kept exactly, whatever the rounding. Truncation never moves more than r times the difference, so at r <= 1 / (2 x axes) no item passes its neighbours (the explicit form's stability). An amount is an integer counted in whole units of its base unit, as a kernel's integer is; D, the cell's extent dx (the lattice's spacing) and dt are the game's. Source: Crank, The Mathematics of Diffusion, 2nd ed., ch. 8; Patankar, Numerical Heat Transfer and Fluid Flow, ch. 4. Checked against a reference in numpy, which has its own tests.
- Source
- Crank, The Mathematics of Diffusion, ch. 8; Patankar, Numerical Heat Transfer and Fluid Flow, ch. 4.
- Proof
- curve, invariant, limit.
diffusion.trait
namespace lusine.laws;
integer<_> diffusion_flux(integer<_> mine, integer<_> theirs, float rate)
{
let flux = rate * (theirs - mine);
return flux < 0.0 * flux ? ceil(flux) : floor(flux);
}
integer<_> diffusion(integer<_> amount, float<m2/s> diffusivity, float<m> extent, float<s> dt)
{
let rate = diffusivity * dt / (extent * extent);
return sum(diffusion_flux(mine = amount, theirs = adjacent.amount, rate));
}
Fick's first law in the explicit finite-volume form, whole units moving between neighbours. The flux through each face is diffusion_flux, r = D dt / dx² times the difference, truncated to a whole unit toward zero; an item's inflow is that kernel applied to each of its neighbours and summed, sum(diffusion_flux(mine = amount, theirs = adjacent.amount, rate)), and the binding adds it, board.amount += .... Each face is whole before the sum and the face's law is antisymmetric, so the level's total never changes by a unit.
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
F32 = np.float32
The kernel¶
The law as the library states it, and as this page runs it: the face's law and the gather over the faces.
print(trait_run.source('laws/diffusion.trait'))
// lusine.laws.diffusion: diffusion over a level, conserved to the unit. Fick's first law in the explicit finite-volume
// form: the flux through each face between an item and its neighbour is r = D dt / dx^2 times the difference of their
// amounts, truncated to a whole unit toward zero (diffusion_flux, the face's law), and an item's inflow is the sum over
// its faces (adjacent: two a line, four a plane, six a space; none past a wall), bound `board.amount +=
// lusine.laws.diffusion(...)`. Each face's flux is a whole unit before the sum (a kernel applied to
// each member), and the face's law is antisymmetric, diffusion_flux(a, b) = -diffusion_flux(b, a), so what one item
// gains through a face its neighbour loses, and the level's total is kept exactly, whatever the rounding. Truncation
// never moves more than r times the difference, so at r <= 1 / (2 x axes) no item passes its neighbours (the explicit
// form's stability). An amount is an integer counted in whole units of its base unit, as a kernel's integer is; D,
// the cell's extent dx (the lattice's spacing) and dt are the game's.
// Source: Crank, The Mathematics of Diffusion, 2nd ed., ch. 8; Patankar, Numerical Heat Transfer and Fluid Flow, ch. 4.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.laws;
integer<_> diffusion_flux(integer<_> mine, integer<_> theirs, float rate)
{
let flux = rate * (theirs - mine);
return flux < 0.0 * flux ? ceil(flux) : floor(flux);
}
integer<_> diffusion(integer<_> amount, float<m2/s> diffusivity, float<m> extent, float<s> dt)
{
let rate = diffusivity * dt / (extent * extent);
return sum(diffusion_flux(mine = amount, theirs = adjacent.amount, rate));
}
The program¶
A line or a plane of cells, each holding a whole number of moles, walled at its edges; the settings hold D, the cell's extent and the step, and main() adds each cell's inflow to it. steps(amounts, n, D, dx, dt) runs n steps from amounts laid out [x] or [x, y], giving the amounts after each. The planted faults below run the same program with the page's own wrong face in the library's place.
PROGRAM = '''namespace level;
lusine.version version = { { number = 1 } };
struct cell { integer<mol> amount; };
lattice<cell, %d> board(spacing = 1 m, boundary = wall);
host struct { float<m2/s> diffusivity; float<m> extent; float<s> dt; } settings;
%s
void main() { board.amount += %s(board.amount, settings.diffusivity, settings.extent, settings.dt); }
'''
LAW = ('', 'lusine.laws.diffusion')
print(PROGRAM % (2, *LAW))
def steps(amounts, n, D, dx, dt, law=LAW):
a = np.asarray(amounts, dtype=np.int64)
run = trait_run.Program(['laws/diffusion.trait'], PROGRAM % (a.ndim, *law), {'level.board': list(a.shape)})
run.store('level.board', 'amount', a)
for name, value in (('diffusivity', D), ('extent', dx), ('dt', dt)):
run.store('level.settings', name, value)
return run.steps('level.main', n, [('level.board', 'amount')])[0].astype(np.int64)
namespace level;
lusine.version version = { { number = 1 } };
struct cell { integer<mol> amount; };
lattice<cell, 2> board(spacing = 1 m, boundary = wall);
host struct { float<m2/s> diffusivity; float<m> extent; float<s> dt; } settings;
void main() { board.amount += lusine.laws.diffusion(board.amount, settings.diffusivity, settings.extent, settings.dt); }
Curve: a cosine decays as exp(−D k² t)¶
A cosine of whole moles on a line of 64 cells, 10⁸ about a mean of 10⁸, at r = 0.2: its amplitude, read off each step's amounts, follows the diffusion equation's e^(−D k² t) within 1% over 600 steps.
n, L, Dk = 64, 1.0, 1e-3
dx = L / n
x = (np.arange(n) + 0.5) * dx
k = 2 * np.pi / L
a0 = (10 ** 8 + np.round(10 ** 7 * np.cos(k * x))).astype(np.int64)
dt = 0.2 * dx * dx / Dk
def decays(a):
t = dt * np.arange(1, len(a) + 1)
amp = (a.astype(float) - 10 ** 8) @ np.cos(k * x) * 2 / n
return bool(np.max(np.abs(amp / (10 ** 7 * np.exp(-Dk * k * k * t)) - 1)) < 0.01), t, amp
a = steps(a0, 600, Dk, dx, dt)
ok, t, amp = decays(a)
assert ok
fig, ax = plt.subplots(figsize=(7, 3.5))
ax.plot(t, 10 ** 7 * np.exp(-Dk * k * k * t), label='exp(-D k² t)')
ax.plot(t[::20], amp[::20], 'o', ms=3, label='exchange, whole moles')
ax.set_xlabel('t (s)')
ax.set_ylabel('amplitude (mol)')
ax.legend()
plt.show()
Reference: bit for bit, in binary32 and with the reference¶
The kernel is a binary32 model bit for bit over those 600 steps: r rounded as the kernel computes it, fl(fl(D dt) / fl(dx dx)), each face's flux fl(r × fl(theirs − mine)) truncated toward zero, and the faces added as integers. Where r is a dyadic number and the differences are under 2²¹, every product is exact in binary32, so on a plane of random amounts at r = 1/8 the kernel is the float64 reference itself , bit for bit, step after step.
sys.path.insert(0, os.path.join(root, 'kernels'))
from lusine.kernels.laws.diffusion import diffusion as reference
def model(a, D, dx, dt):
a = np.asarray(a, dtype=np.int64)
r = F32(F32(F32(D) * F32(dt)) / F32(F32(dx) * F32(dx)))
out = a.copy()
for axis in range(a.ndim):
lo = [slice(None)] * a.ndim
hi = [slice(None)] * a.ndim
lo[axis], hi[axis] = slice(0, -1), slice(1, None)
flux = np.trunc(F32(r * F32(a[tuple(hi)] - a[tuple(lo)]))).astype(np.int64) # into lo, out of hi
out[tuple(lo)] += flux
out[tuple(hi)] -= flux
return out
def follows(a, a0, step):
want, b = [], a0
for _ in range(len(a)):
b = step(b)
want.append(b)
return bool(np.array_equal(a, np.array(want)))
assert follows(a, a0, lambda b: model(b, Dk, dx, dt))
rng = np.random.default_rng(93)
p0 = rng.integers(0, 10 ** 6, (24, 17))
p = steps(p0, 50, 0.125, 1.0, 1.0)
assert follows(p, p0, lambda b: reference(b, 0.125, 1.0, 1.0))
print('600 steps bit for bit with the binary32 model; 50 steps of a plane bit for bit with the float64 reference')
600 steps bit for bit with the binary32 model; 50 steps of a plane bit for bit with the float64 reference
Invariant: the total to the unit, no item past its neighbours, and the mirror¶
On a walled plane of random amounts at the stability limit, r = 1/4, the total after each of 300 steps is the first's to the mole, and at every step every cell stays between the least and the greatest of itself and its neighbours a step before (the discrete maximum principle): truncation never moves more than r times a difference, and at r = 1 / (2 × axes) that is a mean of the neighbourhood. A field mirrored steps to the mirror of its step.
f0 = rng.integers(0, 10, (16, 12))
def conserved(a, a0):
return bool(np.all(a.reshape(len(a), -1).sum(axis=1) == a0.sum()))
def neighbourhood(b, pick):
"""Each cell's pick (np.minimum or np.maximum) of itself and its neighbours, walls excluded."""
out = b.copy()
for axis in range(b.ndim):
for shift in (1, -1):
moved = np.roll(b, shift, axis=axis)
edge = [slice(None)] * b.ndim
edge[axis] = 0 if shift == 1 else -1
moved[tuple(edge)] = b[tuple(edge)]
out = pick(out, moved)
return out
def bounded(a, a0):
"""Every step, every cell between the least and the greatest of itself and its neighbours a step before."""
before = np.concatenate([a0[None], a[:-1]])
return all(np.all(neighbourhood(b, np.minimum) <= c) and np.all(c <= neighbourhood(b, np.maximum))
for b, c in zip(before, a))
f = steps(f0, 300, 0.25, 1.0, 1.0)
assert conserved(f, f0) and bounded(f, f0)
assert np.array_equal(steps(f0[::-1, ::-1], 30, 0.25, 1.0, 1.0), f[:30, ::-1, ::-1])
print(f'the total {int(f0.sum())} mol after every one of 300 steps; every cell within its neighbourhood every step')
the total 850 mol after every one of 300 steps; every cell within its neighbourhood every step
Limit: a flat level stays flat; two cells share¶
A level with no differences moves nothing; two neighbours with a difference each move toward the other, the total kept.
flat = steps(np.full((8, 5), 12345), 3, 0.25, 1.0, 1.0)
assert np.all(flat == 12345)
two = steps([1000, 0], 1, 0.1, 1.0, 1.0)[0]
assert two.tolist() == [900, 100]
print('flat stays flat; two cells after a step:', two)
flat stays flat; two cells after a step: [900 100]
Planted faults¶
Three wrong faces, written here in Trait and run through the same program in the library's place: the flux rounded to the nearest mole rather than truncated (the total kept, but at r = 1/2 a cell overshoots its neighbours), the flux floored (no longer antisymmetric: the total drifts), and the flux the wrong way (toward the fuller cell). The checks above each catch one.
def wrong(face):
return ('''
integer<_> face(integer<_> mine, integer<_> theirs, float rate) { %s }
integer<_> inflow(integer<_> amount, float<m2/s> diffusivity, float<m> extent, float<s> dt)
{
let rate = diffusivity * dt / (extent * extent);
return sum(face(mine = amount, theirs = adjacent.amount, rate));
}
''' % face, 'inflow')
line = rng.integers(0, 10, 40)
rounded = steps(line, 20, 0.5, 1.0, 1.0, wrong('return rate * (theirs - mine);'))
assert conserved(rounded, line) and not bounded(rounded, line)
floored = steps(f0, 50, 0.25, 1.0, 1.0, wrong('return floor(rate * (theirs - mine));'))
assert not conserved(floored, f0)
backward = steps(a0, 600, Dk, dx, dt, wrong('let f = rate * (mine - theirs); return f < 0.0 * f ? ceil(f) : floor(f);'))
assert not decays(backward)[0]
law_line = steps(line, 20, 0.5, 1.0, 1.0)
assert conserved(law_line, line) and bounded(law_line, line)
print('PLANTED: rounding overshoots at r = 1/2, flooring loses moles, the wrong way grows the cosine; each caught')
PLANTED: rounding overshoots at r = 1/2, flooring loses moles, the wrong way grows the cosine; each caught