Skip to content

Diffusion to the unit, live

Fick's first law in the explicit finite-volume form, whole units moving between neighbours: the flux through each face between a cell and its neighbour is r = D dt / dx² times the difference of their amounts, truncated to a whole unit toward zero, and a cell's inflow is the sum over its faces. Each face's flux is whole before the sum and the face's law is antisymmetric, so what one cell gains its neighbour loses and the total never changes by a unit. Below, a cosine of whole moles, a million about a mean of ten million, is laid down on a line of 64 cells, walled at its ends, and the library's law (lusine.laws.diffusion, bound board.amount += ...) exchanges it at r = 0.2, sixty steps a second, on your GPU as you read. Beside it the same program computes the decay the diffusion equation gives a cosine, A e^(−D k² t), and checks the one against the other at every step.

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 line, the upper strip: each of the 64 cells coloured by its amount, red at the cosine's trough, white at its crests, fading as the cosine decays. Every 600 steps the cosine is laid down again.
  • The strip under it is green while every check holds, and turns red for good once one fails.
  • The line's panel: the chart draws the cosine's amplitude, read off the cells' amounts, and the equation's, against the time since the cosine was laid down, a point a step, the two one over the other. Below it, the time, the two amplitudes, and the line's total.
  • The checks: the chart draws the amplitude's miss of the equation's at each step as a share of what decays allows, 0.1%, on a log scale, so the check holds while it stays below 1. conserved holds while the line's total has been the same to the mole at every step; bounded while every cell, at every step, has stayed 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).

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The diffusivity (0.001 m²/s), the cell's extent (a 64th of a metre), the step (a fifth of dx² / D), the cosine and the run's length are the program's settings, which the host may write; this page offers no control to change them yet.

The program

The line and its colours, the run, the amplitude and the checks, the panels the page draws, and the one procedure, main, run each step: the cosine laid down where a run starts, the exchange, and the checks after it.

// The diffusion law's live page: a cosine of whole moles on a line of cells, walled at its ends, exchanged between
// neighbours by the library's law, beside the decay the diffusion equation gives a cosine, with the page's checks
// computed here on the GPU and shown by the panels the program declares.

namespace diffusion;

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

// The reader's settings: the diffusivity, the cell's extent (the line is a metre of 64 cells), the step (a fifth of
// the explicit form's dx^2 / D), the cosine's mean and amplitude in moles, and the steps a run lasts.
host struct {
  float<m2/s> diffusivity = 0.001 m2/s; float<m> extent = 0.015625 m; float<s> dt = 0.048828125 s;
  integer<mol> mean = 10000000 mol; integer<mol> amplitude = 1000000 mol; integer steps = 600;
} settings;

// The line: each cell its whole moles, its share of the cosine's amplitude (the amount over the mean times the cosine
// at the cell, over half the line), whether it is within its neighbourhood of a step before, and its colour.
struct cell { integer<mol> amount; float<mol> share; bool held; vector<float, 4> colour; };
lattice<cell, 2> board(spacing = 1 m, boundary = wall);

// The run: its steps and its time since the cosine was laid down, and the line's total when it was.
struct { integer step; float<s> time; integer<mol> start_total; } run;

