Standard gravity
64-bitg, used for the liquid head across a valve.
g = 9.80665 m/s²
The exact conventional value, not 9.81.
const G = 9.80665;Resources Numerical accuracy Calculation policy
Dynamic simulation and controller tuning
Simulates tanks, reactors, gas vessels, valves and PID controllers over time and suggests controller settings from a step test.
How the precision is kept
States are 64-bit. The trends are as close as the stiff integrator's step control: by default a relative tolerance of 1e-6 and an absolute tolerance of 1e-9 per step, with the accepted value extrapolated from one full and two half implicit Euler steps. Standard gravity is the exact 9.80665. The tuning rules are published rules of thumb with typed constants.
Calculations under this policy
12 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
process/dynamics.js CadCraft 0.8.2
Tests tools/process-test/test-link.mjs, tools/process-test/test-draw4.mjs
Jacobian by differences
64-bitThe matrix of derivatives of the plant's rates with respect to its states.
Jᵢⱼ = (fᵢ(y + h eⱼ) − fᵢ(y)) ÷ h, h = 1e-7 × max(|yⱼ|, 1e-3)
A one-sided difference, good to about half the digits. It only steers the Newton iteration inside a step; the step's answer is tested against the equations themselves.
// ---- stiff solver ------------------------------------------------------------------------------------------------------------------
function jac(f, t, y, f0) {
const n = y.length, J = Array.from({ length: n }, () => new Array(n).fill(0));
for (let j = 0; j < n; j++) { const h = 1e-7 * Math.max(Math.abs(y[j]), 1e-3), yp = y.slice(); yp[j] += h; const fp = f(t, yp); for (let i = 0; i < n; i++) J[i][j] = (fp[i] - f0[i]) / h; }
return J;
}One implicit Euler step
64-bitThe state after one step of size h.
x − y − h f(t + h, x) = 0, (I − hJ) Δx = −residual, repeated
Newton with the Jacobian held fixed, at most 12 passes, until the largest correction is below 1e-11 of the largest state plus 1e-14. A step that does not get there is rejected by the caller.
function ieStep(f, t, y, h, J) { // one implicit Euler step with the (frozen) Jacobian J; null when Newton fails
const n = y.length, M = Array.from({ length: n }, (_, i) => Array.from({ length: n }, (_, j) => (i === j ? 1 : 0) - h * J[i][j]));
let x = y.slice();
for (let it = 0; it < 12; it++) {
const fx = f(t + h, x), r = x.map((v, i) => v - y[i] - h * fx[i]);
const dx = solveLinear(M.map((row) => row.slice()), r.map((v) => -v)); if (!dx || dx.some((v) => !Number.isFinite(v))) return null;
x = x.map((v, i) => v + dx[i]);
const sc = Math.max(1e-12, ...x.map(Math.abs));
if (Math.max(...dx.map(Math.abs)) < 1e-11 * sc + 1e-14) return x;
}
return null;
}Stiff integrator
64-bitThe state at the end time, with every accepted step on the way.
y_new = 2 y_(two half steps) − y_(one full step); error = max |half − full| ÷ (atol + rtol max(|y|, |y_half|)); h ← h × clamp(0.9 ÷ √error, 0.2, 4)
The error test is on the difference between the two first-order results, and the value kept is their extrapolation. Running out of steps or underflowing the step raises an error, so a result that is returned did reach the end time.
/** y' = f(t, y) from t0 to t1. → { y, t, steps, rejected, points: [[t, y]…] } (every accepted step). `onStep(t, y)` may return a replacement y (events). */
export function stiff(f, t0, y0, t1, { rtol = 1e-6, atol = 1e-9, h0 = null, hmax = Infinity, maxSteps = 500000 } = {}) {
let t = t0, y = y0.slice(), h = Math.min(h0 || (t1 - t0) / 200, hmax), steps = 0, rejected = 0; const points = [[t, y.slice()]];
while (t < t1 - 1e-12 * Math.max(1, Math.abs(t1)) && steps < maxSteps) {
if (t + h > t1) h = t1 - t;
const f0 = f(t, y), J = jac(f, t, y, f0);
const full = ieStep(f, t, y, h, J), half1 = full && ieStep(f, t, y, h / 2, J), half2 = half1 && ieStep(f, t + h / 2, half1, h / 2, J);
if (!half2) { h /= 4; rejected++; steps++; if (h < 1e-15 * Math.max(1, Math.abs(t))) throw new Error('stiff solver step underflow at t = ' + t); continue; }
let err = 0; for (let i = 0; i < y.length; i++) { const sc = atol + rtol * Math.max(Math.abs(y[i]), Math.abs(half2[i])); err = Math.max(err, Math.abs(half2[i] - full[i]) / sc); }
if (err <= 1) { y = half2.map((v, i) => 2 * v - full[i]); t += h; points.push([t, y.slice()]); }
else rejected++;
h = Math.min(hmax, h * Math.min(4, Math.max(0.2, 0.9 / Math.sqrt(err || 1e-12)))); steps++;
}
if (t < t1 - 1e-9 * Math.max(1, Math.abs(t1))) throw new Error('stiff solver ran out of steps before t = ' + t1);
return { y, t, steps, rejected, points };
}Valve characteristic
64-bitThe fraction of full flow at a valve opening.
linear: x; equal percentage: R^(x − 1), 0 at x = 0, R = 50 unless given; quick opening: √x
The opening is clamped to [0, 1] first. One power or square root.
// ---- the plant ------------------------------------------------------------------------------------------------------------------
const charact = (v, x) => { x = Math.min(1, Math.max(0, x)); const c = v.charact || 'linear'; if (c === 'equal') return x <= 0 ? 0 : Math.pow(v.R || 50, x - 1); if (c === 'quick') return Math.sqrt(x); return x; };PID controller output
64-bitThe output of a controller between 0 and 1.
e = ±(pv − sp) ÷ span; out = bias + K_c (e + I ÷ T_i + T_d (e − D_f) ÷ t_f), clamped to [0, 1]
The measurement and set point are scaled by the span, which is taken as 1 when the two span limits are equal. The derivative filter time is 0.1 × T_d plus 1e-9 unless given.
function ctrlOut(c, k, y) {
const span = (k.span[1] - k.span[0]) || 1, pv = (measure(c, y, k.pv) - k.span[0]) / span, sp = (k.sp - k.span[0]) / span;
const I = y[c.idx['ctrl.' + k.id + '.I']], Df = y[c.idx['ctrl.' + k.id + '.D']];
const eSigned = // direct: output rises as the measurement rises; reverse: as it falls
(k.action === 'direct' ? 1 : -1) * (pv - sp);
const tf = k.tf || 0.1 * (k.Td || 0) + 1e-9;
const dterm = k.Td ? (k.Td * (eSigned - Df)) / tf : 0;
const raw = (k.bias ?? 0.5) + k.Kc * (eSigned + (k.Ti ? I / k.Ti : 0) + dterm);
return { e: eSigned, raw, u: k.auto === false ? clamp01(k.out ?? 0) : clamp01(raw) };
}Flow through a valve from a tank
64-bitThe volume flow through a valve from the liquid heads on its two sides.
Δp = ρ g (h_from − h_to + Δz); q = k × f(x) × √(max(Δp, 0) ÷ ρ)
The head difference is formed before it is multiplied by ρ g, and a negative difference gives zero flow. A negative level in the source tank is taken as zero.
if (tank) { const rho = tank.rho; const hTo = v.to && (s.tanks || []).some((q2) => q2.id === v.to) ? levelOf(v.to) : 0; const dp = rho * G * (Math.max(y[X[v.from + '.h']], 0) - hTo + (v.dz || 0)); q = v.k * charact(v, x) * Math.sqrt(Math.max(dp, 0) / rho); outflow[v.from] = (outflow[v.from] || 0) + q; if (v.to) acc(v.to, q, null, null); }
else if (gas) { const P = y[X[v.from + '.P']], Pd = v.to && (s.gas || []).some((q2) => q2.id === v.to) ? y[X[v.to + '.P']] : (v.pOut ?? 1e5); q = v.k * charact(v, x) * Math.sqrt(Math.max(P - Pd, 0)); outflow[v.from] = (outflow[v.from] || 0) + q; if (v.to && (s.gas || []).some((q2) => q2.id === v.to)) acc(v.to, q, null, null); }Tank level and gas pressure rates
64-bitHow fast a tank level and a gas vessel pressure change.
dh/dt = (q_in − q_out) ÷ A; dP/dt = R T ÷ V × (n_in − n_out)
Inflows and outflows are each added up over the connected items and then subtracted. A tank at or below zero does not go lower and one at its limit does not go higher. R is the ten-digit constant from data.js.
for (const tk of s.tanks || []) { const h = y[X[tk.id + '.h']]; let d = ((inflow[tk.id] || 0) - (outflow[tk.id] || 0)) / tk.A; if (h <= 0 && d < 0) d = 0; if (tk.hmax != null && h >= tk.hmax && d > 0) d = 0; dy[X[tk.id + '.h']] = d; }
for (const g of s.gas || []) dy[X[g.id + '.P']] = (RGAS * g.T / g.V) * ((inflow[g.id] || 0) - (outflow[g.id] || 0));Stirred reactor balances
64-bitThe rate of change of each concentration and of the temperature in a constant-volume reactor.
r = k₀ exp(−E ÷ RT) Π C^n; dC/dt = q⁄V (C_in − C) + ν r; dT/dt = q⁄V (T_in − T) − ΔH r ÷ (ρ c_p) − UA (T − T_j) ÷ (ρ c_p V) + Q ÷ (ρ c_p V)
One Arrhenius exponential and one power per species; negative concentrations are taken as zero in the rate. The temperature balance adds four terms of mixed sign.
if (r.reaction) { const rc = r.reaction; rate = rc.k0 * Math.exp(-rc.E / (RGAS * T)); for (const [k, nn] of Object.entries(rc.orders || {})) rate *= Math.pow(Math.max(Cm[k] || 0, 0), nn); }
const Tfeed = q > 0 ? hin[r.id] / q : T, qq = Math.max(q, 0);
for (const k of r.comps) dy[X[r.id + '.C.' + k]] = qq / r.V * ((Cin[r.id] && qq > 0 ? (Cin[r.id][k] || 0) / qq : 0) - Cm[k]) + (r.reaction && r.reaction.nu[k] ? r.reaction.nu[k] * rate : 0);
const uH = mvOf('heater', r.id), Qh = (uH != null ? uH : 1) * (r.Q || 0), uJ = mvOf('jacket', r.id), Tj = uJ != null ? r.Tc + (uJ - 0.5) * (r.TcSpan || 0) : r.Tc;
dy[X[r.id + '.T']] = qq / r.V * (Tfeed - T) + (r.reaction ? (-r.reaction.dH * rate) / (r.rho * r.cp) : 0) - (r.UA ? r.UA * (T - Tj) / (r.rho * r.cp * r.V) : 0) + Qh / (r.rho * r.cp * r.V);Plant simulation
64-bitThe trends of every state and controller output from time zero to the end time.
integrate from each output time to the next (end ÷ 200 apart, or the next scripted event) with rtol 1e-6, atol 1e-9
The clock is set to each target time, not advanced by adding steps, so output times do not drift. Each interval is integrated to the default tolerances, which hold the local error per step, not the error at the end of the run.
/**
* → { t[], series: { '<name>': [..] }, outputs: { '<ctrl id>': [u…] }, steps, rejected, problems }
* opts: { tEnd, dt (output spacing, default tEnd / 200), rtol, atol }
*/
export function simulate(spec0, { tEnd, dt = null, rtol = 1e-6, atol = 1e-9 } = {}) {
const c = compile(spec0), problems = [];
if (!(tEnd > 0)) throw new Error('simulate needs tEnd > 0');
const step = dt || tEnd / 200;
const events = [...(c.spec.script || [])].sort((a, b) => a.t - b.t);
let y = y0Of(c), t = 0, steps = 0, rejected = 0;
const times = [0], rows = [y.slice()], outs = [];
const f = (tt, yy) => derivative(c, tt, yy);
const record = (tt, yy) => { f(tt, yy); outs.push(Object.fromEntries((c.spec.controllers || []).map((k) => [k.id, c._ctr[k.id].u]))); };
record(0, y);
let ei = 0;
while (t < tEnd - 1e-12) {
while (ei < events.length && events[ei].t <= t + 1e-12) { setTarget(c.spec, events[ei].set.target, events[ei].set.value); ei++; }
const next = Math.min(tEnd, t + step, ei < events.length ? events[ei].t : Infinity);
const r = stiff(f, t, y, next, { rtol, atol, h0: (next - t) / 4 });
y = r.y; t = next; steps += r.steps; rejected += r.rejected;
times.push(t); rows.push(y.slice()); record(t, y);
}
const series = {}; c.names.forEach((nme, i) => { series[nme] = rows.map((row) => row[i]); });
for (const tk of c.spec.tanks || []) if (series[tk.id + '.h'].some((v) => v < -1e-9)) problems.push({ id: 'process:dyn.dry:' + tk.id, code: 'dyn.dry', severity: 'error', text: `Tank ${tk.id} ran dry in the simulation.`, about: [{ kind: 'unit', id: tk.id }], source: { kind: 'engine', name: 'process' } });
for (const tk of c.spec.tanks || []) if (tk.hmax != null && series[tk.id + '.h'].some((v) => v >= tk.hmax * (1 - 1e-9))) problems.push({ id: 'process:dyn.overflow:' + tk.id, code: 'dyn.overflow', severity: 'error', text: `Tank ${tk.id} overflowed in the simulation.`, about: [{ kind: 'unit', id: tk.id }], source: { kind: 'engine', name: 'process' } });
const outputs = {}; for (const k of c.spec.controllers || []) outputs[k.id] = outs.map((o) => o[k.id]);
return { t: times, series, outputs, steps, rejected, problems };
}Process model from a step test
64-bitGain, time constant and dead time of a first-order model fitted to a step response.
K = Δy ÷ Δu; τ = 1.5 (t₆₃ − t₂₈); θ = max(0, t₆₃ − τ)
The two times are found by straight-line interpolation between samples at 28.3 % and 63.2 % of the change, the rounded fractions of the method. A response that never crosses a level gives a result that is not a number.
/** Smith's two-point method on a step response (y at time t after a step of size du in the input): K, tau, theta. */
export function fopdtFromStep(t, y, du) {
const y0 = y[0], yf = y[y.length - 1], dy = yf - y0, at = (frac) => { const target = y0 + frac * dy; for (let i = 1; i < t.length; i++) if ((y[i] - target) * (y[i - 1] - target) <= 0 && y[i] !== y[i - 1]) return t[i - 1] + (target - y[i - 1]) / (y[i] - y[i - 1]) * (t[i] - t[i - 1]); return NaN; };
const t28 = at(0.283), t63 = at(0.632), tau = 1.5 * (t63 - t28), theta = Math.max(0, t63 - tau);
return { K: dy / du, tau, theta };
}Controller tuning rules
64-bitP, PI or PID settings from the fitted gain, time constant and dead time.
Ziegler–Nichols PI: K_c = 0.9 τ ÷ (K θ), T_i = 3.33 θ; IMC PI: K_c = (τ + θ⁄2) ÷ (K (λ + θ⁄2)), T_i = τ + θ⁄2; SIMC: K_c = τ ÷ (K (τ_c + θ)), T_i = min(τ, 4 (τ_c + θ))
Rules of thumb with typed constants such as 0.9 and 3.33. The dead time is floored at 1e-9 so the divisions do not fail for a process with none.
export function tune(method, { K, tau, theta }, { type = 'PI', lambda = null } = {}) {
const th = Math.max(theta, 1e-9), r = tau / (K * th), x = th / tau;
if (method === 'zn') return type === 'P' ? { Kc: r } : type === 'PI' ? { Kc: 0.9 * r, Ti: 3.33 * th } : { Kc: 1.2 * r, Ti: 2 * th, Td: 0.5 * th };
if (method === 'cohen-coon') return type === 'P' ? { Kc: r * (1 + x / 3) } : type === 'PI' ? { Kc: r * (0.9 + x / 12), Ti: th * (30 + 3 * x) / (9 + 20 * x) } : { Kc: r * (4 / 3 + x / 4), Ti: th * (32 + 6 * x) / (13 + 8 * x), Td: 4 * th / (11 + 2 * x) };
if (method === 'imc') { const l = lambda ?? Math.max(th, 0.2 * tau); return type === 'PI' ? { Kc: (tau + th / 2) / (K * (l + th / 2)), Ti: tau + th / 2 } : { Kc: (tau + th / 2) / (K * (l + th / 2)), Ti: tau + th / 2, Td: tau * th / (2 * tau + th) }; }
if (method === 'simc') { const tc = lambda ?? th; return { Kc: tau / (K * (tc + th)), Ti: Math.min(tau, 4 * (tc + th)), ...(type === 'PID' ? { Td: 0 } : {}) }; }
throw new Error('unknown tuning method ' + method);
}In detail
- Method
- Implicit Euler with step-doubling (Richardson) extrapolation and a frozen finite-difference Jacobian; valve flow by the square-root law; first-order-plus-dead-time fit by the two-point method; Ziegler–Nichols, Cohen–Coon, IMC and SIMC tuning.
- Number types
- 64-bit numbers in plain arrays.
- Tolerances
- Step error: |two half steps − one full step| ÷ (atol + rtol × max |y|) ≤ 1, with rtol 1e-6 and atol 1e-9 by default. Newton inside a step: max |Δx| < 1e-11 × max |x| + 1e-14, 12 passes. Jacobian step 1e-7 × max(|y|, 1e-3). Step factor between 0.2 and 4. The solver stops with an error when the step falls below 1e-15 × max(1, |t|) or the steps run out.
- Limits
- The simulation's output is spaced at the end time ÷ 200 unless set. A tank is reported dry below −1e-9 m and overflowing at its limit × (1 − 1e-9).
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.