Erosion of a cohesive bed¶
The mass a flow lifts from a bed of cohesive sediment, per area and time.
Partheniades's law: nothing below the bed's critical shear stress, and above it a rate growing with the excess, E = M (tau / tau_c - 1) where tau > tau_c, else 0, M the bed's erosion rate constant (kg per square metre per second) and tau the flow's shear stress on the bed. Source: Partheniades, J. Hydraul. Div. ASCE 91 (1), 105 (1965); Winterwerp and van Kesteren, Introduction to the Physics of Cohesive Sediment in the Marine Environment (2004), section 10.2. Checked against a reference in numpy, which has its own tests.
- Source
- Partheniades, J. Hydraul. Div. ASCE 91 (1), 105 (1965); Winterwerp and van Kesteren (2004), section 10.2.
- Proof
- reference, invariant, limit.
erosion.trait
namespace lusine.laws;
float<kg/m2/s> erosion(float<Pa> stress, float<Pa> critical, float<kg/m2/s> rate)
{
return stress > critical ? rate * (stress / critical - 1.0) : 0.0 * rate;
}
Partheniades's law, the mass a flow lifts from a bed per area and time, nothing below the bed's critical shear stress and growing linearly past it.
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
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24
rng = np.random.default_rng(2026)
The kernel¶
The law as the library states it.
print(trait_run.source('laws/erosion.trait'))
// lusine.laws.erosion: the mass a flow lifts from a bed of cohesive sediment, per area and time.
// Partheniades's law: nothing below the bed's critical shear stress, and above it a rate growing with the excess,
// E = M (tau / tau_c - 1) where tau > tau_c, else 0,
// M the bed's erosion rate constant (kg per square metre per second) and tau the flow's shear stress on the bed.
// Source: Partheniades, J. Hydraul. Div. ASCE 91 (1), 105 (1965); Winterwerp and van Kesteren, Introduction to the
// Physics of Cohesive Sediment in the Marine Environment (2004), section 10.2.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.laws;
float<kg/m2/s> erosion(float<Pa> stress, float<Pa> critical, float<kg/m2/s> rate)
{
return stress > critical ? rate * (stress / critical - 1.0) : 0.0 * rate;
}
The program¶
Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.
print(trait_run.law_program('laws/erosion.trait', 'erosion'))
from lusine.kernels.laws.erosion import erosion as reference
def law(**kw):
return trait_run.law('laws/erosion.trait', 'erosion', **kw)
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float<Pa> stress;
host float<Pa> critical;
host float<kg/m2/s> rate;
float<kg/m2/s> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.laws.erosion(items.stress, items.critical, items.rate); }
Reference: stresses round the threshold¶
Stresses from a tenth to ten times the critical: the kernel within binary32's roundings of the reference.
tc = rng.uniform(0.05, 1.0, 400)
s = tc * rng.uniform(0.1, 10.0, 400)
got, want = law(stress=s, critical=tc, rate=1e-4), reference(s, tc, 1e-4)
scale = 1e-4 * s / tc
err = np.max(np.abs(got - want) / scale)
assert err <= 8 * U, err / U
print(f'within {err / U:.1f} x 2^-24 over {len(s)} stresses, {int((s <= tc).sum())} of them below the threshold')
within 2.3 x 2^-24 over 400 stresses, 35 of them below the threshold
Invariant: the threshold and the line¶
Nothing is lifted at or below the critical stress, and past it the flux grows by M / tau_c per pascal.
s = np.linspace(0.0, 1.0, 41)
e = law(stress=s, critical=0.2, rate=1e-4)
assert np.all(e[s <= 0.2] == 0.0)
assert np.allclose(np.diff(e[s > 0.2]), 1e-4 / 0.2 * 0.025, rtol=1e-5)
print('zero to the threshold, then 5e-4 kg/(m^2 s) a pascal')
zero to the threshold, then 5e-4 kg/(m^2 s) a pascal
Limit: twice the critical stress¶
At twice the critical stress the flux is the rate constant M.
e = law(stress=np.array([0.4, 2.0]), critical=np.array([0.2, 1.0]), rate=1e-4)
assert np.allclose(e, 1e-4, rtol=4 * U)
print('M at twice the threshold')
M at twice the threshold