Resources Numerical accuracy Calculation policy
Runs the matrix products, dot products and vector updates of the simulation solvers on several threads, and the conjugate-gradient solves built on them.
Every vector is a Float64Array on shared memory and all arithmetic is 64-bit. Each output entry of a matrix product is written by one thread, with its terms added in element or row order. Dot products are summed in fixed chunks of 4096 entries and the chunk sums are then added in order, so the result does not depend on the number of threads. The accuracy of a solve is set by its stopping tolerance on the relative residual.
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.
Tests tools/partcraft-test/test-threads.mjs
Each thread sums its share of a dot product in chunks of 4096 entries.
P_c = Σ aᵢ·bᵢ over chunk c
Each chunk is a plain running sum of 64-bit products. The chunk boundaries are fixed by the vector length, not by the thread count.
// 0 dot: P ← chunks of Σ a·b I: a, b, n
(t, T) => { const a = B[I[0]], b = B[I[1]]; chunked(I[2], t, T, (i0, i1) => { let s = 0; for (let i = i0; i < i1; i++) s += a[i] * b[i]; return s; }); },Adds the chunk sums of a dot product in chunk order on the coordinating thread.
a·b = Σ P_c
The partials are a Float64Array and are added in order with no compensation. Because chunks and order are fixed, the result is the same for any thread count; it can differ from one running sum over the whole vector by rounding.
sum(nc) { let s = 0; for (let c = 0; c < nc; c++) s += this.P[c]; return s; }Advances the solution and the residual by one step and sums the new residual's square.
x += α·p; r −= α·q; P_c = Σ rᵢ²
The residual is updated recursively and its norm is taken from the updated values; the kernel never recomputes b − Kx. Each entry is touched by one thread.
// 1 cg update: x += α p, r −= α q; P ← chunks of r·r I: x, p, r, q, n F: α
(t, T) => { const x = B[I[0]], p = B[I[1]], r = B[I[2]], q = B[I[3]], a = F[0]; chunked(I[4], t, T, (i0, i1) => { let s = 0; for (let i = i0; i < i1; i++) { x[i] += a * p[i]; r[i] -= a * q[i]; s += r[i] * r[i]; } return s; }); },y = A·x for a matrix stored row by row, with the rows split between the threads.
y_r = Σ_k v_k · x[col_k]
Each row is one running 64-bit sum in stored column order, written by one thread, so the product equals the one-thread loop bit for bit.
// 3 CSR y = A x over a row partition I: rowPtr, cols, vals, x, y, part
(t, T) => { const rp = B[I[0]], cl = B[I[1]], v = B[I[2]], x = B[I[3]], y = B[I[4]], [r0, r1] = part(I[5], t); for (let r = r0; r < r1; r++) { let s = 0; for (let k = rp[r]; k < rp[r + 1]; k++) s += v[k] * x[cl[k]]; y[r] = s; } },y = K·x for a voxel model without assembling K: each node gathers the contributions of its elements.
y_n = Σ over elements at n of s_e · (K_e · x_e) rows of n, plus the 3 × 3 penalty block of a constrained node
For each node the element contributions are added in element order, each as a 24-term 64-bit sum scaled by the element's stiffness factor. The element number and corner are packed in one integer and split with >> 3 and & 7; those are indexes, not measured values. An element with scale 0 is skipped.
// 5 voxel K x, gathered per node over its elements in element order (= fea.js applyK's scatter, bitwise); explicit levels
// (mg.js) take one 24×24 matrix per element; an optional penalty block per constrained node is added last (mg level 0)
// I: nodePtr, nodeElem (m·8 + corner), node0, off, Ke, scale|−1, x, y, nodes(nn), explicit(0|1), consIdx|−1, pen|−1, nodePart F: kScale
(t, T) => {
const np = B[I[0]], ne = B[I[1]], n0 = B[I[2]], off = B[I[3]], K = B[I[4]], sc = I[5] >= 0 ? B[I[5]] : null, x = B[I[6]], y = B[I[7]], expl = I[9] === 1;
const ci = I[10] >= 0 ? B[I[10]] : null, pen = I[11] >= 0 ? B[I[11]] : null, kS = F[0], [a0, a1] = part(I[12], t);
const xe = new Float64Array(24);
for (let n = a0; n < a1; n++) {
let y0 = 0, y1 = 0, y2 = 0;
for (let k = np[n]; k < np[n + 1]; k++) {
const m = ne[k] >> 3, c = ne[k] & 7, s = expl ? 1 : (sc ? sc[m] * kS : kS);
if (s === 0) continue;
const b0 = n0[m];
for (let q = 0; q < 8; q++) { const d = 3 * (b0 + off[q]); xe[3 * q] = x[d]; xe[3 * q + 1] = x[d + 1]; xe[3 * q + 2] = x[d + 2]; }
const base = expl ? m * 576 : 0, r = base + 3 * c * 24;
let t0 = 0, t1 = 0, t2 = 0;
for (let j = 0; j < 24; j++) { const v = xe[j]; t0 += K[r + j] * v; t1 += K[r + 24 + j] * v; t2 += K[r + 48 + j] * v; }
if (expl) { y0 += t0; y1 += t1; y2 += t2; } else { y0 += s * t0; y1 += s * t1; y2 += s * t2; }
}
if (ci && ci[n] >= 0) { const p = 9 * ci[n], d = 3 * n, a = x[d], b = x[d + 1], cc = x[d + 2]; y0 += pen[p] * a + pen[p + 1] * b + pen[p + 2] * cc; y1 += pen[p + 3] * a + pen[p + 4] * b + pen[p + 5] * cc; y2 += pen[p + 6] * a + pen[p + 7] * b + pen[p + 8] * cc; }
y[3 * n] = y0; y[3 * n + 1] = y1; y[3 * n + 2] = y2;
}
},Multiplies each node's three values by its 3 × 3 block, or adds a damped correction ω·M·(b − r).
dst = M·src or dst += ω · M · (b − src)
Three 3-term 64-bit sums per node; each node belongs to one thread.
// 6 3×3 blocks per node: dst = M src (+= when I[4] = 1: x += ω M (b − r) for the mg smoother) I: M, src, dst, nn, mode, b|−1 F: ω
(t, T) => {
const M = B[I[0]], src = B[I[1]], dst = B[I[2]], [n0, n1] = span(I[3], t, T), mode = I[4];
if (mode === 0) for (let nd = n0; nd < n1; nd++) { const d = 3 * nd, m = 9 * nd, x = src[d], y = src[d + 1], w = src[d + 2]; dst[d] = M[m] * x + M[m + 1] * y + M[m + 2] * w; dst[d + 1] = M[m + 3] * x + M[m + 4] * y + M[m + 5] * w; dst[d + 2] = M[m + 6] * x + M[m + 7] * y + M[m + 8] * w; }
else { const b = B[I[5]], om = F[0]; for (let nd = n0; nd < n1; nd++) { const d = 3 * nd, m = 9 * nd, e0 = b[d] - src[d], e1 = b[d + 1] - src[d + 1], e2 = b[d + 2] - src[d + 2]; dst[d] += om * (M[m] * e0 + M[m + 1] * e1 + M[m + 2] * e2); dst[d + 1] += om * (M[m + 3] * e0 + M[m + 4] * e1 + M[m + 5] * e2); dst[d + 2] += om * (M[m + 6] * e0 + M[m + 7] * e1 + M[m + 8] * e2); } }
},Sets the values of inactive nodes to zero and applies each constrained node's 3 × 3 projector.
v_n = 0 (inactive); v_n ← P_n · v_n (constrained)
Zeroing is exact; the projection is three 3-term 64-bit sums per constrained node from a Float64Array of projectors.
// 7 project (fea.js projector): inactive nodes 0, constrained nodes v ← P v I: v, active, nn, consNodes, P, ncons
(t, T) => {
const v = B[I[0]], act = B[I[1]], [n0, n1] = span(I[2], t, T);
for (let n = n0; n < n1; n++) if (!act[n]) { v[3 * n] = 0; v[3 * n + 1] = 0; v[3 * n + 2] = 0; }
if (I[5] > 0) { const cn = B[I[3]], P = B[I[4]], [m0, m1] = span(I[5], t, T); for (let m = m0; m < m1; m++) { const d = 3 * cn[m], p = 9 * m, x = v[d], y = v[d + 1], z = v[d + 2]; v[d] = P[p] * x + P[p + 1] * y + P[p + 2] * z; v[d + 1] = P[p + 3] * x + P[p + 4] * y + P[p + 5] * z; v[d + 2] = P[p + 6] * x + P[p + 7] * y + P[p + 8] * z; } }
},Carries a residual from a fine grid to the grid of half the resolution.
r_c(N) = Σ over the 27 fine neighbors of w · r_f; w = product per axis of 1 (same node) or ½ (between nodes)
The weights are products of 1 and ½, so they are exact in binary. Fine nodes are visited in increasing index order for each coarse node, and inactive fine nodes are left out.
// 9 mg restriction gathered per coarse node (fine nodes in increasing order = mg.js restrict's scatter order)
// I: rf, rc, fineActive, FNX, FNY, FNZ (fine node counts per axis), CNX, CNY, CNZ (coarse node range by t)
(t, T) => {
const rf = B[I[0]], rc = B[I[1]], act = B[I[2]], FX = I[3], FY = I[4], FZ = I[5], CX = I[6], CY = I[7], CZ = I[8], [c0, c1] = span(CX * CY * CZ, t, T);
for (let N = c0; N < c1; N++) {
const I0 = N % CX, J0 = Math.floor(N / CX) % CY, K0 = Math.floor(N / (CX * CY));
let s0 = 0, s1 = 0, s2 = 0;
for (let k = 2 * K0 - 1; k <= 2 * K0 + 1; k++) {
if (k < 0 || k >= FZ) continue; const wk = k & 1 ? 0.5 : 1;
for (let j = 2 * J0 - 1; j <= 2 * J0 + 1; j++) {
if (j < 0 || j >= FY) continue; const wj = j & 1 ? 0.5 : 1;
for (let i = 2 * I0 - 1; i <= 2 * I0 + 1; i++) {
if (i < 0 || i >= FX) continue;
const n = i + FX * (j + FY * k); if (!act[n]) continue;
const w = (i & 1 ? 0.5 : 1) * wj * wk;
s0 += w * rf[3 * n]; s1 += w * rf[3 * n + 1]; s2 += w * rf[3 * n + 2];
}
}
}
rc[3 * N] = s0; rc[3 * N + 1] = s1; rc[3 * N + 2] = s2;
}
},Adds a coarse-grid correction to the fine grid by linear interpolation.
x_f(n) += w · Σ x_c over the 1, 2, 4 or 8 coarse nodes around n; w = ½ per axis in which n lies between coarse nodes
The weights are exact powers of ½. The halving of the fine index (i >> 1) and the odd test (i & 1) act on grid indexes only.
// 10 mg prolongation: xf += P xc per fine node (the order of mg.js eachTransfer) I: xc, xf, fineActive, FNX, FNY, FNZ, CNX, CNY
(t, T) => {
const xc = B[I[0]], xf = B[I[1]], act = B[I[2]], FX = I[3], FY = I[4], FZ = I[5], CX = I[6], CY = I[7], [f0, f1] = span(FX * FY * FZ, t, T);
for (let n = f0; n < f1; n++) {
if (!act[n]) continue;
const i = n % FX, j = Math.floor(n / FX) % FY, k = Math.floor(n / (FX * FY)), I0 = i >> 1, J0 = j >> 1, K0 = k >> 1, oi = i & 1, oj = j & 1, ok = k & 1;
const w = (oi ? 0.5 : 1) * (oj ? 0.5 : 1) * (ok ? 0.5 : 1);
for (let c = 0; c <= ok; c++) for (let b = 0; b <= oj; b++) for (let a = 0; a <= oi; a++) { const N = (I0 + a) + CX * ((J0 + b) + CY * (K0 + c)); xf[3 * n] += w * xc[3 * N]; xf[3 * n + 1] += w * xc[3 * N + 1]; xf[3 * n + 2] += w * xc[3 * N + 2]; }
}
},y = s·(G ⊗ I₃)·x from per-element blocks, used by the buckling study.
y_n = s · Σ over elements at n of Σ_b G_e[a, b] · x_b
Gathered per node in element order with 64-bit sums. The element number is recovered with Math.floor of an index division, not from a measured value.
// 11 y = s·(G ⊗ I₃) x from per-element npe×npe blocks (buckling.js geometricOperator), gathered per node in element order
// I: nodePtr, nodeElem (e·npe + a), conn, G (ne·npe²), x, y, npe, nodePart F: s
(t, T) => {
const np = B[I[0]], ne = B[I[1]], conn = B[I[2]], G = B[I[3]], x = B[I[4]], y = B[I[5]], npe = I[6], s = F[0], [a0, a1] = part(I[7], t), q = npe * npe;
for (let n = a0; n < a1; n++) {
let y0 = 0, y1 = 0, y2 = 0;
for (let k = np[n]; k < np[n + 1]; k++) {
const e = Math.floor(ne[k] / npe), a = ne[k] - e * npe, g0 = e * q + a * npe, c0 = e * npe;
let t0 = 0, t1 = 0, t2 = 0;
for (let b = 0; b < npe; b++) { const g = G[g0 + b], nb = 3 * conn[c0 + b]; t0 += g * x[nb]; t1 += g * x[nb + 1]; t2 += g * x[nb + 2]; }
y0 += s * t0; y1 += s * t1; y2 += s * t2;
}
y[3 * n] = y0; y[3 * n + 1] = y1; y[3 * n + 2] = y2;
}
},Solves K·u = f for a voxel stress model on the thread pool.
res = ‖r‖ ÷ ‖b‖; stop when res < tol; α = (r·z) ÷ (p·Kp); β = (r·z)_new ÷ (r·z)_old
Stops when the relative residual is under tol (1e-8 by default) or after max(2000, 600 × ∛elements) iterations, and reports converged only in the first case. A zero right-hand side returns zeros at once, so the division by ‖b‖ is guarded; the loop also leaves when p·Kp is not positive. The solution is copied out as a Float64Array.
/**
* fea.js pcg on the pool: the same projected PCG (same operator, preconditioner and stopping test; dot products summed in chunks).
* mg: an mg.js preconditioner (threaded when built with the pool); else M: node-block Jacobi (9 per node). → { u (a plain copy), … }
*/
export function voxelPCG(P, model, cons, f, { x0 = null, tol = 1e-8, maxIter = 0, scale = null, mg = null, M = null, onProgress = null, cancelled = null, label = 'Solving' } = {}) {
const n = 3 * model.nn, proj = projectOperator(P, model, cons), applyK = hexOperator(P, model, { scale });
const u = P.f64(n), r = P.f64(n), z = P.f64(n), p = P.f64(n), q = P.f64(n), b = P.f64(n), tmp = mg ? P.f64(n) : null, Ms = mg ? null : P.share(M);
try {
if (x0) u.set(x0);
proj(u);
b.set(f); proj(b);
const bn = Math.sqrt(P.dot(b, b));
if (bn === 0) return { u: new Float64Array(n), iterations: 0, residual: 0, converged: true, threads: P.n };
applyK(u, q);
P.run(K.AMINUS, [b._pid, q._pid, n]); P.copy(r, q); // r = b − K u
proj(r);
const prec = (src, dst) => {
if (mg) { P.copy(tmp, src); proj(tmp); mg.apply(tmp, dst); return proj(dst); }
P.run(K.BLOCK3, [Ms._pid, src._pid, dst._pid, model.nn, 0, -1]);
return proj(dst);
};
prec(r, z); P.copy(p, z);
let rz = P.dot(r, z), rn = P.dot(r, r);
const nc = P.partials(n);
const cap = maxIter || Math.max(2000, 20 * Math.cbrt(model.elems.length) * 30);
let it = 0, res = 1, t0 = Date.now();
for (; it < cap; it++) {
res = Math.sqrt(rn) / bn;
if (res < tol) break;
if (onProgress && Date.now() - t0 > 200) { t0 = Date.now(); onProgress(Math.min(0.99, Math.log(res) / Math.log(tol)), `${label}: iteration ${it}, residual ${res.toExponential(1)}`); } // i18n-ignore: progress label, never shown
if (cancelled && (it & 15) === 0 && cancelled()) throw new Error('Cancelled');
applyK(p, q); proj(q);
const pq = P.dot(p, q);
if (!(pq > 0)) break;
const a = rz / pq;
P.run(K.CG, [u._pid, p._pid, r._pid, q._pid, n], [a]); rn = P.sum(nc);
prec(r, z);
const rz2 = P.dot(r, z);
const beta = rz2 / rz; rz = rz2;
P.run(K.XPBY, [p._pid, z._pid, n], [beta]);
}
return { u: Float64Array.from(u), iterations: it, residual: res, converged: res < tol, threads: P.n };
} finally { P.free([u, r, z, p, q, b, tmp, Ms === M ? null : Ms].filter(Boolean)); }
}Solves A·x = b for an assembled sparse matrix on the thread pool.
res = ‖r‖ ÷ ‖b‖; stop when res < tol; α = (r·z) ÷ (p·Ap)
Stops when the relative residual is under tol (1e-10 by default) or after max(1000, 100 × round(∛n)) iterations. The preconditioner is applied on one thread. A zero right-hand side returns zeros, and the loop leaves when p·Ap is not positive. No other file of the released app calls this function.
/** thermal.js csrPcg on the pool: the matrix products and vector updates threaded; the preconditioner M (IC(0)) stays serial. */
export function csrPCG(P, pat, val, b, { x0 = null, tol = 1e-10, maxIter = 0, cancelled = null } = {}, M) {
const n = pat.n, A = csrOperator(P, pat, val);
const x = P.f64(n), r = P.f64(n), z = P.f64(n), p = P.f64(n), q = P.f64(n), bs = P.share(b);
try {
const bn = Math.sqrt(P.dot(bs, bs));
if (bn === 0) return { x: new Float64Array(n), iterations: 0, residual: 0, converged: true };
if (x0) x.set(x0);
A(x, q);
P.copy(r, q); P.run(K.AMINUS, [bs._pid, r._pid, n]);
M.apply(r, z); P.copy(p, z);
let rz = P.dot(r, z), rn = P.dot(r, r);
const nc = P.partials(n);
const cap = maxIter || Math.max(1000, 10 * Math.round(Math.cbrt(n)) * 10);
let it = 0, res = 1;
for (; it < cap; it++) {
res = Math.sqrt(rn) / bn;
if (res < tol) break;
if (cancelled && (it & 31) === 0 && cancelled()) throw new Error('Cancelled');
A(p, q);
const pq = P.dot(p, q);
if (!(pq > 0)) break;
const a = rz / pq;
P.run(K.CG, [x._pid, p._pid, r._pid, q._pid, n], [a]); rn = P.sum(nc);
M.apply(r, z);
const rz2 = P.dot(r, z);
const beta = rz2 / rz; rz = rz2;
P.run(K.XPBY, [p._pid, z._pid, n], [beta]);
}
return { x: Float64Array.from(x), iterations: it, residual: res, converged: res < tol };
} finally { P.free([x, r, z, p, q, bs === b ? null : bs].filter(Boolean)); }
}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.