Erosion of a cohesive bed, live

Partheniades's law, the mass a flow lifts from a bed of cohesive sediment per area and time: nothing below the bed's critical shear stress τc, and above it a rate growing with the excess, E = M (τ / τc − 1), M the bed's erosion rate constant and τ the flow's shear stress on the bed. Below, a bed with a critical stress of 0.2 Pa and a rate constant of 0.005 kg/(m² s) is put under a shear stress swept from nothing to 1 Pa, 0.0125 Pa a step, twenty steps a second, and at each the library's law (lusine.laws.erosion) gives the mass lifted on your GPU as you read.

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 bed: the chart draws the mass lifted against the shear stress, a point a step over one whole sweep: flat at nothing to the critical stress, then a straight line. Below it, the stress and the mass lifted there.
• The checks: the slope from a step's stress less to the stress now, beside M / τc, the line's; threshold, true while nothing has been lifted at or below the critical stress; a line past it, while that slope has stayed within 0.01% of M / τc at every step past the threshold; M at twice, while the mass lifted at twice the critical stress has been the rate constant itself.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The critical stress, the rate constant and the sweep are the program's settings, which the host may write; this page offers no control to change them yet.

The program

The sweep, the law's answers and the checks, the panels the page draws, and the one procedure, main, run each step: the sweep moves on, the law is asked at the stress, a step's stress less and twice the critical stress, and the checks are taken.

// Partheniades's law's live page: the shear stress a flow puts on a cohesive bed swept from nothing to five times the
// bed's critical stress, the library's erosion at each, with the page's checks computed here on the GPU and shown by
// the panels the program declares.

namespace erosion;

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

// The reader's settings: the bed's critical shear stress and erosion rate constant, and the sweep's step.
host struct {
  float<Pa> critical = 0.2 Pa; float<kg/m2/s> rate = 0.005 kg/m2/s; float<Pa> stride = 0.0125 Pa;
} settings;

// The sweep: the step it is at, from 0 to 80 and round again, and the flow's shear stress on the bed there.
struct { integer step; float<Pa> stress; } sweep;

// The mass lifted at the stress, and the slope from a step's stress less to it, against the line's, M / tau_c.
struct { float<kg/m2/s> lifted; float<kg/(m2 s Pa)> slope; float<kg/(m2 s Pa)> expected; } law;

// The checks: nothing lifted at or below the critical stress; past it, a line of slope M / tau_c, within 0.01% at every
// step; and at twice the critical stress the rate constant itself, M.
struct { bool threshold = true; bool line = true; bool twice = 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);
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the bed", chart = "erosion", quantity = law.lifted, steps = 81, against = sweep.stress,
                    label = "E" } },
  lusine.value { { panel = "the bed", quantity = sweep.stress, label = "shear stress" } },
  lusine.value { { panel = "the bed", quantity = law.lifted, label = "mass lifted" } },
  lusine.value { { panel = "the checks", quantity = law.slope, label = "slope" } },
  lusine.value { { panel = "the checks", quantity = law.expected, label = "M / tau_c" } },
  lusine.value { { panel = "the checks", quantity = checks.threshold, label = "threshold" } },
  lusine.value { { panel = "the checks", quantity = checks.line, label = "a line past it" } },
  lusine.value { { panel = "the checks", quantity = checks.twice, label = "M at twice" } },
};

void main()
{
  sweep.step = sweep.step >= 80 ? 0 : sweep.step + 1;
  sweep.stress = float(sweep.step) * settings.stride;
  law.lifted = lusine.laws.erosion(sweep.stress, settings.critical, settings.rate);
  law.expected = settings.rate / settings.critical;
  law.slope = sweep.stress - settings.stride > settings.critical
              ? (law.lifted - lusine.laws.erosion(stress = sweep.stress - settings.stride, settings.critical,
                                                  settings.rate)) / settings.stride
              : law.slope;
  checks.threshold = checks.threshold && (sweep.stress > settings.critical || law.lifted == 0.0 kg/m2/s);
  checks.line = checks.line
                && (sweep.stress - settings.stride <= settings.critical
                    || abs(law.slope / law.expected - 1.0) <= 0.0001);
  checks.twice = checks.twice
                 && abs(lusine.laws.erosion(stress = 2.0 * settings.critical, settings.critical, settings.rate)
                        / settings.rate - 1.0) <= 0.000001;
  lamps.colour = lamp_colour(holds = checks.threshold && checks.line && checks.twice);
}

The law's other proof (the mass lifted against a reference over stresses round the threshold) is on its proof page.