The site’s own reference code for this law, in TypeScript. It computes the numbers above when the page runs it, and it is what the page uses wherever a FrankenSim result is not available.
driftDiffusionFrames1d · src/physics/reference/diffusion/driftDiffusion.ts · this site’s source, as built · sha256:af1cc1a77a820260026aa7f774161b211ece08205f7c21eac6447d1ab98ce6d1
This function computes the listed outputs when it runs.
In words
An explicit finite-volume solve of the drift-diffusion equation on a one-dimensional grid with zero-flux walls. It returns every frame's density, the cell width, the drift velocity, the Peclet number that says which of the two transports dominates, and the stability ratio the step was accepted under. At the defaults of this laboratory the Peclet number is 3.90, so drift is a few times stronger than diffusion, and the stability ratio is 0.0052, far inside its limit of 1.
Implementation
export function driftDiffusionFrames1d(
p: DriftDiffusionParameters,
budget = REFERENCE_FTCS_BUDGET,
): Computation<DriftDiffusionFrames> {
const {
cells,
width,
frames,
stepsPerFrame,
dt,
kickDiffusivity,
mobility,
force,
temperature,
profile,
set,
} = p;
if (!integer(cells, 3) || !integer(frames, 1) || !integer(stepsPerFrame, 1)) {
return invalid(
["cells", "frames", "stepsPerFrame"],
"Use at least three cells and positive whole-number frame and step counts.",
);
}
if (
![width, dt, kickDiffusivity, mobility, force, temperature].every(Number.isFinite) ||
width <= 0 ||
dt <= 0 ||
kickDiffusivity < 0 ||
mobility <= 0 ||
temperature <= 0
) {
return invalid(
["width", "dt", "kickDiffusivity", "mobility", "force", "temperature"],
"Use positive width, dt, mobility, and temperature, nonnegative kick diffusivity, and a finite force.",
);
}
const D = kickDiffusivity;
const u = mobility * force;
const dx = width / cells;
if (!Number.isFinite(u) || !Number.isFinite(1 / dx) || dx <= 0) {
return invalid(
["width", "cells", "mobility", "force"],
"The drift speed and cell scale must be representable.",
);
}
const work = cells * (frames - 1) * stepsPerFrame;
if (work > budget.workUnits || cells * frames * 8 > budget.allocationBytes) {
return {
kind: "outcome",
outcome: {
outcome: "budget-exhausted",
message: executionOutcomeRegistry["budget-exhausted"].message,
retry: executionOutcomeRegistry["budget-exhausted"].retry,
requested: { workUnits: work, allocationBytes: cells * frames * 8 },
allowed: { ...budget },
},
};
}
const speeds = transportSpeeds(D, u, dx);
const rate = (speeds.right + speeds.left) / dx;
const sigma = dt * rate;
if (!Number.isFinite(rate) || !Number.isFinite(sigma)) {
return invalid(
["dt", "width", "kickDiffusivity", "force"],
"The explicit transport rate must be finite.",
);
}
if (sigma > 1) {
const dtMax = 1 / rate;
const driftOnly = D === 0;
return {
kind: "refused",
refusal: makeRefusal(
driftOnly ? "drift-cfl-exceeded" : "drift-diffusion-unstable",
{
parameterIds: driftOnly ? ["dt", "force"] : ["dt"],
capabilityId: "diffusion.driftDiffusionFrames1d",
},
{
details: driftOnly
? { courant: sigma, limit: 1, dtMax }
: { ratio: sigma, limit: 1, dtMax },
rankedRepairs: [
{
label: "Reduce the time step to the explicit scheme's positivity limit.",
action: { parameterId: "dt", value: dtMax },
},
{
label: "Use a coarser spatial grid.",
action: { parameterId: "cells", value: Math.max(3, Math.floor(cells / 2)) },
},
...(driftOnly
? [
{
label: "Use a weaker force.",
action: { parameterId: "force", value: force / 2 },
},
]
: []),
],
},
),
};
}
const field = initialProfile(cells, dx, width, D, u, force, temperature, profile, set);
if (!field) {
return invalid(
["profile"],
"The initial profile must have a finite, positive, representable mass.",
);
}
const values = new Float64Array(cells * frames);
values.set(field);
// Reuse numerical buffers. Face objects are needed only for the accepted final frame.
const flux = new Float64Array(cells + 1);
const next = new Float64Array(cells);
for (let frame = 1; frame < frames; frame++) {
for (let s = 0; s < stepsPerFrame; s++) {
const step = advance(field, flux, next, dx, dt, speeds);
if (step.kind === "refused") return step;
}
values.set(field, frame * cells);
}
const pe = D === 0 ? (u === 0 ? 0 : Math.sign(u) * Number.POSITIVE_INFINITY) : u / (D / dx);
return {
kind: "accepted",
data: Object.freeze({
values,
shape: Object.freeze([frames, cells] as const),
dx,
peclet: pe,
sigma,
driftVelocity: u,
// Diagnostics describe this frame, not the field before its last time step.
faceFlux: Object.freeze(faceFlux(field, speeds, D, u)),
boundary: "zero-flux",
ownerId: "diffusion.driftDiffusionFrames1d",
modelIdentity: D === 0 ? DRIFT_ONLY_NO_KICKS : "drift-diffusion",
}),
};
}