// The cosine's amplitude read off the line, the diffusion equation's for it, A exp(-D k^2 t) with k = 2 pi / L, and
// the line's total; and the checks: the amplitude within 0.1% of the equation's at every step (its miss also as a
// share of that, which the page charts), the total the same to the mole at every step, and every cell at every step
// between the least and the greatest of itself and its neighbours a step before (the discrete maximum principle).
struct {
  float<mol> amplitude; float<mol> closed_form; integer<mol> total; float error_share;
  bool decays = true; bool conserved = true; bool bounded = 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 cosine laid down on the line, once round it: the mean plus the amplitude times cos(k x) at each cell's centre,
// in whole moles; the amount as it was where the run goes on.
integer<mol> laid(integer<mol> amount, integer<mol> mean, integer<mol> amplitude, bool again)
{
  let x = float(count(along(-x))) + 0.5;
  let n = float(count(along(-x)) + count(along(+x)) + 1);
  if (again) return mean + floor(float(amplitude) * cos(6.2831853 * x / n) + 0.5 mol);
  return amount;
}

// The cell's part of the cosine's amplitude: its amount over the mean, times the cosine at the cell, over half the
// line's cells.
float<mol> share(integer<mol> amount, integer<mol> mean)
{
  let x = float(count(along(-x))) + 0.5;
  let n = float(count(along(-x)) + count(along(+x)) + 1);
  return float(amount - mean) * cos(6.2831853 * x / n) * 2.0 / n;
}

bool holds_within(integer<mol> amount)
{
  let before = previous(amount);
  let least = min(previous(adjacent.amount));
  let greatest = max(previous(adjacent.amount));
  return amount >= (before < least ? before : least) && amount <= (before > greatest ? before : greatest);
}

float<mol> amplitude(table<float<mol>> share) { return sum(share); }
integer<mol> total(table<integer<mol>> amount) { return sum(amount); }
bool all_held(table<bool> held) { return all(held); }

float<mol> decayed(integer<mol> amplitude, float<m2/s> diffusivity, float<m> extent, float<s> time)
{
  let k = 6.2831853 / (64.0 * extent);
  return float(amplitude) * exp(-diffusivity * k * k * time);
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the line", chart = "the cosine's amplitude", quantity = checks.amplitude, steps = 600,
                    against = run.time, label = "exchanged" } },
  lusine.series { { panel = "the line", chart = "the cosine's amplitude", quantity = checks.closed_form, steps = 600,
                    against = run.time, label = "A exp(-D k^2 t)" } },
  lusine.value { { panel = "the line", quantity = run.time, label = "time" } },
  lusine.value { { panel = "the line", quantity = checks.amplitude, label = "amplitude" } },
  lusine.value { { panel = "the line", quantity = checks.closed_form, label = "the equation's" } },
  lusine.value { { panel = "the line", quantity = checks.total, label = "total" } },
  lusine.series { { panel = "the checks", chart = "share of the tolerance", quantity = checks.error_share, steps = 600,
                    log = true, label = "amplitude error" } },
  lusine.value { { panel = "the checks", quantity = checks.decays, label = "decays" } },
  lusine.value { { panel = "the checks", quantity = checks.conserved, label = "conserved" } },
  lusine.value { { panel = "the checks", quantity = checks.bounded, label = "bounded" } },
};

void main()
{
  // Every run's steps the cosine is laid down again, and the run's clock and total start again with it.
  board.amount = laid(board.amount, settings.mean, settings.amplitude, again = run.step == 0);
  run.start_total = run.step == 0 ? total(board.amount) : run.start_total;
  run.time = run.step == 0 ? 0.0 s : run.time;
  board.amount += lusine.laws.diffusion(board.amount, settings.diffusivity, settings.extent, settings.dt);
  board.held = holds_within(board.amount);
  checks.bounded = checks.bounded && (run.step == 0 || all_held(board.held));
  run.time = run.time + settings.dt;
  board.share = share(board.amount, settings.mean);
  checks.amplitude = amplitude(board.share);
  checks.closed_form = decayed(settings.amplitude, settings.diffusivity, settings.extent, run.time);
  checks.error_share = abs(checks.amplitude / checks.closed_form - 1.0) / 0.001;
  checks.decays = checks.decays && checks.error_share <= 1.0;
  checks.total = total(board.amount);
  checks.conserved = checks.conserved && checks.total == run.start_total;
  run.step = run.step + 1 >= settings.steps ? 0 : run.step + 1;
  board.colour = lusine.display.ramp(value = float(board.amount), low = float(settings.mean - settings.amplitude),
                                     high = float(settings.mean + settings.amplitude));
  lamps.colour = lamp_colour(holds = checks.decays && checks.conserved && checks.bounded);
}

The law's other checks (bit for bit against a binary32 model and the float64 reference, a plane at the stability limit, a mirrored field, and planted faults the checks each catch) are on its proof page.