diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 34ea1d4..3baf87d 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -31,3 +31,6 @@ jobs: - run: node witness-map/table.mjs - run: node witness-map/table.mjs 1e-6 cryo - run: node notes/bh-lensing.mjs + - run: node turbulence/table.mjs quick + - run: node turbulence/edges.mjs quick + - run: node turbulence/edges.mjs 256 diff --git a/CORRECTIONS.md b/CORRECTIONS.md index 2fb953e..0c545e0 100644 --- a/CORRECTIONS.md +++ b/CORRECTIONS.md @@ -75,11 +75,11 @@ Claims made during the exploration that did not survive it, kept so the record s ## 16 October 2026: the gravitational atom (issue #4) -12. **Atoms as little solar systems.** Claimed (README, and issue #4): "Atoms are not little solar systems ... Except - in one case, where they are", the one place where that picture "becomes literally true". Too strong. In the - gravitational atom the light particles form a wave cloud with hydrogen's shapes and levels; nothing in it moves - on an orbit. What carries over from a solar system is the force, gravity, and for a heavy hole the size: in - Baumann, Chia and Porto's 2019 example the 211 cloud is most likely found 1.4 AU out. The README now says +12. **Atoms as little solar systems.** Claimed (README, and issue #4): "Atoms are not little solar systems ... + Except in one case, where they are", the one place where that picture "becomes literally true". Too strong. In + the gravitational atom the light particles form a wave cloud with hydrogen's shapes and levels; nothing in it + moves on an orbit. What carries over from a solar system is the force, gravity, and for a heavy hole the size: + in Baumann, Chia and Porto's 2019 example the 211 cloud is most likely found 1.4 AU out. The README now says gravity builds a real atom. 13. **"The same energy levels as a hydrogen atom."** Means the same pattern, falling as 1/n², not the same energies. In the same example the levels are 2×10¹⁹ times closer together than hydrogen's, and the pattern @@ -99,16 +99,16 @@ Claims made during the exploration that did not survive it, kept so the record s ## 16 October 2026: black hole and wave dark matter lensing (issue #1) 15. **Superradiance as a second route.** Issue #1 said a spinning black hole amplifying very light particles into a - cloud was not ruled out as a way for a galaxy's central black hole to shape wave dark matter where lensing looks. - For a 10⁻²² eV particle and a 10⁸ to 10⁹ M☉ black hole it is ruled out: the cloud would take 10²⁶ to 10³⁴ years - to grow, and at its radius the galaxy, not the black hole, holds the particles. Any cloud these black holes can - grow within the age of the universe, at any particle mass, peaks within about 0.1 pc. See + cloud was not ruled out as a way for a galaxy's central black hole to shape wave dark matter where lensing + looks. For a 10⁻²² eV particle and a 10⁸ to 10⁹ M☉ black hole it is ruled out: the cloud would take 10²⁶ to 10³⁴ + years to grow, and at its radius the galaxy, not the black hole, holds the particles. Any cloud these black + holes can grow within the age of the universe, at any particle mass, peaks within about 0.1 pc. See [notes/bh-lensing.md](notes/bh-lensing.md). ## 18 October 2026: the black hole lensing note revisited (issue #28) -16. **SDP.81 as the anomaly lens with a bound on its black hole.** Claimed (notes/bh-lensing.md, issue #1): "Only one - anomaly lens, SDP.81, has even a bound on its black hole". The note paired a 10^(8.96±0.12) M☉ dark lump at +16. **SDP.81 as the anomaly lens with a bound on its black hole.** Claimed (notes/bh-lensing.md, issue #1): "Only + one anomaly lens, SDP.81, has even a bound on its black hole". The note paired a 10^(8.96±0.12) M☉ dark lump at 6.9σ (Hezaveh et al. 2016) with a black hole above 3×10⁸ M☉ from the missing central image. That no longer holds. Stacey et al. (A&A, 2025) re-analysed the same ALMA data. They find the lump most likely spurious, produced by poorly fitting multipoles, and no evidence for either previously reported lump, only for extra @@ -135,3 +135,43 @@ Claims made during the exploration that did not survive it, kept so the record s 2.79 µm diameter they used, slip-corrected Stokes drag gives 175 µs at 2.75 kPa against their 147.3 µs; `tests/quantum-paths.mjs` records that as a disagreement. (The first draft on #25 used the nominal 3 µm, which hid the gap, and called the trapped path smooth, then rough, then smooth; review there caught both.) + +## 18 October 2026: turbulence (issue #19) + +19. **Issue #19's numbers, written from recall.** Checked against search text, and in #34's review against the + papers (see [turbulence](turbulence/)), two did not survive and five needed qualifying. + - Did not survive: "σ = 0.72 S^0.5 km/s". Solomon et al.'s fit is 1.0 S^0.5; with their effective radius + R = 1.92 S it becomes 0.72 R^0.5 (this study's conversion), so 0.72 goes with R. The 0.07 pc sonic scale it + gave is a radius, 0.15 pc across. + - Did not survive: Larson's σ ∝ L^0.38 against Solomon et al.'s 0.5 as "a change of cascade". Over most of + Larson's range his relation is supersonic (non-thermal Mach 1.0 at 0.1 pc to 19 at 100 pc), and the four + objects he plots at 0.1 to 0.2 pc scatter from Mach 0.44 (TMC2) to 2.45 (L1407/L1257), with only TMC1 and TMC2 + subsonic in his own words, so his 0.38 is not the subsonic cascade's. That the difference lies in his sample + and size definition is this study's inference; Larson read it as supersonic turbulence, or subsonic turbulence + in warmer gas. + - Qualified: "the spectral break near 0.6 Hz". V/(2π d_i) gives 0.63 Hz, but the fluctuations are elongated + along the field, so the frequency carries sin θ_BV: 0.46 Hz at the 47° Parker angle. Observed breaks spread + from 0.1 to 1 Hz (Chen et al. 2014), and Leamon et al.'s intervals average near 0.38 Hz, so both values lie + inside the spread and the data do not choose. + - Qualified: "a kinetic-range slope near −2.8 (Bruno and Carbone)". −2.8 is one measurement (Alexandrova et al. + 2009); Bruno and Carbone give about −7/3 as typical (as read on p. 146 in #34's re-review), and measured + slopes run from about −2 to −4 (as Alexandrova et al. and Bruno and Carbone summarise them; Smith et al. 2006 + find −1 to −4). + - Qualified: "the Kolmogorov scale is about a millimetre in the atmosphere and upper ocean". Right as an order + of magnitude, with a factor of 10 either way: 0.3 mm in strong wind to over a centimetre on a quiet night, + and 2 to 3 mm at the one upper-ocean station checked. + - Qualified: "the Kolmogorov constant (~1.5)". Fair for the 3D spectrum. Sreenivasan's (1995) compilation gives + 0.53 (scatter 0.055) for the 1D spectrum, 1.62 ± 0.17 in 3D. + - Qualified: Nastrom and Gage's k⁻³ range "above ~500 km". The transition is quoted between 400 and 1,000 km, + most often near 600. +20. **"A fractional-Brownian estimate gives 7/6" for cloud edges.** The wrong formula for single loops without + holes. 7/6 = 2 − H is the dimension of all the contours at one level together, the relation Stutzki et al. + (1998) take from Voss, so the issue followed the literature. A perimeter-area relation over separate clouds, if + holes are not counted, measures single loops, (3 − H)/2 (Kondev and Henley 1995, a conjecture exact only at + H = 0 and 1), which is 13/12 for projected Kolmogorov density. Computed on Gaussian surfaces, the loops follow + (3 − H)/2: across six spectra their dimension falls with the effective roughness H_eff at a slope of −0.53 + ± 0.05 (−0.44 against the nominal H), not −1. They also sit 0.02 to 0.04 above it, depending on where H_eff is + measured, for reasons not found. The issue's conclusion survives: projected Kolmogorov density gives 1.15 there, + not 1.36. Whether the clouds' 1.36 matches the laboratory's 1.3 to 1.4 for a shared reason is not shown. The + laboratory measures sections of turbulent surfaces (7/3 in 3D for the interface), the clouds a projection of + supersonic gas, and the two only meet near 4/3. diff --git a/README.md b/README.md index 152e37a..627d845 100644 --- a/README.md +++ b/README.md @@ -26,6 +26,16 @@ with Claude (Anthropic), for [Foundation0](https://github.com/Foundation-0/f0). - **Cosmic structure is fractal only up to a point.** Galaxies cluster in a fractal-like way on small scales, but beyond about 80 to 300 million parsecs (h⁻¹ Mpc, 80 to 300) they are spread almost evenly (Scrimgeour et al., 2012). The nesting stops, the way a coastline's roughness stops at the size of molecules. +- **Turbulence repeats, but not exactly.** Between the size where a flow is stirred and the size where viscosity + takes over, the same energy flux passes through every scale. In a cascade model it holds to 0.2% over the first + nine of 21 shells (2.4 decades), and a third-order law that follows from it to 0.1%. The bursts do not scale + exactly: the distribution's tails fatten from scale to scale (intermittency). Where the cascade ends depends on + the medium: about a millimetre in air and water, give or take a factor of 10; the ion scales (about 100 km) in + the solar wind, where collisions are too rare for viscosity and a steeper kinetic cascade takes over. In + molecular clouds the supersonic range ends at the sonic scale (about 0.1 pc), and a shallower subsonic cascade + continues. The fractal edges of the Taurus clouds (dimension 1.36) agree in number with laboratory turbulent + surfaces, but the two measure different things, and across clouds the dimension runs from about 1.2 to 1.6 + ([turbulence](turbulence/)). Some repetition across scales has to be tuned, and some sets its own range. At a critical point, such as water's liquid-vapour critical point or the Curie temperature of a magnet whose spins point only up or down, fluctuations diff --git a/tests/turbulence-edges.mjs b/tests/turbulence-edges.mjs new file mode 100644 index 0000000..ddc6aed --- /dev/null +++ b/tests/turbulence-edges.mjs @@ -0,0 +1,129 @@ +// Checks the contour-loop machinery behind the cloud-edge question (turbulence/contours.js): the FFT, the random +// surfaces, the marching squares, and which dimension a perimeter-area fit measures. Run: npm test. +// Each check states the source, the expected value, the model's value and the tolerance. +import { createRequire } from 'node:module'; +const require = createRequire(import.meta.url); +const C = require('../turbulence/contours.js'); + +let fails = 0; +function check(name, ok, detail, src) { + console.log(`${ok ? 'PASS' : 'FAIL'} ${name}\n ${detail} ${src}`); + if (!ok) fails++; +} +const e1 = x => x.toExponential(1), f3 = x => x.toFixed(3), near = (a, b, t) => Math.abs(a - b) <= t; +const SEED = 2026; + +// 1. The pieces. FFT against the direct sum, Parseval and the round trip. +const rnd = C.gaussian(C.stream(SEED, 0, 1 << 20)), n1 = 64, xr = Float64Array.from({ length: n1 }, rnd), xi = Float64Array.from({ length: n1 }, rnd); +const yr = Float64Array.from(xr), yi = Float64Array.from(xi); C.fft(yr, yi, false); +let dft = 0; +for (let k = 0; k < n1; k++) { + let sr = 0, si = 0; + for (let t = 0; t < n1; t++) { const th = -2 * Math.PI * k * t / n1; sr += xr[t] * Math.cos(th) - xi[t] * Math.sin(th); si += xr[t] * Math.sin(th) + xi[t] * Math.cos(th); } + dft = Math.max(dft, Math.abs(sr - yr[k]), Math.abs(si - yi[k])); +} +const L2 = 128, N2 = L2 * L2, zr = Float64Array.from({ length: N2 }, rnd), zi = Float64Array.from({ length: N2 }, rnd), wr = Float64Array.from(zr), wi = Float64Array.from(zi); +C.fft2(wr, wi, L2, false); +let ex = 0, eX = 0; for (let i = 0; i < N2; i++) { ex += zr[i] ** 2 + zi[i] ** 2; eX += wr[i] ** 2 + wi[i] ** 2; } +C.fft2(wr, wi, L2, true); +let trip = 0; for (let i = 0; i < N2; i++) trip = Math.max(trip, Math.abs(wr[i] - zr[i]), Math.abs(wi[i] - zi[i])); +check('FFT equals the direct DFT sum (N = 64); 2D FFT keeps Parseval and inverts itself (128 x 128)', dft < 1e-12 && Math.abs(ex - eX / N2) / ex < 1e-12 && trip < 1e-12, + `DFT ${e1(dft)}, Parseval ${e1(Math.abs(ex - eX / N2) / ex)}, round trip ${e1(trip)} (each under 1e-12)`, 'definition of the DFT'); + +// The random numbers: each pair of fields takes its own consecutive block of the generator's cycle. +// That holds only if the normals use exactly two uniforms per pair, which is checked too. +const s0 = C.stream(7, 0, 4), s1 = C.stream(7, 1, 4), b0 = [...Array(5)].map(() => s0()), b1 = s1(); +let threw = false; try { C.stream(7, 2, 2 ** 31); } catch (e) { threw = true; } +let draws = 0; const counted = C.stream(7, 0, 64), g10 = C.gaussian(() => (draws++, counted())); for (let i = 0; i < 10; i++) g10(); +check('Random streams: block 1 starts where block 0 ends, overlapping blocks are refused, and 10 normals use 10 uniforms', + b0[4] === b1 && threw && draws === 10, `5th draw of block 0 = 1st draw of block 1: ${b0[4] === b1}; overlapping request refused: ${threw}; uniforms for 10 normals: ${draws}`, + 'mulberry32 advances its state by a fixed odd step'); + +// The surfaces, against a naive synthesis: the same normals, the filter |k|^(-beta/2), an O(N^2) inverse DFT. +function naive(L, beta, pair) { + const N = L * L, g = C.gaussian(C.stream(SEED, pair, 2 * N)), Wr = new Float64Array(N), Wi = new Float64Array(N), kof = i => (i <= L / 2 ? i : i - L); + for (let v = 0; v < N; v++) { Wr[v] = g(); Wi[v] = g(); } + for (let j = 0; j < L; j++) for (let i = 0; i < L; i++) { const k = Math.hypot(kof(i), kof(j)), F = k ? Math.pow(k, -beta / 2) : 0; Wr[j * L + i] *= F; Wi[j * L + i] *= F; } + const fr = new Float64Array(N), fi = new Float64Array(N); + for (let y = 0; y < L; y++) for (let x = 0; x < L; x++) { + let sr = 0, si = 0; + for (let j = 0; j < L; j++) for (let i = 0; i < L; i++) { const th = 2 * Math.PI * (i * x + j * y) / L, c = Math.cos(th), s = Math.sin(th); sr += Wr[j * L + i] * c - Wi[j * L + i] * s; si += Wr[j * L + i] * s + Wi[j * L + i] * c; } + fr[y * L + x] = sr; fi[y * L + x] = si; + } + return [C.normalise(fr), C.normalise(fi)]; +} +let synth = 0; +for (const [L, beta, pair] of [[16, 3, 0], [16, 2.4, 3]]) { const a = C.gaussianFields(L, beta, SEED, pair), b = naive(L, beta, pair); for (let k = 0; k < 2; k++) for (let i = 0; i < L * L; i++) synth = Math.max(synth, Math.abs(a[k][i] - b[k][i])); } +const slopes = [2.4, 3, 11 / 3].map(beta => C.gaussianFields(512, beta, SEED, 0).map(f => { + const sp = C.radialSpectrum(f, 512), xs = [], ys = []; sp.k.forEach((k, i) => { if (k >= 2 && k <= 128) { xs.push(Math.log(k)); ys.push(Math.log(sp.P[i])); } }); return [beta, -C.lsq(xs, ys).slope]; +})).flat(); +check('Surfaces equal a naive synthesis (L = 16), and their power spectra recover beta (L = 512, 2 <= |k| <= 128)', synth < 1e-10 && slopes.every(([b, s]) => near(s, b, 0.06)), + `largest difference ${e1(synth)} (under 1e-10); fitted beta ${slopes.map(([b, s]) => `${f3(s)} for ${f3(b)}`).join(', ')} (+- 0.06)`, 'construction'); + +// Marching squares on shapes with known answers. +const L3 = 128, r0 = 20.3, circle = (cx, cy) => { const f = new Float64Array(L3 * L3); for (let y = 0; y < L3; y++) for (let x = 0; x < L3; x++) { const dx = ((x - cx + L3 / 2) % L3 + L3) % L3 - L3 / 2, dy = ((y - cy + L3 / 2) % L3 + L3) % L3 - L3 / 2; f[y * L3 + x] = r0 - Math.hypot(dx, dy); } return f; }; +const mid = C.contours(circle(63.4, 64.7), L3, 0), corner = C.contours(circle(0.4, 127.7), L3, 0); +check('A cone cut at r0 = 20.3: one loop with P = 2 pi r0, A = pi r0^2, R_g = r0; the same circle across both seams of the torus is identical', + mid.P.length === 1 && near(mid.P[0] / (2 * Math.PI * r0), 1, 0.003) && near(mid.A[0] / (Math.PI * r0 * r0), 1, 0.003) && near(mid.Rg[0] / r0, 1, 0.003) && + corner.P.length === 1 && corner.nWrap === 0 && Math.abs(corner.P[0] - mid.P[0]) < 1e-6 && Math.abs(corner.A[0] - mid.A[0]) < 1e-3, + `P/2pi r0 ${(mid.P[0] / (2 * Math.PI * r0)).toFixed(5)}, A/pi r0^2 ${(mid.A[0] / (Math.PI * r0 * r0)).toFixed(5)}, R_g/r0 ${(mid.Rg[0] / r0).toFixed(5)} (1 +- 0.003); across the seams P ${corner.P[0].toFixed(4)} vs ${mid.P[0].toFixed(4)}`, 'a circle'); +const L4 = 64, band = new Float64Array(L4 * L4), diag = new Float64Array(L4 * L4); +for (let y = 0; y < L4; y++) for (let x = 0; x < L4; x++) { band[y * L4 + x] = Math.sin(2 * Math.PI * (y + 0.37) / L4) + 0.2 * Math.sin(2 * Math.PI * x / L4); diag[y * L4 + x] = Math.sin(2 * Math.PI * (x + y + 0.37) / L4); } +const cb = C.contours(band, L4, 0), cd = C.contours(diag, L4, 0); +check('Contours around the torus are found and dropped (bands, and diagonal lines)', cb.nWrap === 2 && cb.P.length === 0 && cd.nWrap === 2 && cd.P.length === 0 && near(cd.wrapLength, 2 * Math.SQRT2 * L4, 1e-6), + `bands ${cb.nWrap} around, ${cb.P.length} closed; diagonals ${cd.nWrap} around, length ${cd.wrapLength.toFixed(4)} (2 sqrt2 L = ${(2 * Math.SQRT2 * L4).toFixed(4)})`, 'construction'); +const sad = (hi, lo, bd) => { const L5 = 8, f = new Float64Array(L5 * L5).fill(lo); if (bd) { f[2 * L5 + 3] = hi; f[3 * L5 + 2] = hi; } else { f[2 * L5 + 2] = hi; f[3 * L5 + 3] = hi; } const c = C.contours(f, L5, 0); return { n: c.P.length, ok: near(C.aboveArea(f, L5, 0), c.A.reduce((a, x) => a + x, 0), 1e-12) }; }; +const sd = [sad(3, -1, false), sad(1, -3, false), sad(3, -1, true), sad(1, -3, true)]; +check('Saddle cells: two diagonal highs are joined when the cell average is above the level, split when below; areas agree cell by cell', + sd[0].n === 1 && sd[1].n === 2 && sd[2].n === 1 && sd[3].n === 2 && sd.every(s => s.ok), `loops ${sd.map(s => s.n).join(', ')} (expected 1, 2, 1, 2)`, 'construction'); + +// Random surfaces: loop areas add up to the area above the level (contours() itself throws if a walk fails to close or +// reuses an edge); the loops do not change under a periodic shift, a quarter turn, or turning the surface upside down. +const rows = []; +for (const beta of [2.4, 3.4]) for (let p = 0; p < 3; p++) for (const f of C.gaussianFields(256, beta, SEED, p)) { + const c = C.contours(f, 256, 0); + if (c.nWrap) continue; // a contour around the torus has no inside, so its area is not defined + const d = (C.aboveArea(f, 256, 0) - c.A.reduce((a, x) => a + x, 0)) / 65536; + rows.push(Math.abs(d - Math.round(d))); +} +check("Green's theorem: on random surfaces without contours around the torus, loop areas add up to the area above the level", + rows.length >= 6 && rows.every(d => d < 1e-9), `${rows.length} of 12 surfaces at L = 256 compared; largest mismatch ${e1(Math.max(...rows))} of L^2, modulo L^2 (under 1e-9)`, "Green's theorem"); +const g = C.gaussianFields(256, 2.4, SEED, 9)[0], key = c => c.P.map((p, i) => [p, Math.abs(c.A[i])]).sort((a, b) => a[0] - b[0]); +const turn = new Float64Array(g.length), shift = new Float64Array(g.length); +for (let y = 0; y < 256; y++) for (let x = 0; x < 256; x++) { turn[x * 256 + (255 - y)] = g[y * 256 + x]; shift[((y + 77) & 255) * 256 + ((x + 31) & 255)] = g[y * 256 + x]; } +const base = key(C.contours(g, 256, 0)), others = [turn, shift, g.map(v => -v)].map(h => key(C.contours(h, 256, 0))); +const diff = Math.max(...others.map(o => o.length === base.length ? Math.max(...o.map((q, i) => Math.max(Math.abs(q[0] - base[i][0]), Math.abs(q[1] - base[i][1])))) : Infinity)); +check('Loops unchanged by a periodic shift, a quarter turn and turning the surface upside down', diff < 1e-8, `${base.length} loops; largest change in P or A ${e1(diff)} (under 1e-8)`, 'symmetry'); + +// 2. Which dimension. Single loops have dimension (3 - H)/2 (Kondev & Henley); all the contours together, 2 - H +// (Mandelbrot). L = 512, 16 surfaces for each of beta = 2.4 and 3. On a finite grid the roughness differs from +// (beta - 2)/2, so H_eff, measured over the same range of sizes, is used for the comparison. +const est = beta => { const recs = []; for (let p = 0; p < 8; p++) for (const f of C.gaussianFields(512, beta, SEED, p)) recs.push(C.fieldRecord(f, 512, { yardMax: 16 })); return { th: C.theory(beta), e: C.estimate(recs, 512, 16) }; }; +const lo = est(2.4), hi = est(3); +check('Loops sit nearer (3 - H)/2 and the whole level set nearer 2 - H, at beta = 2.4 and 3', + [lo, hi].every(({ th, e }) => Math.abs(e.loops.D - th.loop) < Math.abs(e.loops.D - th.set) && Math.abs(e.yard.D - th.set) < Math.abs(e.yard.D - th.loop)), + [lo, hi].map(({ th, e }) => `beta ${th.H * 2 + 2}: loops ${f3(e.loops.D)}, level set ${f3(e.yard.D)} against ${f3(th.loop)} and ${f3(th.set)}`).join('; '), + 'Kondev & Henley, Phys. Rev. Lett. 74, 4580 (1995); Mandelbrot'); +// The loops sit about 0.04 above (3 - H_eff)/2 at larger grids (node turbulence/edges.mjs); on this grid the offset +// is tested against that value with its own jackknife error. +const offs = [lo, hi].map(({ th, e }) => { + const v = e.loops.loo.map((d, i) => d - (3 - e.HeffLoop.loo[i]) / 2), n = v.length, m = v.reduce((a, x) => a + x, 0) / n; + return { beta: th.H * 2 + 2, off: e.loops.D - (3 - e.HeffLoop.D) / 2, se: Math.sqrt((n - 1) / n * v.reduce((a, x) => a + (x - m) ** 2, 0)) }; +}); +check('Each loop dimension sits at (3 - H_eff)/2 plus the 0.04 offset found on larger grids, within 3 jackknife errors', + offs.every(o => Math.abs(o.off - 0.04) < 3 * o.se), + offs.map(o => `beta ${o.beta}: offset ${f3(o.off)} +- ${f3(o.se)}`).join('; ') + ' (this seed; a constant D fails)', + 'Kondev & Henley 1995 (their x_l = 1/2 is a conjecture, exact only at H = 0 and 1); the offset from edges.mjs at L = 1024 to 4096'); +// Both spectra use the same random numbers, field by field, so the fall and the rise are differenced field by field +// inside a paired jackknife. +const pj = fn => { const n = lo.e.loops.loo.length, v = Array.from({ length: n }, (_, i) => fn(i)), m = v.reduce((a, x) => a + x, 0) / n; return Math.sqrt((n - 1) / n * v.reduce((a, x) => a + (x - m) ** 2, 0)); }; +const dDi = i => lo.e.loops.loo[i] - hi.e.loops.loo[i], dHi = i => hi.e.HeffLoop.loo[i] - lo.e.HeffLoop.loo[i]; +const dD = lo.e.loops.D - hi.e.loops.D, dH = hi.e.HeffLoop.D - lo.e.HeffLoop.D; +const seD = pj(dDi), seHalf = pj(i => dDi(i) - dHi(i) / 2), seAll = pj(i => dDi(i) - dHi(i)); +check('From beta 2.4 to 3 the loop dimension falls, by half the rise in H_eff (single loops), not by all of it (the level set)', + dD > 2 * seD && Math.abs(dD - dH / 2) < 2 * seHalf && Math.abs(dD - dH) > 3 * seAll, + `fall ${f3(dD)} +- ${f3(seD)} (paired jackknife, over 2 errors); H_eff rises by ${f3(dH)}: fall minus half of it ${f3(dD - dH / 2)} +- ${f3(seHalf)} (within 2 errors), ` + + `minus all of it ${f3(dD - dH)} +- ${f3(seAll)} (beyond 3; verified on this seed, about 80% of seeds clear it at this size)`, 'Kondev & Henley 1995'); + +console.log(fails ? `\n${fails} check(s) failed` : '\nall checks pass'); +process.exit(fails ? 1 : 0); diff --git a/tests/turbulence.mjs b/tests/turbulence.mjs new file mode 100644 index 0000000..d36c9c3 --- /dev/null +++ b/tests/turbulence.mjs @@ -0,0 +1,119 @@ +// Checks the turbulence model against exact relations, published numbers and its own long run. Run: npm test. +// Each check states the source, the published or expected value, the model's value and the tolerance. +import { createRequire } from 'node:module'; +const require = createRequire(import.meta.url); +const T = require('../turbulence/model.js'); + +let fails = 0; +const fmt = x => (x === 0 ? '0' : x.toExponential(2)); +function within(name, got, want, factor, src) { + const ok = got > 0 && got / want <= factor && want / got <= factor; + console.log(`${ok ? 'PASS' : 'FAIL'} ${name}\n model ${fmt(got)} expected ${fmt(want)} (within x${factor}) ${src}`); + if (!ok) fails++; +} +function close(name, got, want, tol, src) { + const ok = Math.abs(got - want) <= tol; + console.log(`${ok ? 'PASS' : 'FAIL'} ${name}\n model ${got.toFixed(4)} expected ${want.toFixed(4)} (+- ${tol}) ${src}`); + if (!ok) fails++; +} +function check(name, ok, detail, src) { + console.log(`${ok ? 'PASS' : 'FAIL'} ${name}\n ${detail} ${src}`); + if (!ok) fails++; +} +const SAB = "Sabra model of L'vov et al., Phys. Rev. E 58, 1811 (1998)"; + +// 1. The shell model conserves what it cascades: with a + b + c = 0, energy and a helicity-like second invariant +// (inviscid, unforced, 10 shells, t = 200, dt = 1e-3). The opposite sign convention for c does not conserve energy. +const good = T.conservation({ N: 10, dt: 1e-3, T: 200, qs: [-2] }); +check('Energy and helicity conserved to integrator accuracy when a + b + c = 0 (a = 1, b = c = -1/2)', + Math.abs(good.dE) < 1e-9 && Math.abs(good.dH[0]) < 1e-8, `max |dE|/E = ${fmt(good.dE)}, max |dH|/norm = ${fmt(good.dH[0])}`, + SAB + ', Secs. II A and III A (as read from the paper in #34\'s review); derivation in model.js'); +const bad = T.conservation({ N: 10, c: 0.5, dt: 1e-3, T: 200, stopAt: 0.5 }); +check('With a + b - c = 0 instead (c = +1/2), energy is not conserved', bad.tStop !== null && bad.tStop < 5, + `E changes by more than 50% by t = ${bad.tStop === null ? 'never' : bad.tStop.toFixed(2)}`, 'model'); + +// 2. The cascade: in a steady state the same energy flux passes through every shell between the forcing and the +// dissipation, and a third-order correlator scales exactly as 1/k_n (the shell model's 4/5 law). +// CI run: 18 shells, nu = 1e-6, force 0.005(1 + i) on shell 0, 4,800 time units (about 860 large-eddy times). +const res = T.simulate(T.RUNS.CI), a = T.analyse(res), av = a.av, eps = a.eps, k = res.model.k, P = a.plateau; +console.log(` (CI run: ${res.steps} steps, eps = ${fmt(eps)}, Re = U L/nu = ${fmt(a.Re)}, k_eta = (eps/nu^3)^(1/4) at shell ${a.nEta.toFixed(2)})`); +close('Mean dissipation over mean injection', a.dis / eps, 1, 0.01, 'model (energy balance)'); +const shape = [0, 1, 2, 3, 4, 5, 6].map(n => av.Pi[n] / P); +check('The energy flux is the same through every shell from the forcing to n = 6, and equals the injection', + Math.abs(P / eps - 1) <= 0.012 && shape.every(x => Math.abs(x - 1) <= 0.006), + `plateau P/eps = ${(P / eps).toFixed(4)} (+- 0.012); Pi_n/P = ${shape.map(x => x.toFixed(4)).join(' ')} (+- 0.006)`, 'exact flux from the triads (model.js)'); +const law = [1, 2, 3, 4, 5].map(n => [3 * k[n + 2] * av.Sig[n] / P, 1 - Math.pow(-0.5, n)]); +check("Shell-model 4/5 law: 3 k_n = eps (1 - (-1/2)^n), helicity injected at eps on shell 0", + law.every(([g, w]) => Math.abs(g - w) <= 0.01), law.map(([g, w]) => `${g.toFixed(4)}/${w.toFixed(4)}`).join(' ') + ' (model/exact, +- 0.01)', + SAB + ', eq. (16) with Sigma_n = -S_3 (as read from the paper in #34\'s review); derived in model.js'); + +// 3. But the statistics do not repeat exactly. Kolmogorov 1941 (K41) predicts structure-function exponents p/3; the +// cascade is intermittent, so zeta_2 > 2/3 and zeta_6 < 2. Exponents from extended self-similarity of the +// flux-based structure functions <|Pi_n/k_n|^(p/3)> over shells 3..8 (the plain <|u_n|^p> ripple from shell to shell). +const z = p => a.zeta.essF[p - 1], ze = p => a.err.essF[p - 1]; +check('Anomalous scaling: zeta_2 above 2/3 and zeta_6 below 2, each by more than 5 batch-means errors', + z(2) - 2 / 3 > 5 * ze(2) && 2 - z(6) > 5 * ze(6), `zeta_2 = ${z(2).toFixed(3)} +- ${ze(2).toFixed(4)}, zeta_6 = ${z(6).toFixed(3)} +- ${ze(6).toFixed(4)}`, + 'K41: Kolmogorov, Dokl. Akad. Nauk SSSR 30 (1941)'); +const LONG = { 2: [0.719, 0.01], 4: [1.255, 0.015], 6: [1.716, 0.03] }; +for (const p of [2, 4, 6]) close(`zeta_${p}, CI run against the long run`, z(p), LONG[p][0], LONG[p][1], 'this model, node turbulence/table.mjs (20,000 time units)'); +const LVOV = { 2: [0.720, 0.008], 4: [1.256, 0.012], 6: [1.706, 0.015] }; +check("Published Sabra exponents: zeta_2, zeta_4, zeta_6 within twice L'vov et al.'s errors", + [2, 4, 6].every(p => Math.abs(z(p) - LVOV[p][0]) <= 2 * LVOV[p][1]), + [2, 4, 6].map(p => `zeta_${p} ${z(p).toFixed(3)} against ${LVOV[p][0]} +- ${LVOV[p][1]}`).join('; ') + ' (theirs: 34 shells, forcing on shells 0-1, <|u|^q>)', + SAB + ', Table I (as read from the paper in #34\'s review)'); +const KUCHLER = [[0.693, 0.011], [0.698, 0.011]]; +check('The shell model is more intermittent than measured turbulence: its zeta_2 sits over 1.5 errors above both published values', + KUCHLER.every(([v, e]) => z(2) - v > 1.5 * e), + `model ${z(2).toFixed(3)}; measured ${KUCHLER.map(([v, e]) => `${v} +- ${e} (${((z(2) - v) / e).toFixed(1)} errors)`).join(' and ')} at R_lambda 2000-6000; She-Leveque ${T.sheLeveque(2).toFixed(3)}`, + 'Kuchler, Bodenschatz & Bewley, arXiv:2006.10993v2 (2020), Fig. 3 (the abstract still says +- 0.003), and arXiv:2106.01795v2, Fig. 5, as read in #34\'s review; She & Leveque, PRL 72, 336 (1994)'); + +// 4. Where the cascade stops: K41 puts the end at k_eta = (eps/nu^3)^(1/4) up to an O(1) factor it does not predict. +check('The flux falls to half, and <|u_n|^2> to half its power law, within 1.5 shells (a factor 2.8) of k_eta', + Math.abs(a.nd - a.nEta) <= 1.5 && Math.abs(a.ndSpec - a.nEta) <= 1.5, + `flux half-point ${a.nd.toFixed(2)}, spectrum half-point ${a.ndSpec.toFixed(2)}, k_eta at ${a.nEta.toFixed(2)} (shells)`, 'K41 dissipation scale'); + +// 5. Real flows: the Kolmogorov scale eta = (nu^3/eps)^(1/4) for published dissipation rates. +const etaMm = (nu, e) => T.kolmogorov(nu, e).eta * 1e3; +const ocean = [1e-8, 1e-7].map(e => etaMm(T.NU.seawater, e)), air = [0.4, 0.05, 1e-7].map(e => etaMm(T.NU.air, e)); +check('Upper-ocean mixed layer: eps 1e-8 to 1e-7 W/kg, as at one station, puts eta at 1.8 to 3.3 mm; 1 mm needs 1e-6', + ocean[0] > 3 && ocean[0] < 3.5 && ocean[1] > 1.7 && ocean[1] < 2, `eta = ${ocean.map(x => x.toFixed(2)).join(' and ')} mm; 1e-6 W/kg gives ${etaMm(T.NU.seawater, 1e-6).toFixed(2)} mm`, + 'arithmetic; eps from Yang et al., Sensors 22, 1167 (2022), one 25-hour station, as read in #34\'s review'); +check('Air: eta from 0.3 mm (eps 0.4 m^2/s^3, a cold front) to 1.4 cm (1e-7, a quiet night)', + air[0] > 0.28 && air[0] < 0.32 && air[2] > 13 && air[2] < 14.5, `eta = ${air.map(x => x.toFixed(2)).join(', ')} mm for eps = 0.4, 0.05, 1e-7`, + 'arithmetic; eps from Balsley, Tjernstrom & Svensson, AMS 18th Symp. Boundary Layers and Turbulence (2008), as seen in search text'); + +// 6. The solar wind is too dilute for viscosity: protons cross 0.2 to 1 AU between collisions, so the MHD cascade ends at +// the ion scales instead, and a steeper kinetic one continues to the electrons. Typical wind at 1 AU: n = 5 cm^-3, +// V = 400 km/s, B = 5 nT, T = 1e5 K (slow-wind protons are cooler, a few 1e4 K). +const W = T.WIND, di = T.inertialLength(W.n), rho = T.gyroradius(W.T, W.B, 2), col = T.protonCollisions(W), cool = T.protonCollisions({ ...W, T: 4e4 }); +within('Proton inertial length c/omega_pi at n = 5 cm^-3 (km)', di / 1e3, 228 / Math.sqrt(5), 1.002, 'NRL Plasma Formulary, 228 km / sqrt(n per cm^3), as seen in search text'); +check('Collisionless: the proton mean free path is over 10^5 inertial lengths, at 1e5 K and at 4e4 K', col.mfp / di > 1e5 && cool.mfp / di > 1e5, + `mean free path ${(col.mfp / T.C.AU).toFixed(2)} AU at 1e5 K (collision time ${(col.tau / 86400).toFixed(0)} days), ${(cool.mfp / T.C.AU).toFixed(2)} AU at 4e4 K; ${fmt(cool.mfp / di)} d_i at least (ln Lambda = 25)`, + 'arithmetic; Braginskii, Rev. Plasma Phys. 1, 205 (1965), ion collision time'); +const psi = T.parkerAngle(W.V), fdi = T.taylorFrequency(W.V, di, psi), frho = T.taylorFrequency(W.V, rho, psi); +const fdi0 = T.taylorFrequency(W.V, di); +check('The ion inertial length swept past the spacecraft, with and without the Parker-spiral angle, lands in the observed 0.1-1 Hz', + [fdi, fdi0].every(f => f >= 0.1 && f <= 1), + `V sin(${(psi * 180 / Math.PI).toFixed(0)} deg)/(2 pi d_i) = ${fdi.toFixed(3)} Hz, V/(2 pi d_i) = ${fdi0.toFixed(3)} Hz; gyroradius ${frho.toFixed(2)} Hz. The spread does not choose between them`, + 'Chen et al., GRL 41, 8081 (2014), as confirmed in #34\'s review; Sahraoui et al. 2009 (0.4 Hz, one interval, a proton-gyroscale break) and Leamon et al. 1998 for context'); + +// 7. Molecular clouds: the cascade is supersonic at large scales, and the supersonic range ends at the sonic scale, +// where the non-thermal velocity falls to the sound speed (0.19 km/s at 10 K); a shallower subsonic cascade +// continues below. Larson's data reach about 0.1 pc; Solomon et al.'s giant clouds are pc-sized and larger, so their +// sonic size is an extrapolation of over a decade. Each relation is read with the thermal width its own dispersion +// holds: the gas's sigma_s = 0.32 km/s (3D) in Larson's, CO's thermal width in Solomon et al.'s. +const sL = T.sonicScale(T.LARSON81), sS = 2 * T.sonicScale(T.SOLOMON87), mL = T.machAt(T.LARSON81, 0.1); +check('Sonic sizes from Larson (maximum size) and Solomon et al. (diameter) within a factor of 2 of the 0.10 pc filament width', + [sL, sS].every(x => x / 0.1 < 2 && 0.1 / x < 2), `Larson ${sL.toFixed(3)} pc, Solomon et al. ${sS.toFixed(3)} pc; filaments 0.10 +- 0.03 pc (FWHM)`, + 'Larson, MNRAS 194, 809 (1981), pp. 813-814; Solomon et al., ApJ 319, 730 (1987), p. 737; Arzoumanian et al., A&A 529, L6 (2011); all as read in #34\'s review'); +check("At Larson's smallest scales his relation is about sonic: non-thermal Mach number 0.8 to 1.25 at 0.1 pc", + mL > 0.8 && mL < 1.25, `Mach ${mL.toFixed(2)} at 0.1 pc (${T.machAt(T.LARSON81, 100).toFixed(0)} at 100 pc)`, + 'this study\'s tolerance on Larson (1981), eq. (1)'); +const objs = T.LARSON_OBJECTS.map(o => ({ ...o, mach: T.machOf(T.LARSON81, o.sigma) })); +check("Larson's small objects: TMC1 and TMC2 subsonic, as he says, and the others supersonic, in his convention", + objs.every(o => (o.mach < 1) === o.subsonic), + objs.map(o => `${o.name} (${o.size} pc, ${o.sigma} km/s) Mach ${o.mach.toFixed(2)}`).join('; '), + 'Larson (1981), Table 1 (pp. 812-813) and p. 820, as read in #34\'s re-review'); + +console.log(fails ? `\n${fails} check(s) failed` : '\nall checks pass'); +process.exit(fails ? 1 : 0); diff --git a/turbulence/README.md b/turbulence/README.md new file mode 100644 index 0000000..005af35 --- /dev/null +++ b/turbulence/README.md @@ -0,0 +1,487 @@ +# Turbulence: a pattern that repeats over many decades, and the sizes where it stops + +Status: **exploration, not peer reviewed, no new result claimed.** Issue +[#19](https://github.com/Foundation-0/fractals/issues/19). + +A turbulent flow looks statistically alike at every size between the scale where energy is stirred in and the +scale where it turns into heat. This study asks how exact that repetition is, using a cascade model that runs in +a few seconds. It then asks where the repetition ends in air, water, the solar wind and molecular clouds, and +whether the fractal edges of clouds really match laboratory turbulence. + +## The cascade + +Kolmogorov's 1941 picture (K41): energy fed in at a large scale L passes to smaller and smaller eddies at a rate +ε per unit mass. The rate is the same at every scale in between. Viscosity ν takes the energy only at + +``` +η = (ν³/ε)^(1/4) +``` + +so the range of repeating scales is L/η ~ Re^(3/4), with Re the Reynolds number. In between, the velocity +difference across a distance r depends only on ε and r. Its p-th moment then grows as r^(p/3). One result is +exact: the third moment of the longitudinal velocity difference is −(4/5) ε r in the inertial range (Kolmogorov's +4/5 law). The energy spectrum is C ε^(2/3) k^(−5/3). The constant is about 1.5 for the 3D spectrum. For the 1D +longitudinal spectrum, Sreenivasan's (1995) compilation of experiments gives 0.53 with a scatter of 0.055, which +isotropy turns into 1.62 ± 0.17 in 3D, so "about 1.5" is fair once its spectrum is named. + +## A cascade to run: the Sabra shell model + +The Sabra model (L'vov, Podivilov, Pomyalov, Procaccia and Vandembroucq 1998) keeps one complex velocity u_n per +octave of wavenumber, k_n = 2^n, coupled to its nearest neighbours in scale: + +``` +du_n/dt = i (a k_{n+1} u_{n+2} u*_{n+1} + b k_n u_{n+1} u*_{n−1} − c k_{n−1} u_{n−1} u_{n−2}) − ν k_n² u_n + f_n +``` + +The nonlinear terms move energy between shells without creating it if a + b + c = 0. A second quantity, Σ (a/c)^n +|u_n|², is then conserved too; with a = 1 and b = c = −1/2 it plays the part of helicity, as in three dimensions +(shell models are reviewed by Biferale 2003). The run here forces shell 0 and lets viscosity act on all shells: 18 +shells and about 860 large-eddy times in the test, 21 shells and about 3,600 in the table. Both are seeded and +bit-for-bit reproducible. + +**What repeats exactly: the flux.** In a steady state the same energy flux passes through every shell from the +forcing down to the dissipation range. In the long run it stays within 0.2% of ε from shell 0 to shell 8. A +third-order correlation obeys the model's own version of the 4/5 law: + +``` +3 k_n ⟨Im u*_{n−1} u*_n u_{n+1}⟩ = ε − δ (−1/2)^n +``` + +Here δ is the helicity flux, equal to ε when only shell 0 is forced. This is derived in `model.js` from the +equation above; it is L'vov et al.'s own eq. (16), with Σ_n = −S₃. The model matches it to better than 0.1% in +shells 1 to 8 of the long run. Over 32 seeds of the test's shorter run the largest deviation in shells 1 to 5 is +typically 0.0014, at most 0.006. + +**What does not repeat exactly: the exponents.** If each shell were a scaled copy of the last, the p-th moment +would fall as k_n^(−p/3). It does not: + +| p | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | +|---|---|---|---|---|---|---|---|---| +| Sabra (this model) | 0.395 | 0.719 | 1 | 1.255 | 1.493 | 1.716 | 1.93 | 2.12 | +| K41, p/3 | 0.333 | 0.667 | 1 | 1.333 | 1.667 | 2 | 2.333 | 2.667 | +| She–Leveque | 0.364 | 0.696 | 1 | 1.280 | 1.538 | 1.778 | 2.001 | 2.211 | +| L'vov et al. 1998, moments of u | | 0.720 | | 1.256 | | 1.706 | | | +| L'vov et al. 1998, flux-based | | 0.719 | | 1.249 | 1.477 | 1.691 | 1.893 | | + +The model's exponents come from the flux-based moments ⟨|Π_n/k_n|^(p/3)⟩, scaled against the third (extended +self-similarity, Benzi et al. 1993), over shells 3 to 11. L'vov et al.'s Table I comes from 34 shells forced at +random on shells 0 and 1, without helicity. This model agrees with their moments of u within their errors (± 0.008, +0.012 and 0.015 for p = 2, 4 and 6) and sits 0.006 to 0.034 above their flux-based column for p = 4 to 7; they also +give ζ₂ = 0.728 ± 0.002 by another method. Statistical errors are 0.0002 for p = 2 and 0.004 for p = 6; the fit +range moves them by up to 0.002 and 0.006. The plain moments ⟨|u_n|^p⟩ ripple from shell to shell. Their fitted +exponents shift by several hundredths, up to about 0.1 for p = 6, with the fit range, so they are printed but not +used. + +At small scales the velocity is more often near zero, and more often very large, than a scaled copy would make it. +The tails of its distribution fatten: the flatness ⟨|u|⁴⟩/⟨|u|²⟩² grows toward small scales as k^(2ζ₂ − ζ₄) = k^0.18 +here, from 1.8 at shell 2 to 5.7 at shell 10. That is intermittency, and it is the main way real turbulence breaks +exact self-similarity. The shell model has more of it than measured flows. Its ζ₂ is 0.720 in the test's run (0.719 +in the long run). A wind tunnel at Taylor-scale Reynolds numbers of 2,000 to 6,000 gives ζ₂ = 0.693 ± 0.011 after +correcting for the decay of the turbulence (Küchler, Bodenschatz and Bewley 2020, Fig. 3 of the revised preprint; +its abstract still gives ± 0.003), close to She and Leveque's 0.696, and the same group's later analysis 0.698 +± 0.011 (Fig. 5 of arXiv:2106.01795v2). The model sits about two of the measurement's errors above them (2.0 to 2.4 +for the test's run). That group's published letter reads their high-Reynolds-number statistics as a universal form +rather than a pure power law (Küchler, Bewley and Bodenschatz 2023). + +**Where it stops.** The flux falls to half about one shell before k_η = (ε/ν³)^(1/4), and the spectrum falls to +half its power law within half a shell of k_η. Over four viscosities the cascade widens with the Reynolds number: + +| ν | Re = UL/ν | Shell of k_η | Flux half-point | Spectrum half-point | Width, decades | (3/4) log₁₀ Re | +|---|---|---|---|---|---|---| +| 10⁻⁵ | 1.7×10⁴ | 10.15 | 9.04 | 9.75 | 2.72 | 3.18 | +| 10⁻⁶ | 1.8×10⁵ | 12.65 | 11.60 | 12.55 | 3.49 | 3.94 | +| 10⁻⁷ | 1.8×10⁶ | 15.15 | 14.30 | 15.15 | 4.31 | 4.69 | +| 10⁻⁸ | 1.8×10⁷ | 17.64 | 16.95 | 17.64 | 5.10 | 5.44 | + +The width grows by 0.78 to 0.79 decades per decade of Re. The cut can be placed where the flux has fallen to a +quarter, a half or three quarters of its plateau, or where the spectrum has fallen to a half or a quarter of its +power law (at three quarters the spectrum cut falls near the end of the fit range at ν = 10⁻⁵ and is left out). +Over eight seeds and these five cuts the slopes run from 0.770 to 0.794; the flux half-point gives 0.787 with a +seed-to-seed SD of 0.003, so an error of 0.001 on the mean (`node turbulence/seeds.mjs`). K41 gives 3/4: the flux +half-point drifts toward k_η by about half a shell over three decades. Intermittency accounts for the difference. +In the multifractal picture the cut-off goes as Re^(1/(1 + h*)), with h* = dζ_p/dp at p = 3, the scaling of the +velocity that carries the flux (Paladin and Vulpiani 1987; Frisch 1995; in this form, Boffetta, Mazzino and +Vulpiani 2008). With the model's exponents h* ≈ (ζ₄ − ζ₂)/2 = 0.27, which gives Re^0.79. Viscous loss ν k² S₂(k) +equal to ε gives Re^(1/(2 − ζ₂)) = Re^0.78. On seed noise alone the flux half-point favours the multifractal value, +but the two stay unseparated: the cut's definition moves the eight-seed means from 0.781 to 0.791, and the +predictions rest on estimates of h* good to about 0.01. Boffetta et al.'s eq. (II.5) ties the spectrum's cut to +h*(2) ≈ (ζ₃ − ζ₁)/2 = 0.30, Re^0.77; the measured spectrum cuts (0.781 and 0.786) sit above that. A local Reynolds +number of 1 on the rms velocity, S₂(k)^(1/2) = ν k, gives Re^(2/(2 + ζ₂)) = Re^0.74, below K41 and below every +measured slope. + +## Real cascades, and what stops them + +| Flow | Repeats down to | Decades | What ends it | +|---|---|---|---| +| Wind tunnel, R_λ ≈ 6,000 | η | about 5 | viscosity | +| Atmospheric surface layer, R_λ 1,600 to 6,600 (Gulitski et al. 2007, Part 1) | η = 0.3 to 0.5 mm in wind | 4.2 to 5.1 | viscosity | +| Upper ocean | η = 1 to 3 mm, more in quiet water | | viscosity | +| Solar wind | the ion scales, about 100 km (the MHD range) | | the ions' own motion, not collisions | +| Interstellar electron density | at least 2×10⁶ m near the Sun | 6.7 measured | not seen locally; 50 to 200 km inferred on strongly scattered lines of sight | +| Molecular clouds, supersonic range | the sonic scale, about 0.1 pc | | the motions become subsonic; a shallower subsonic cascade continues | + +The decades for the tunnel and the surface layer are L/η = (0.15 R_λ²)^(3/4), with Pope's conventions for L. The +scaling range proper is narrower: dissipation starts near 60 η, and the largest eddies are not yet self-similar. + +### Air and water + +The Kolmogorov scale is often quoted as about a millimetre. It moves with ε as ε^(−1/4): + +| Medium | ε (W/kg) | η | Where | +|---|---|---|---| +| Air | 0.4 | 0.30 mm | surface layer, cold front | +| Air | 0.05 | 0.51 mm | the same, before the front | +| Air | 10⁻⁷ | 1.4 cm | quiet nocturnal boundary layer | +| Seawater | 10⁻⁷ | 1.8 mm | upper mixed layer, by day | +| Seawater | 10⁻⁸ | 3.3 mm | upper mixed layer, quieter | + +The ε values are from Balsley, Tjernström and Svensson (2008) for air and Yang et al. (2022), one 25-hour station, +for the ocean. A millimetre needs ε ≈ 3×10⁻³ W/kg in air and 10⁻⁶ in the sea, which an actively stirred mixed layer +reaches. So "about a millimetre" is the right order in both, with a spread of a factor of 10 either way. + +### Atmosphere: two power laws + +Aircraft spectra of wind and temperature (Nastrom and Gage 1985, over 6,000 commercial flights near the +tropopause) follow k^(−5/3) from a few kilometres to about 300 to 400 km, and k^(−3) at larger scales. The +transition is quoted between about 400 and 1,000 km, most often near 600 km. Only the slopes are universal; the +amplitudes vary with latitude, season and altitude. A two-dimensional cascade forced at one scale would put the +two ranges the other way round, k^(−3) on the small-scale side (as Kraichnan 1967 conjectured), which is why the +mesoscale −5/3 is still argued over. A 2025 preprint (DeWitt and Garrett) reads global sonde data as one +anisotropic cascade with no transition. + +### Solar wind: cut by the ions, not by viscosity + +The solar wind is too dilute for viscosity to end its cascade. A proton at 1 AU travels 0.2 to 1 AU between +Coulomb collisions (Braginskii's collision time: 49 days at 10⁵ K, 12 days at the 4×10⁴ K of cooler, slow wind). +That is 10⁵ to 10⁶ times the proton inertial length c/ω_pi. The MHD cascade instead ends where the protons' own +motion decouples them from the field, and a steeper, kinetic cascade carries on to the electron scales, about a +kilometre: + +| Typical wind: n = 5 cm⁻³, V = 400 km/s, B = 5 nT, T = 10⁵ K (β_p = 0.69) | Length | Swept past: no angle | at the Parker angle, 47° | +|---|---|---|---| +| Proton inertial length c/ω_pi | 102 km | 0.63 Hz | 0.46 Hz | +| Proton gyroradius, v_th = √(2kT/m) | 85 km | 0.75 Hz | 0.55 Hz | +| Proton gyroradius, v_th = √(kT/m) | 60 km | 1.06 Hz | 0.78 Hz | + +A spacecraft sees structures carried past at the wind speed (Taylor's hypothesis). The fluctuations are elongated +along the field, so a size ℓ across the field shows up at V sin θ_BV/(2πℓ), with θ_BV the angle between the field +and the flow; at 1 AU the Parker spiral puts it near 47° for a 400 km/s wind. Observed breaks fall between 0.1 and +1 Hz across conditions, at the larger of the two ion lengths (Chen et al. 2014). Sahraoui et al. (2009) find 0.4 Hz +in one fast, high-β interval and read it as the proton gyroscale. Leamon et al.'s 33 intervals, as replotted by +Leamon et al. (2000), run from about 0.12 to 0.87 Hz with a mean near 0.38 Hz; a fit through the origin puts the +break at 0.686 V sin Ψ/(2π d_i), an unconstrained one at 0.15 Hz + 0.45 V sin Ψ/(2π d_i), and their fit is best for +an orientation near 65°. The inertial length gives 0.46 Hz at the Parker angle and 0.63 Hz without it; both lie +inside the observed spread, which does not choose between them. At β near 1 the two ion lengths are nearly equal, +so the break cannot tell them apart either. + +At larger scales (lower frequencies) than the break the magnetic spectrum falls as k^(−5/3) and the velocity +spectrum as k^(−3/2) (Podesta, Roberts and Goldstein 2007). At smaller scales the kinetic range's slope is not +universal: −2.8 (Alexandrova et al. 2009) and −2.3 (Sahraoui et al. 2009). Smith et al. (2006) find −1 to −4, with +means of −2.0 in magnetic clouds and −2.6 on open field lines; Alexandrova et al. and Bruno and Carbone (2013) +summarise the range as about −2 to −4, and Bruno and Carbone give about −7/3 as typical. + +### The interstellar medium + +Interstellar electron-density fluctuations follow one Kolmogorov-like power law from about 2×10⁶ m to 10¹³ m, +with a three-dimensional spectral index near 3.7 against Kolmogorov's 11/3 (Armstrong, Rickett and Spangler 1995). +That is 6.7 decades. No inner scale is seen near the Sun, where Armstrong et al. bound it below about 10⁸ m; +on strongly scattered lines of sight Spangler and Gwinn (1990) infer 50 to 200 km. Chepurnov and Lazarian (2010) +extend the power law to 10¹⁷ m, about 3 pc, but that extension assumes the Kolmogorov form rather than fitting it. + +### Molecular clouds: the sonic scale + +Cloud turbulence is supersonic on large scales. Its velocity dispersion falls with size, and where the non-thermal +part reaches the sound speed, 0.19 km/s at 10 K, the supersonic range ends. A shallower, subsonic cascade carries on +below it (Federrath et al. 2021). Two classic fits give: + +| Fit | Dispersion | Size | Sonic size | +|---|---|---|---| +| Larson (1981): σ = 1.10 L^0.38 km/s | 3D, all internal motions, thermal ones included | maximum projected size L | 0.099 pc | +| Solomon et al. (1987): σ = 1.0 S^0.5 km/s, 0.72 R^0.5 in this study's conversion | 1D, from ¹²CO line widths | effective radius R = 1.92 S | 0.074 pc radius, 0.147 pc across | + +Each dispersion already holds a thermal part, and each is read in its own convention. Larson's σ is the total 3D rms +velocity of all internal motions, including the gas's thermal σ_s = 0.32 km/s at 10 K (his pp. 813–814), so the +non-thermal 1D dispersion is √((σ² − σ_s²)/3). Solomon et al.'s is the width of the ¹²CO line, which holds only +CO's own thermal width, 0.054 km/s. Larson's own sonic size, 0.04 pc (his p. 820), uses another criterion, σ = σ_s; +this relation gives 0.039 pc for it. Within its errors, (1.0 ± 0.1) S^(0.5 ± 0.05), Solomon et al.'s fit puts the +sonic size between 0.08 and 0.24 pc across. + +Larson's data reach about 0.1 pc, so his sonic size is a short extrapolation, and there his relation is about +sonic (non-thermal Mach 1.0). He calls TMC1 and TMC2 subsonic, the only two in his Table 1: TMC2 (0.1 pc) is his +smallest plotted object, at Mach 0.44 in his convention, and TMC1 (0.2 pc) is at 0.74, while L1407/L1257 (0.12 pc) +and L1544/L63 (0.18 pc) are supersonic, at Mach 2.45 and 1.74. Solomon et al.'s clouds are giant molecular clouds, a +parsec and larger, so theirs is an extrapolation of more than a decade. Both sit near the width of interstellar +filaments, 0.10 ± 0.03 pc in Arzoumanian et al. (2011) and a median of 0.10 pc with an interquartile range of +0.07 pc over 599 filaments in Arzoumanian et al. (2019), and the size of coherent cores, +inside which the line width stops growing with size (Goodman et al. 1998). Pineda et al. (2010) see the dispersion +rise by a factor of 2 within 0.04 pc at a core's edge. In simulations the transition is smooth over a factor of +about 3 in scale (Federrath et al. 2021). + +The fits differ in slope: Larson's 0.38 against Solomon et al.'s 0.5, with Heyer and Brunt (2004) at +0.49 ± 0.15 cloud by cloud. In a cascade the dispersion grows as ℓ^(ζ₂/2). That gives 1/3 for K41, 0.36 for the +shell model and 0.35 for She–Leveque, all incompressible, and 1/2 for a field of shocks. Federrath et al.'s +simulations give 0.49 ± 0.01 in the supersonic range and 0.39 ± 0.02 in the subsonic one. Larson's 0.38 is close +to the subsonic value, but his clouds are supersonic over most of his range (non-thermal Mach 1.0 at 0.1 pc to 19 at +100 pc). Larson read his slope as supersonic turbulence, or subsonic turbulence in warmer gas; Solomon et al. read +0.5 as virial equilibrium. That the difference lies in Larson's sample and size definition, and that the supersonic +data favour about 0.5, is this study's inference. + +### Superfluid helium + +In superfluid helium, vorticity comes in quantised lines with circulation h/m₄ = 9.97×10⁻⁸ m²/s. Above the +spacing between lines, bundles of them carry a classical Kolmogorov cascade. Maurer and Tabeling (1998) measured +k^(−5/3) spectra both above and below the superfluid transition, using a pressure probe as a proxy for velocity. +Below that spacing the energy travels along single lines as Kelvin waves (Barenghi, Skrbek and Sreenivasan +2014). There the repetition stops at a length set by quantum mechanics. + +## Cloud edges: is the agreement with the laboratory real? + +Cloud edges in projected CO maps of the Taurus field follow a perimeter-area relation P ∝ A^(D/2) with +D = 1.36 ± 0.02, the same from about 100 pc down to 0.1 pc (Falgarone, Phillips and Walker 1991). The error is a +combined-fit error. Their maps one by one give 1.46 ± 0.18, 1.32 ± 0.14 and 1.38 ± 0.08 (their Fig. 6a), which agree +with each other, as they say, though they add that error bars are difficult to estimate. Other data for the Taurus +complex gave 1.51 (Stutzki 1993, Fig. 5, on maps compiled by Zimmermann et al.). Rain and cloud areas in the +atmosphere give 1.35 ± 0.05 (Lovejoy 1982). Laser-sheet sections through the turbulent/non-turbulent interfaces of +laboratory jets and wakes give box-counting dimensions of 1.3 to 1.4, sections of surfaces of dimension 2.3 to 2.4 +(Sreenivasan and Meneveau 1986), which Sreenivasan, Ramshankar and Meneveau (1989) explain as 7/3; sections of +iso-concentration and iso-dissipation surfaces give 1.36 ± 0.05 as well (Sreenivasan and Meneveau 1986 and +Meneveau 1989, as Falgarone et al. quote them). The match between clouds and laboratory is often cited. + +A fractional-Brownian estimate depends on which dimension the relation measures. For a Gaussian random surface +with power spectrum k^(−β), so H = (β − 2)/2: + +- all the contours at one level together have dimension 2 − H (Mandelbrot); +- a single contour loop has dimension (3 − H)/2 (Kondev and Henley 1995; a conjecture, x_l = 1/2, exact only at + H = 0 and 1). + +A perimeter-area relation over separate clouds measures single loops, if holes are not counted. `contours.js` checks +this by making the surfaces and measuring both. The table below uses 32 surfaces of 2048 × 2048 points per spectrum +(`node turbulence/edges.mjs`): + +| β | H | (3 − H)/2 | 2 − H | Loops, P ∝ R_g^D | Loops, P ∝ A^(D/2) | All contours together | +|---|---|---|---|---|---|---| +| 2.4 | 0.20 | 1.400 | 1.800 | 1.429 ± 0.012 | 1.430 ± 0.013 | 1.716 | +| 8/3 | 0.33 | 1.333 | 1.667 | 1.357 ± 0.017 | 1.364 ± 0.022 | 1.623 | +| 2.8 | 0.40 | 1.300 | 1.600 | 1.335 ± 0.013 | 1.346 ± 0.017 | 1.570 | +| 3 | 0.50 | 1.250 | 1.500 | 1.275 ± 0.025 | 1.294 ± 0.025 | 1.485 | +| 3.4 | 0.70 | 1.150 | 1.300 | 1.219 ± 0.040 | 1.215 ± 0.039 | 1.309 | +| 11/3 | 0.83 | 1.083 | 1.167 | 1.150 ± 0.025 | 1.145 ± 0.020 | 1.209 | + +R_g is a loop's radius of gyration, fitted over 8 to 128 grid spacings. The errors combine a jackknife over +surfaces with the spread over fit windows. The last column comes from the total contour length measured on +coarser and coarser grids. With holes counted (outline and holes against net area) the loops give 1.65 at β = 2.4 +and 1.52 at 2.8. + +- **Loops follow Kondev and Henley.** On a finite periodic grid the roughness the loops see differs from + (β − 2)/2: higher for flat spectra (0.24 against 0.20 at β = 2.4), lower for steep ones (0.76 against 0.83 at + 11/3). Measured over the same range of sizes, the loop dimension falls with this effective roughness H_eff at a + slope of −0.53 ± 0.05. The six spectra share their random numbers, so the error combines a joint jackknife + (± 0.040) with the spread over fit windows (± 0.031). Single loops predict −1/2; the whole level set −1, which is + 9 errors away. At 4096 points per side the slope is −0.532 ± 0.022, 21 errors from −1; at 1024 it is + −0.47 ± 0.11, 5 errors from −1. +- **An offset.** Loops sit above (3 − H_eff)/2 at every spectrum: by 0.039, 0.038 and 0.035 on average at 1,024, + 2,048 and 4,096 points per side (`node turbulence/edges.mjs 1024`, `4096`), with H_eff measured over the loops' + own range of sizes. Measured over separations two and four times larger it is 0.029 and 0.017 (at 2048), and #34's + review found about 0.014 with a covering-box size in place of R_g. So its size depends on the set-up, it barely + shrinks with the grid, and its origin is open. Kondev, Henley and Salinas's own fits (their Table I, as read in + #34's review) overshoot the formula only at H = 0.6 and above, the smoother surfaces, so they do not explain it. +- **The whole level set follows Mandelbrot:** its dimension tracks 2 − H_eff to within 0.02. +- **What does not change the answer:** + - fitting perimeter against size or against area; + - fitting only the outlines of clouds (loops with the cloud inside), which agree with all loops to within 0.02. +- **Not tested: the level.** The fields have zero mean and their median sits within 0.2 standard deviations of + it, so the median column (which gives the same loops) does not probe the level; contours away from the mean + were not computed. +- **Counting holes changes it.** Adding a cloud's holes to its perimeter and taking its net area gives a larger + number. Falgarone et al. take the area as the number of map positions inside the contour and the perimeter as + the number of positions externally adjacent to it (their p. 195). Holes are not mentioned; read literally, the + definition counts hole borders and nets the area. + +So, on Gaussian surfaces: + +- **Projected Kolmogorov density** (β = 11/3) gives 1.15, not 1.36. Single loops predict 13/12; the 7/6 sometimes + quoted is the formula for all contours together. +- **CO maps** have projected spectra near β = 2.8 (Stutzki et al. 1998: 2.8 for the Polaris Flare, 2.7 for another + field). Single loops there give 1.34 to 1.35 without holes and 1.52 ± 0.06 with them. Falgarone et al. fitted + perimeter against area. Without holes, their 1.36 ± 0.02 corresponds to β = 2.70 (2.57 to 2.86 with this + computation's own error added); the formula alone, without the loops' offset, puts it at 2.56. With holes, as + their definition reads, it corresponds to β = 3.18 (3.00 to 3.36). So if holes were left out, the agreement with + CO-map spectra leans on the offset; if they were counted, the edges imply a steeper spectrum than the maps'. +- **Measured values on real maps run from about 1.2 to 1.6.** Falgarone et al. list 1.26 ± 0.04 (Bazell and + Désert), 1.2 to 1.55 (Vogelaar et al.) and 1.18 to 1.28 (Dickman et al.). Sánchez et al. (2005) summarise + 1.2 ≲ D ≲ 1.5 and measure Orion A in ¹³CO at 1.32 ± 0.02; they count hole borders and report that including or + excluding holes did not change their result. On these surfaces their procedure still gives 1.31 without holes + and 1.40 with them at β = 2.8 (as #34's re-review computed), so the difference probably lies in the fields: + projected particle clouds against Gaussian surfaces. Stutzki et al. (1998) find d ≈ 1.6 for the Polaris Flare, + which they take as d = 2 − H after Voss, and typically 1.3 to 1.5 for other clouds. Counting holes, or whole + contour sets, raises the number, so it can account for the high end only. Values below the single loops' 1.34 + at β = 2.8 need steeper spectra (single loops give 1.29 at β = 3) or coarser maps (Sánchez et al. find D falls + as the number of pixels falls). +- **The calibration behind 2 − H.** Vogelaar and Wakker (1994) generated fractional-Brownian maps with contour + dimensions of 1.2 to 1.6 and recovered them to within 0.05 with a per-object perimeter-area code: 4-connected + pixel objects, a perimeter summed from pixel sides, with the diagonal instead of two sides at corner pixels + (p. 559), which, read literally, counts hole borders, and all levels pooled. Their text gives no mapping from D to + H and does not discuss holes; 2 − H is Stutzki et al.'s reading of that calibration. On these surfaces at 2048, + where 2 − H = 1.6 and 1.5 (β = 2.8 and 3), single loops give 1.35 and 1.29 without holes and 1.52 and 1.41 with + them, while Vogelaar and Wakker's Fig. 5 reads about 1.62 for a generated 1.58 and 1.48 for 1.50. A pixel + estimator in their style, but counting plain pixel sides without the corner rule, gives 1.42 and 1.37 with holes + (L = 256, β = 2.84 and 3; #34's re-review); a diagonal at convex corners adds about 0.03 (about 1.45 and 1.39), + and also cutting concave corners takes off about 0.015. No reading reproduces their calibration at D = 1.4 to + 1.6, and the cause (the generator, the pooling of levels, the map size) is left open. + +The laboratory values are sections of two kinds of surface: the interface between turbulent and quiet fluid, +explained as 7/3 (Sreenivasan, Ramshankar and Meneveau 1989), and iso-scalar surfaces inside the turbulence. An +interface bounds one region, so a single loop on a slice through a Kolmogorov field (β = 8/3, H = 1/3) is a natural +Gaussian analogue. That gives 1.36 here, and (3 − 1/3)/2 = 4/3 equals the 7/3 surface's section. But the section +dimension 1 + H and (3 − H)/2 coincide only at H = 1/3, and the interface is not a level set of a Gaussian field. +Iso-scalar surfaces have a dimension that falls as the level moves away from the mean: about 8/3 at the mean +(Constantin, Procaccia and Sreenivasan 1991), whose section 5/3 = 2 − H is the whole-level-set value. Iyer et al. +(2020) find 8/3 an upper bound reached there, a plateau near 2.35 at 1.5 standard deviations from the mean (which +they say matches Sreenivasan and Meneveau's experiments), and 2 by about three standard deviations. Cloud maps are +projections of supersonic gas whose spectrum is shallower than Kolmogorov's; the laboratory values come from +sections through incompressible turbulence. Both sit near 4/3 for different reasons, so their agreement does not +show shared physics. + +Real column densities are not Gaussian: they are close to log-normal and intermittent. Sánchez et al. (2005) also +find that the projected perimeter dimension depends on the resolution and on the cloud's 3D dimension. So this +settles which formula a perimeter-area fit over single loops measures, not what the clouds are. + +## Run it + +``` +node turbulence/table.mjs # everything above but the seed spreads, with the long run and the scan (about a minute) +node turbulence/seeds.mjs # the seed spreads: 32 seeds of the test's run and 8 of the scan (about 4 minutes) +node turbulence/table.mjs quick # the same with the test's short run, no scan (a few seconds) +node turbulence/edges.mjs # the cloud-edge table, 2048 points per side (about 45 s on 4 cores) +node turbulence/edges.mjs 4096 # any power of 2; 4096 takes about 4 minutes on 4 cores +node turbulence/edges.mjs quick # a smoke test of the same code at 256 points per side (about a second) +npm test # includes tests/turbulence.mjs and tests/turbulence-edges.mjs (about 5 s) +``` + +`model.js` is the cascade and the scales; `contours.js` is the random surfaces and their contours. Both load in +Node (`require`) and in a browser (`window.TURB`, `window.CONTOURS`). + +## Checks + +`tests/turbulence.mjs`: + +| Check | Expected | Model | +|---|---|---| +| Inviscid, unforced: energy and helicity with a + b + c = 0 | conserved | drift 6×10⁻¹¹ and 5×10⁻¹⁰ over t = 200 | +| The other sign convention, a + b − c = 0 | not conserved | E off by 50% at t = 0.9 | +| Dissipation over injection (model) | 1 ± 0.01 | 0.9999 | +| Flux through shells 0 to 6, against the injection and the plateau | 1 ± 0.012, flat to ± 0.006 | 1.0005, flat to 0.003 | +| Shell-model 4/5 law, shells 1 to 5 (L'vov et al.'s eq. 16; derived in `model.js`) | 1 − (−1/2)^n ± 0.01 | within 0.0004 (typically 0.0014 over seeds) | +| Anomalous scaling: ζ₂ − 2/3 and 2 − ζ₆ | over 5 standard errors | over 80 and 50 | +| ζ₂, ζ₄, ζ₆ in the test's run against the long run | 0.719, 1.255, 1.716 | 0.720, 1.254, 1.717 | +| ζ₂, ζ₄, ζ₆ against L'vov et al.'s Table I (moments of u) | within twice their errors | 0.000, 0.002, 0.011 apart | +| Model's ζ₂ against the measured ones (Küchler et al.) | over 1.5 errors above 0.693 ± 0.011 and 0.698 ± 0.011 | 0.720: 2.4 and 2.0 errors | +| Cascade end against k_η (flux and spectrum half-points) | within 1.5 shells | −1.03 and −0.07 | +| Upper-ocean η for ε = 10⁻⁸ to 10⁻⁷ W/kg (arithmetic) | 1.8 to 3.3 mm | 1.84, 3.28 mm | +| Air η for ε = 0.4 to 10⁻⁷ (arithmetic) | 0.3 mm to 1.4 cm | 0.30, 13.6 mm | +| Proton inertial length at 5 cm⁻³ (NRL formulary) | 228/√5 km | 101.8 km | +| Solar-wind proton mean free path at 10⁵ K and 4×10⁴ K, in inertial lengths (arithmetic) | over 10⁵ | 1.7×10⁶ and 2.7×10⁵ | +| V/(2π d_i), with and without sin θ_BV at the Parker angle, against observed ion breaks (Chen et al. 2014) | 0.1 to 1 Hz | 0.457 and 0.625 Hz | +| Sonic sizes from Larson and Solomon et al. against filament widths | within a factor of 2 of 0.10 pc | 0.099 and 0.147 pc | +| Larson's relation at 0.1 pc: non-thermal Mach number (this study's tolerance) | 0.8 to 1.25 | 1.01 | +| Larson's Table 1: TMC1 and TMC2 subsonic, L1407/L1257 and L1544/L63 supersonic, in his convention | as he says | 0.74, 0.44; 2.45, 1.74 | + +`tests/turbulence-edges.mjs`: + +| Check | Expected | Model | +|---|---|---| +| FFT against the direct DFT; 2D Parseval and round trip | equal | within 4×10⁻¹³ | +| Random streams: consecutive blocks, overlap refused, two uniforms per two normals | | yes | +| Surfaces against a naive synthesis; fitted spectral slopes | equal; β ± 0.06 | 2×10⁻¹⁴; within 0.012 | +| Circle: P, A, R_g; the same circle across the torus seams | 2πr₀, πr₀², r₀; identical | within 0.05% | +| Contours around the torus found and dropped | 2 and 2 | 2 and 2 | +| Saddle cells joined or split by the cell average | 1, 2, 1, 2 loops | 1, 2, 1, 2 | +| Green's theorem: loop areas add up to the area above the level (10 of 12 surfaces) | exact | 2×10⁻¹⁵ | +| Loops under a shift, a quarter turn, upside down | unchanged | within 7×10⁻¹³ | +| β = 2.4 and 3: loops nearer (3 − H)/2, level set nearer 2 − H | yes | yes | +| β = 2.4 and 3: the loops' offset above (3 − H_eff)/2 against the 0.04 found on larger grids | within 3 jackknife errors | 0.004 ± 0.018 and 0.030 ± 0.018 (this seed) | +| Fall of the loop dimension from β = 2.4 to 3 (paired jackknife) | over 2 errors; within 2 errors of half the rise in H_eff (0.118), beyond 3 from all of it (0.236) | 0.092 ± 0.024 (this seed; about 80% of seeds clear the last at this size) | + +## What the model leaves out + +- The shell model is a caricature of the Navier–Stokes cascade: one mode per octave, local interactions only, no + geometry. It has an exact flux and an anomalous scaling of its own, but its exponents are not those of real + flows, and nothing here carries its numbers over to them. +- Kolmogorov scales use published dissipation rates; the η values are arithmetic, not measurements. +- The solar-wind numbers are for one typical wind; the Parker angle stands in for the field-flow angle of each + interval. +- Larson's fit is extrapolated a little below its data to reach the sonic scale, Solomon et al.'s by more than a + decade. Larson's smallest objects were observed in HC₃N, HC₅N, NH₃ and H₂CO, not CO; his σ already holds the + gas's thermal motions. +- The multifractal cut-off is checked second-hand (Boffetta, Mazzino and Vulpiani 2008, citing Paladin and Vulpiani + 1987); Frisch's book was not reached. Kraichnan's ordering is from his abstract, where it is a conjecture. +- The cloud-edge computation uses Gaussian surfaces. Real column densities are close to log-normal and + intermittent, and observed maps add noise, beams and optical depth. The loops' offset above (3 − H_eff)/2 + depends on the set-up and is not explained. +- Sources were first read through search-result text. #34's review and re-review then read many of them; the + Sources list marks what was read where, and values taken from those reviews are marked so in the tests. + +## Sources + +DOIs and titles were confirmed from search results (publisher, ADS or arXiv pages, or citing papers' reference +lists), and checked against Crossref in #34's re-review; entries without one have none that was found. Each entry +says how it was read: in the paper during #34's review or re-review, in its abstract only, or through search-result +text only. + +- A. N. Kolmogorov, *The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers*, Dokl. Akad. Nauk SSSR 30 (1941); translation Proc. R. Soc. Lond. A 434, 9 (1991). Search text only. https://doi.org/10.1098/rspa.1991.0075 +- A. N. Kolmogorov, *Dissipation of energy in the locally isotropic turbulence*, Dokl. Akad. Nauk SSSR 32, 16 (1941); translation Proc. R. Soc. Lond. A 434, 15 (1991). Search text only. https://doi.org/10.1098/rspa.1991.0076 +- V. S. L'vov, E. Podivilov, A. Pomyalov, I. Procaccia and D. Vandembroucq, *Improved shell model of turbulence*, Phys. Rev. E 58, 1811 (1998). Read in #34's review. https://doi.org/10.1103/PhysRevE.58.1811 (https://arxiv.org/abs/chao-dyn/9803025) +- L. Biferale, *Shell models of energy cascade in turbulence*, Annu. Rev. Fluid Mech. 35, 441 (2003). Search text only. https://doi.org/10.1146/annurev.fluid.35.101101.161122 +- Z.-S. She and E. Leveque, *Universal scaling laws in fully developed turbulence*, Phys. Rev. Lett. 72, 336 (1994). Search text only. https://doi.org/10.1103/PhysRevLett.72.336 +- C. Küchler, E. Bodenschatz and G. P. Bewley, *Scaling in decaying turbulence at high Reynolds numbers*, preprint (2020; v2 2022). Read in #34's review. https://arxiv.org/abs/2006.10993 +- C. Küchler, G. P. Bewley and E. Bodenschatz, *Universal velocity statistics in decaying turbulence*, Phys. Rev. Lett. 131, 024001 (2023). Read in #34's review. https://doi.org/10.1103/PhysRevLett.131.024001 (https://arxiv.org/abs/2106.01795) +- G. Paladin and A. Vulpiani, *Degrees of freedom of turbulence*, Phys. Rev. A 35, 1971 (1987). Abstract only; the full text was not reached. https://doi.org/10.1103/PhysRevA.35.1971 +- G. Boffetta, A. Mazzino and A. Vulpiani, *Twenty-five years of multifractals in fully developed turbulence: a tribute to Giovanni Paladin*, J. Phys. A 41, 363001 (2008). Read in #34's review. https://doi.org/10.1088/1751-8113/41/36/363001 (https://arxiv.org/abs/0809.0196) +- U. Frisch, *Turbulence: The Legacy of A. N. Kolmogorov* (Cambridge University Press, 1995). Not reached. https://doi.org/10.1017/CBO9781139170666 +- R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Massaioli and S. Succi, *Extended self-similarity in turbulent flows*, Phys. Rev. E 48, R29 (1993). Abstract only (#34's re-review). https://doi.org/10.1103/PhysRevE.48.R29 +- K. R. Sreenivasan, *On the universality of the Kolmogorov constant*, Phys. Fluids 7, 2778 (1995). Search text only. https://doi.org/10.1063/1.868656 +- S. B. Pope, *Turbulent Flows* (Cambridge University Press, 2000). Search text only. https://doi.org/10.1017/CBO9780511840531 +- G. Gulitski et al., *Velocity and temperature derivatives in high-Reynolds-number turbulent flows in the atmospheric surface layer. Part 1. Facilities, methods and some general results*, J. Fluid Mech. 589, 57 (2007). Abstract only (#34's re-review); the source of R_λ 1,600 to 6,600. https://doi.org/10.1017/S0022112007007495 +- G. Gulitski et al., *Velocity and temperature derivatives in high-Reynolds-number turbulent flows in the atmospheric surface layer. Part 2. Accelerations and related matters*, J. Fluid Mech. 589, 83 (2007). Abstract only (#34's re-review). https://doi.org/10.1017/S0022112007007501 +- B. B. Balsley, M. Tjernström and G. Svensson, *Turbulence in the nocturnal boundary layer: highly-structured, strongly-variable, but ubiquitous*, AMS 18th Symposium on Boundary Layers and Turbulence, 18B.1 (2008). Search text only. https://ams.confex.com/ams/18BLT/techprogram/session_22073.htm +- W. Yang et al., *Direct measurements of turbulence in the upper western Pacific North Equatorial Current over a 25-h period*, Sensors 22, 1167 (2022). Read in #34's review. https://doi.org/10.3390/s22031167 +- G. D. Nastrom and K. S. Gage, *A climatology of atmospheric wavenumber spectra of wind and temperature observed by commercial aircraft*, J. Atmos. Sci. 42, 950 (1985). Search text only. https://doi.org/10.1175/1520-0469(1985)042<0950:ACOAWS>2.0.CO;2 +- R. H. Kraichnan, *Inertial ranges in two-dimensional turbulence*, Phys. Fluids 10, 1417 (1967). Abstract only. https://doi.org/10.1063/1.1762301 +- T. DeWitt and T. Garrett, *Global sonde datasets do not support a mesoscale transition in the turbulent energy cascade*, preprint (2025). Search text only. https://arxiv.org/abs/2510.23625 +- J. W. Armstrong, B. J. Rickett and S. R. Spangler, *Electron density power spectrum in the local interstellar medium*, ApJ 443, 209 (1995). Read in #34's review. https://doi.org/10.1086/175515 +- S. R. Spangler and C. R. Gwinn, *Evidence for an inner scale to the density turbulence in the interstellar medium*, ApJ 353, L29 (1990). Abstract only (#34's re-review, from the ADS scan). https://doi.org/10.1086/185700 +- A. Chepurnov and A. Lazarian, *Extending the big power law in the sky with turbulence spectra from Wisconsin Hα Mapper data*, ApJ 710, 853 (2010). Read in #34's review. https://doi.org/10.1088/0004-637X/710/1/853 (https://arxiv.org/abs/0905.4413) +- S. I. Braginskii, *Transport processes in a plasma*, in Reviews of Plasma Physics 1, 205 (Consultants Bureau, 1965). Search text only. +- NRL Plasma Formulary (Naval Research Laboratory). Search text only. https://www.nrl.navy.mil/ppd/content/nrl-plasma-formulary/ +- R. J. Leamon, C. W. Smith, N. F. Ness, W. H. Matthaeus and H. K. Wong, *Observational constraints on the dynamics of the interplanetary magnetic field dissipation range*, J. Geophys. Res. 103, 4775 (1998). Abstract only. https://doi.org/10.1029/97JA03394 +- R. J. Leamon et al., *MHD-driven kinetic dissipation in the solar wind and corona*, ApJ 537, 1054 (2000). Read in #34's review. https://doi.org/10.1086/309059 +- F. Sahraoui et al., *Evidence of a cascade and dissipation of solar-wind turbulence at the electron gyroscale*, Phys. Rev. Lett. 102, 231102 (2009). Read in #34's review. https://doi.org/10.1103/PhysRevLett.102.231102 +- O. Alexandrova et al., *Universality of solar-wind turbulent spectrum from MHD to electron scales*, Phys. Rev. Lett. 103, 165003 (2009). Read in #34's review. https://doi.org/10.1103/PhysRevLett.103.165003 (https://arxiv.org/abs/0906.3236) +- C. H. K. Chen, L. Leung, S. Boldyrev, B. A. Maruca and S. D. Bale, *Ion-scale spectral break of solar wind turbulence at high and low beta*, Geophys. Res. Lett. 41, 8081 (2014). Read in #34's review. https://doi.org/10.1002/2014GL062009 +- C. W. Smith, K. Hamilton, B. J. Vasquez and R. J. Leamon, *Dependence of the dissipation range spectrum of interplanetary magnetic fluctuations on the rate of energy cascade*, ApJ 645, L85 (2006). Read in #34's re-review. https://doi.org/10.1086/506151 +- R. Bruno and V. Carbone, *The solar wind as a turbulence laboratory*, Living Rev. Solar Phys. 10, 2 (2013). Read in #34's re-review. https://doi.org/10.12942/lrsp-2013-2 +- J. J. Podesta, D. A. Roberts and M. L. Goldstein, *Spectral exponents of kinetic and magnetic energy spectra in solar wind turbulence*, ApJ 664, 543 (2007). Search text only. https://doi.org/10.1086/519211 +- R. B. Larson, *Turbulence and star formation in molecular clouds*, MNRAS 194, 809 (1981). Read in #34's review and re-review. https://doi.org/10.1093/mnras/194.4.809 +- P. M. Solomon, A. R. Rivolo, J. Barrett and A. Yahil, *Mass, luminosity, and line width relations of Galactic molecular clouds*, ApJ 319, 730 (1987). Read in #34's review. https://doi.org/10.1086/165493 +- M. H. Heyer and C. M. Brunt, *The universality of turbulence in Galactic molecular clouds*, ApJ 615, L45 (2004). Read in #34's review. https://doi.org/10.1086/425978 (https://arxiv.org/abs/astro-ph/0409420) +- C. Federrath, R. S. Klessen, L. Iapichino and J. R. Beattie, *The sonic scale of interstellar turbulence*, Nature Astronomy 5, 365 (2021). Read in #34's review. https://doi.org/10.1038/s41550-020-01282-z (https://arxiv.org/abs/2011.06238) +- A. A. Goodman, J. A. Barranco, D. J. Wilner and M. H. Heyer, *Coherence in dense cores. II. The transition to coherence*, ApJ 504, 223 (1998). Search text only. https://doi.org/10.1086/306045 +- J. E. Pineda et al., *Direct observation of a sharp transition to coherence in dense cores*, ApJ 712, L116 (2010). Search text only. https://doi.org/10.1088/2041-8205/712/1/L116 (https://arxiv.org/abs/1002.2946) +- D. Arzoumanian et al., *Characterizing interstellar filaments with Herschel in IC 5146*, A&A 529, L6 (2011). Read in #34's review. https://doi.org/10.1051/0004-6361/201116596 (https://arxiv.org/abs/1103.0201) +- D. Arzoumanian et al., *Characterizing the properties of nearby molecular filaments observed with Herschel*, A&A 621, A42 (2019). Read in #34's re-review. https://doi.org/10.1051/0004-6361/201832725 (https://arxiv.org/abs/1810.00721) +- J. Maurer and P. Tabeling, *Local investigation of superfluid turbulence*, Europhys. Lett. 43, 29 (1998). Search text only. https://doi.org/10.1209/epl/i1998-00314-9 +- C. F. Barenghi, L. Skrbek and K. R. Sreenivasan, *Introduction to quantum turbulence*, PNAS 111 (Suppl. 1), 4647 (2014). Read in #34's review. https://doi.org/10.1073/pnas.1400033111 (https://arxiv.org/abs/1404.1909) +- E. Falgarone, T. G. Phillips and C. K. Walker, *The edges of molecular clouds: fractal boundaries and density structure*, ApJ 378, 186 (1991). Read in #34's review and re-review. https://doi.org/10.1086/170419 +- J. Stutzki, *The small scale structure of molecular clouds*, Rev. Mod. Astron. 6, 209 (1993), Fig. 5. Read in #34's review. +- J. Stutzki, F. Bensch, A. Heithausen, V. Ossenkopf and M. Zielinsky, *On the fractal structure of molecular clouds*, A&A 336, 697 (1998). Read in #34's review. +- M. G. R. Vogelaar and B. P. Wakker, *Measuring the fractal structure of interstellar clouds*, A&A 291, 557 (1994). Read in #34's re-review. +- N. Sánchez, E. J. Alfaro and E. Pérez, *The fractal dimension of projected clouds*, ApJ 625, 849 (2005). Read in #34's review and re-review. https://doi.org/10.1086/429553 (https://arxiv.org/abs/astro-ph/0501573) +- S. Lovejoy, *Area-perimeter relation for rain and cloud areas*, Science 216, 185 (1982). Search text only. https://doi.org/10.1126/science.216.4542.185 +- K. R. Sreenivasan and C. Meneveau, *The fractal facets of turbulence*, J. Fluid Mech. 173, 357 (1986). Abstract only (#34's review), which gives no numbers; its 1.36 ± 0.05 for iso-concentration and iso-dissipation sections is as Falgarone et al. quote it, jointly with Meneveau 1989. https://doi.org/10.1017/S0022112086001209 +- C. Meneveau, Ph.D. thesis, Yale University (1989). Not reached; cited as Falgarone et al. 1991 quote it. +- K. R. Sreenivasan, R. Ramshankar and C. Meneveau, *Mixing, entrainment and fractal dimensions of surfaces in turbulent flows*, Proc. R. Soc. Lond. A 421, 79 (1989). Search text only. https://doi.org/10.1098/rspa.1989.0004 +- P. Constantin, I. Procaccia and K. R. Sreenivasan, *Fractal geometry of isoscalar surfaces in turbulence: theory and experiments*, Phys. Rev. Lett. 67, 1739 (1991). Abstract only; it gives no number. https://doi.org/10.1103/PhysRevLett.67.1739 +- K. P. Iyer, J. Schumacher, K. R. Sreenivasan and P. K. Yeung, *Fractal iso-level sets in high-Reynolds-number scalar turbulence*, Phys. Rev. Fluids 5, 044501 (2020). Read in #34's review and re-review. https://doi.org/10.1103/PhysRevFluids.5.044501 (https://arxiv.org/abs/1910.00607) +- B. B. Mandelbrot, *The Fractal Geometry of Nature* (W. H. Freeman, 1982). Search text only. +- J. Kondev and C. L. Henley, *Geometrical exponents of contour loops on random Gaussian surfaces*, Phys. Rev. Lett. 74, 4580 (1995). Search text only. https://doi.org/10.1103/PhysRevLett.74.4580 +- J. Kondev, C. L. Henley and D. G. Salinas, *Nonlinear measures for characterizing rough surface morphologies*, Phys. Rev. E 61, 104 (2000). Read in #34's review. https://doi.org/10.1103/PhysRevE.61.104 (https://arxiv.org/abs/cond-mat/9907229) diff --git a/turbulence/contours.js b/turbulence/contours.js new file mode 100644 index 0000000..315f688 --- /dev/null +++ b/turbulence/contours.js @@ -0,0 +1,448 @@ +/* Contour loops of Gaussian random surfaces: which dimension a perimeter-area relation measures. No dependencies. + * + * Cloud edges in projected maps follow a perimeter-area relation P ~ A^(D/2) with D ~ 1.36 (Falgarone, Phillips & + * Walker, ApJ 378, 186 (1991)). For a Gaussian random surface with 2D power spectrum P(k) ~ |k|^-beta, so Hurst + * exponent H = (beta - 2)/2, two dimensions are quoted: + * all the contours at one level together D = 2 - H (Mandelbrot) + * one contour loop D = (3 - H) / 2 (Kondev & Henley, Phys. Rev. Lett. 74, 4580 (1995); a + * conjecture, x_l = 1/2, exact only at H = 0 and 1) + * This module makes the surfaces and measures both, so a perimeter-area fit can be compared with each. + * + * Field: complex white noise drawn in k space on an L x L periodic grid, multiplied by |k|^(-beta/2) (k in units of + * 2 pi/L, k = 0 removed, every other lattice mode kept), inverse FFT (own radix-2). The filter is real and even in k, + * so the real and imaginary parts of the result are two independent real fields with the same spectrum. Each is set to + * zero mean and unit variance. Random numbers: mulberry32, whose state advances by a fixed odd constant per draw, so + * pair p of a run takes its own consecutive stretch of the generator's cycle and no two pairs share a draw. + * Contours: marching squares, crossing points by linear interpolation along cell edges; a vertex is "above" when + * f > level. A saddle cell (diagonal corners above, the other two below) joins the above corners when the average of + * its four corners is above the level. Loops are walked with unwrapped coordinates, above side on the left, so the + * shoelace area is positive around an above island. A loop whose walk returns displaced by a multiple of L winds + * around the torus; it is counted and dropped from the loop statistics. + * Per loop: length P, signed area A (positive around an above island), and radius of gyration R_g of the polyline + * as a uniform wire. + */ +(function (root) { + 'use strict'; + + // ---------- random numbers ---------- + const STEP = 0x6D2B79F5; + function mulberry32(state) { // the i-th draw uses state + i STEP (mod 2^32); returns uniform [0, 1) + let a = state >>> 0; + return function () { + a = (a + STEP) >>> 0; + let t = a; + t = Math.imul(t ^ (t >>> 15), t | 1); + t ^= t + Math.imul(t ^ (t >>> 7), t | 61); + return ((t ^ (t >>> 14)) >>> 0) / 4294967296; + }; + } + // Block `block` of `size` draws from the cycle that starts at `seed`: disjoint from every other block while all + // blocks used together stay under 2^32 draws. + function stream(seed, block, size) { + if ((block + 1) * size > 4294967296) throw new Error('random stream blocks would overlap'); + return mulberry32((seed + Math.imul((block * size) >>> 0, STEP)) >>> 0); + } + function gaussian(rand) { // Box-Muller, both values used: exactly two draws per pair of normals + let spare = 0, has = false; + return function () { + if (has) { has = false; return spare; } + const r = Math.sqrt(-2 * Math.log(1 - rand())), th = 2 * Math.PI * rand(); + spare = r * Math.sin(th); has = true; return r * Math.cos(th); + }; + } + + // ---------- FFT ---------- + const plans = new Map(); + function plan(n) { + let p = plans.get(n); if (p) return p; + if (!(n >= 1) || (n & (n - 1))) throw new Error('FFT length must be a power of 2, got ' + n); + const bits = Math.round(Math.log2(n)), rev = new Uint32Array(n), h = Math.max(1, n >> 1); + const c = new Float64Array(h), s = new Float64Array(h); + for (let i = 0; i < n; i++) { let r = 0; for (let b = 0; b < bits; b++) r |= ((i >> b) & 1) << (bits - 1 - b); rev[i] = r; } + for (let i = 0; i < h; i++) { c[i] = Math.cos(2 * Math.PI * i / n); s[i] = Math.sin(2 * Math.PI * i / n); } + p = { n, rev, c, s }; plans.set(n, p); return p; + } + // In place, iterative radix-2. Forward: X_k = sum_n x_n e^(-2 pi i k n / N). Inverse includes the 1/N. + function fft(re, im, inverse) { + const n = re.length, { rev, c, s } = plan(n), sg = inverse ? 1 : -1; + for (let i = 0; i < n; i++) { + const j = rev[i]; + if (j > i) { let t = re[i]; re[i] = re[j]; re[j] = t; t = im[i]; im[i] = im[j]; im[j] = t; } + } + for (let size = 2; size <= n; size <<= 1) { + const half = size >> 1, step = n / size; + for (let i0 = 0; i0 < n; i0 += size) { + for (let j = 0, k = 0; j < half; j++, k += step) { + const a = i0 + j, b = a + half, wr = c[k], wi = sg * s[k]; + const tr = re[b] * wr - im[b] * wi, ti = re[b] * wi + im[b] * wr; + re[b] = re[a] - tr; im[b] = im[a] - ti; re[a] += tr; im[a] += ti; + } + } + } + if (inverse) { const q = 1 / n; for (let i = 0; i < n; i++) { re[i] *= q; im[i] *= q; } } + } + // 2D, row-major (index j L + i, i along x). Rows by fft(); columns by the same butterflies applied to whole rows. + function fft2(re, im, L, inverse) { + for (let j = 0; j < L; j++) fft(re.subarray(j * L, j * L + L), im.subarray(j * L, j * L + L), inverse); + const { rev, c, s } = plan(L), sg = inverse ? 1 : -1, tmp = new Float64Array(L); + for (let j = 0; j < L; j++) { + const r = rev[j]; + if (r > j) for (const A of [re, im]) { tmp.set(A.subarray(j * L, j * L + L)); A.copyWithin(j * L, r * L, r * L + L); A.set(tmp, r * L); } + } + for (let size = 2; size <= L; size <<= 1) { + const half = size >> 1, step = L / size; + for (let j0 = 0; j0 < L; j0 += size) { + for (let q = 0, k = 0; q < half; q++, k += step) { + const a = (j0 + q) * L, b = a + half * L, wr = c[k], wi = sg * s[k]; + for (let x = 0; x < L; x++) { + const br = re[b + x], bi = im[b + x], tr = br * wr - bi * wi, ti = br * wi + bi * wr; + re[b + x] = re[a + x] - tr; im[b + x] = im[a + x] - ti; re[a + x] += tr; im[a + x] += ti; + } + } + } + } + if (inverse) { const q = 1 / L, N = L * L; for (let i = 0; i < N; i++) { re[i] *= q; im[i] *= q; } } + } + + // ---------- Gaussian random surfaces ---------- + let filterKey = '', filterArr = null; + function filter(L, beta) { // |k|^(-beta/2), zero at k = 0 + const key = L + ':' + beta; if (key === filterKey) return filterArr; + const N = L * L, F = new Float64Array(N), e = -beta / 4; + for (let j = 0; j < L; j++) { + const ky = j <= L / 2 ? j : j - L; + for (let i = 0; i < L; i++) { const kx = i <= L / 2 ? i : i - L, k2 = kx * kx + ky * ky; F[j * L + i] = k2 ? Math.pow(k2, e) : 0; } + } + filterKey = key; filterArr = F; return F; + } + function normalise(f) { + const N = f.length; let m = 0; for (let i = 0; i < N; i++) m += f[i]; m /= N; + let v = 0; for (let i = 0; i < N; i++) { f[i] -= m; v += f[i] * f[i]; } + const q = 1 / Math.sqrt(v / N); for (let i = 0; i < N; i++) f[i] *= q; + return f; + } + // Pair `pair` of independent fields from the cycle starting at `seed` (2 L^2 draws per pair). The DFT of complex + // white noise is again i.i.d. complex Gaussian, so the noise is drawn directly in k space. ws = { re, im } reuses + // buffers between calls. + function gaussianFields(L, beta, seed, pair, ws) { + const N = L * L, re = ws && ws.re && ws.re.length === N ? ws.re : new Float64Array(N), im = ws && ws.im && ws.im.length === N ? ws.im : new Float64Array(N); + const g = gaussian(stream(seed, pair, 2 * N)); + for (let i = 0; i < N; i++) { re[i] = g(); im[i] = g(); } + const F = filter(L, beta); + for (let i = 0; i < N; i++) { re[i] *= F[i]; im[i] *= F[i]; } + fft2(re, im, L, true); + if (ws) { ws.re = re; ws.im = im; } + return [normalise(re), normalise(im)]; + } + // Radially averaged power per mode in shells of width 1, 1 <= |k| <= L/2: mean |k| and mean power. + function radialSpectrum(f, L) { + const N = L * L, re = Float64Array.from(f), im = new Float64Array(N); + fft2(re, im, L, false); + const K = L / 2, sk = new Float64Array(K + 1), sp = new Float64Array(K + 1), n = new Float64Array(K + 1); + for (let j = 0; j < L; j++) { + const ky = j <= L / 2 ? j : j - L; + for (let i = 0; i < L; i++) { + const kx = i <= L / 2 ? i : i - L, k = Math.sqrt(kx * kx + ky * ky), b = Math.round(k); + if (b < 1 || b > K) continue; + const v = j * L + i; sk[b] += k; sp[b] += re[v] * re[v] + im[v] * im[v]; n[b]++; + } + } + const k = [], P = []; + for (let b = 1; b <= K; b++) if (n[b]) { k.push(sk[b] / n[b]); P.push(sp[b] / n[b] / (N * N)); } + return { k, P }; + } + + // ---------- small statistics ---------- + function lsq(x, y) { // least squares y = intercept + slope x + const n = x.length; let mx = 0, my = 0; for (let i = 0; i < n; i++) { mx += x[i]; my += y[i]; } mx /= n; my /= n; + let sxx = 0, sxy = 0; for (let i = 0; i < n; i++) { sxx += (x[i] - mx) ** 2; sxy += (x[i] - mx) * (y[i] - my); } + const b = sxy / sxx; return { slope: b, intercept: my - b * mx, n }; + } + function wlsq(x, y, s) { // weighted least squares with errors s: slope, its error, intercept, chi^2 + let S = 0, Sx = 0, Sy = 0, Sxx = 0, Sxy = 0; + for (let i = 0; i < x.length; i++) { const w = 1 / (s[i] * s[i]); S += w; Sx += w * x[i]; Sy += w * y[i]; Sxx += w * x[i] * x[i]; Sxy += w * x[i] * y[i]; } + const d = S * Sxx - Sx * Sx, slope = (S * Sxy - Sx * Sy) / d, intercept = (Sxx * Sy - Sx * Sxy) / d; + let chi2 = 0; for (let i = 0; i < x.length; i++) chi2 += ((y[i] - intercept - slope * x[i]) / s[i]) ** 2; + return { slope, slopeErr: Math.sqrt(S / d), intercept, chi2 }; + } + function select(a, k) { // quickselect: a[k] becomes the k-th smallest + let lo = 0, hi = a.length - 1; + while (hi > lo) { + const mid = (lo + hi) >> 1; + if (a[mid] < a[lo]) { const t = a[mid]; a[mid] = a[lo]; a[lo] = t; } + if (a[hi] < a[lo]) { const t = a[hi]; a[hi] = a[lo]; a[lo] = t; } + if (a[hi] < a[mid]) { const t = a[hi]; a[hi] = a[mid]; a[mid] = t; } + const p = a[mid]; let i = lo, j = hi; + while (i <= j) { + while (a[i] < p) i++; while (a[j] > p) j--; + if (i <= j) { const t = a[i]; a[i] = a[j]; a[j] = t; i++; j--; } + } + if (k <= j) hi = j; else if (k >= i) lo = i; else break; + } + return a[k]; + } + function median(f) { const a = Float64Array.from(f), n = a.length, hi = select(a, n >> 1); if (n & 1) return hi; let lo = -Infinity; for (let i = 0; i < n >> 1; i++) if (a[i] > lo) lo = a[i]; return (lo + hi) / 2; } + + // ---------- marching squares ---------- + // Cell (i, j) has corners a = (i, j), b = (i+1, j), c = (i+1, j+1), d = (i, j+1) (bits 0-3 of the case) and sides + // 0 bottom (a-b), 1 right (b-c), 2 top (d-c), 3 left (a-d). Edge ids: the horizontal edge from vertex v to its right + // neighbour is 2v, the vertical edge from v to the vertex above is 2v + 1. EXIT[case * 4 + entry side] = exit side. + const EXIT = new Int8Array(64).fill(-1); + for (let cs = 0; cs < 16; cs++) { + const a = cs & 1, b = (cs >> 1) & 1, c = (cs >> 2) & 1, d = (cs >> 3) & 1; + const sides = [a !== b, b !== c, d !== c, a !== d].map((x, s) => (x ? s : -1)).filter(s => s >= 0); + if (sides.length === 2) { EXIT[cs * 4 + sides[0]] = sides[1]; EXIT[cs * 4 + sides[1]] = sides[0]; } + } + const CUT_BD = [1, 0, 3, 2]; // saddle: segments bottom-right and top-left, cutting off corners b and d (a, c joined) + const CUT_AC = [3, 2, 1, 0]; // saddle: segments left-bottom and right-top, cutting off corners a and c (b, d joined) + const joinedAt = level => (a, b, c, d) => (a + b + c + d) * 0.25 > level; + + // contours(f, L, level, opts) -> loops. opts.clouds: also group loops into connected above-level regions (clouds with + // their holes); opts.lengthOnly: totals only, no per-loop arrays. + function contours(f, L, level, opts) { + opts = opts || {}; + const N = L * L, M = L - 1; + if (L & M) throw new Error('L must be a power of 2'); + const above = new Uint8Array(N); + for (let v = 0; v < N; v++) above[v] = f[v] > level ? 1 : 0; + const seen = new Uint8Array(2 * N), joined = joinedAt(level), lengthOnly = !!opts.lengthOnly; + const out = { L, level, P: [], A: [], Rg: [], lab: [], nWrap: 0, wrapLength: 0, wrapLab: [], totalLength: 0, crossedEdges: 0, segments: 0, maxClosure: 0 }; + + function trace(i0, j0, kind, av, v0) { + const e0 = 2 * v0 + kind; + let I, J, sin, x0, y0, lab; + if (kind === 0) { + const vr = j0 * L + ((i0 + 1) & M); x0 = i0 + (level - f[v0]) / (f[vr] - f[v0]); y0 = j0; lab = av ? v0 : vr; + if (av) { I = i0; J = j0; sin = 0; } else { I = i0; J = j0 - 1; sin = 2; } + } else { + const vu = ((j0 + 1) & M) * L + i0; x0 = i0; y0 = j0 + (level - f[v0]) / (f[vu] - f[v0]); lab = av ? v0 : vu; + if (!av) { I = i0; J = j0; sin = 3; } else { I = i0 - 1; J = j0; sin = 1; } + } + seen[e0] = 1; + let qx = 0, qy = 0, P = 0, S = 0, Sx = 0, Sy = 0, S2 = 0, n = 0; + for (;;) { + const ci = I & M, cj = J & M, ci1 = (ci + 1) & M, cj1 = (cj + 1) & M; + const va = cj * L + ci, vb = cj * L + ci1, vc = cj1 * L + ci1, vd = cj1 * L + ci; + const cs = above[va] | (above[vb] << 1) | (above[vc] << 2) | (above[vd] << 3); + const ex = cs === 5 || cs === 10 ? ((cs === 5) === joined(f[va], f[vb], f[vc], f[vd]) ? CUT_BD : CUT_AC)[sin] : EXIT[cs * 4 + sin]; + let e, x, y; + if (ex === 0) { e = 2 * va; x = I + (level - f[va]) / (f[vb] - f[va]); y = J; } + else if (ex === 1) { e = 2 * vb + 1; x = I + 1; y = J + (level - f[vb]) / (f[vc] - f[vb]); } + else if (ex === 2) { e = 2 * vd; x = I + (level - f[vd]) / (f[vc] - f[vd]); y = J + 1; } + else if (ex === 3) { e = 2 * va + 1; x = I; y = J + (level - f[va]) / (f[vd] - f[va]); } + else throw new Error('marching squares: no exit from cell case ' + cs); + x -= x0; y -= y0; + const dx = x - qx, dy = y - qy, l = Math.sqrt(dx * dx + dy * dy); + if (!lengthOnly) { + const mx = 0.5 * (x + qx), my = 0.5 * (y + qy); + S += qx * y - x * qy; Sx += l * mx; Sy += l * my; S2 += l * (mx * mx + my * my + l * l / 12); + } + P += l; qx = x; qy = y; n++; + if (e === e0) break; + if (seen[e]) throw new Error('contour walk reached an edge already used by another loop'); + seen[e] = 1; + if (ex === 0) { J--; sin = 2; } else if (ex === 1) { I++; sin = 3; } else if (ex === 2) { J++; sin = 0; } else { I--; sin = 1; } + } + out.segments += n; out.totalLength += P; + const wx = Math.round(qx / L), wy = Math.round(qy / L); + out.maxClosure = Math.max(out.maxClosure, Math.abs(qx - wx * L), Math.abs(qy - wy * L)); + if (wx || wy) { out.nWrap++; out.wrapLength += P; out.wrapLab.push(lab); return; } + if (lengthOnly) return; + const cx = Sx / P, cy = Sy / P; + out.P.push(P); out.A.push(S / 2); out.Rg.push(Math.sqrt(Math.max(0, S2 / P - cx * cx - cy * cy))); out.lab.push(lab); + } + + for (let j = 0; j < L; j++) { + for (let i = 0; i < L; i++) { + const v = j * L + i, av = above[v]; + if (av !== above[j * L + ((i + 1) & M)]) { out.crossedEdges++; if (!seen[2 * v]) trace(i, j, 0, av, v); } + if (av !== above[((j + 1) & M) * L + i]) { out.crossedEdges++; if (!seen[2 * v + 1]) trace(i, j, 1, av, v); } + } + } + if (opts.clouds) out.clouds = clouds(f, above, L, level, out, joined); + return out; + } + + // Connected above-level regions under the same saddle rule. Each loop belongs to the region on its left. A bounded + // region ("cloud") has exactly one positive (outer) loop and no loop around the torus; its perimeter is the sum over + // its loops (outer plus holes) and its area the sum of their signed areas. + function clouds(f, above, L, level, loops, joined) { + const N = L * L, M = L - 1, parent = new Int32Array(N); + for (let v = 0; v < N; v++) parent[v] = v; + const find = x => { while (parent[x] !== x) { parent[x] = parent[parent[x]]; x = parent[x]; } return x; }; + const union = (a, b) => { a = find(a); b = find(b); if (a !== b) { if (a < b) parent[b] = a; else parent[a] = b; } }; + for (let j = 0; j < L; j++) { + for (let i = 0; i < L; i++) { + const va = j * L + i, vb = j * L + ((i + 1) & M), vd = ((j + 1) & M) * L + i, vc = ((j + 1) & M) * L + ((i + 1) & M); + if (above[va] && above[vb]) union(va, vb); + if (above[va] && above[vd]) union(va, vd); + const cs = above[va] | (above[vb] << 1) | (above[vc] << 2) | (above[vd] << 3); + if ((cs === 5 || cs === 10) && joined(f[va], f[vb], f[vc], f[vd])) { if (cs === 5) union(va, vc); else union(vb, vd); } + } + } + const reg = new Map(); + const get = r => { let g = reg.get(r); if (!g) { g = { P: 0, A: 0, nPos: 0, wrap: false }; reg.set(r, g); } return g; }; + for (let k = 0; k < loops.P.length; k++) { + const g = get(find(loops.lab[k])); g.P += loops.P[k]; g.A += loops.A[k]; + if (loops.A[k] > 0) g.nPos++; + } + for (const lab of loops.wrapLab) get(find(lab)).wrap = true; + const res = { P: [], A: [], dropped: 0 }; + for (const g of reg.values()) { + if (g.wrap || g.nPos !== 1) { res.dropped++; continue; } + res.P.push(g.P); res.A.push(g.A); + } + return res; + } + + // Area above the level, cell by cell, with the same interpolation and saddle rule; independent of the walk. + function aboveArea(f, L, level) { + const M = L - 1, join = joinedAt(level); let area = 0; + for (let j = 0; j < L; j++) { + for (let i = 0; i < L; i++) { + const va = j * L + i, vb = j * L + ((i + 1) & M), vd = ((j + 1) & M) * L + i, vc = ((j + 1) & M) * L + ((i + 1) & M); + const fa = f[va], fb = f[vb], fc = f[vc], fd = f[vd]; + const cs = (fa > level) | ((fb > level) << 1) | ((fc > level) << 2) | ((fd > level) << 3); + if (cs === 0) continue; if (cs === 15) { area += 1; continue; } + const t = (p, q) => (level - p) / (q - p); + const tb = (cs & 1) !== ((cs >> 1) & 1) ? t(fa, fb) : 0, tr = ((cs >> 1) & 1) !== ((cs >> 2) & 1) ? t(fb, fc) : 0; + const tt = ((cs >> 3) & 1) !== ((cs >> 2) & 1) ? t(fd, fc) : 0, tl = (cs & 1) !== ((cs >> 3) & 1) ? t(fa, fd) : 0; + const triA = 0.5 * tb * tl, triB = 0.5 * (1 - tb) * tr, triC = 0.5 * (1 - tr) * (1 - tt), triD = 0.5 * tt * (1 - tl); + switch (cs) { + case 1: area += triA; break; case 2: area += triB; break; case 4: area += triC; break; case 8: area += triD; break; + case 14: area += 1 - triA; break; case 13: area += 1 - triB; break; case 11: area += 1 - triC; break; case 7: area += 1 - triD; break; + case 3: area += 0.5 * (tl + tr); break; case 12: area += 0.5 * ((1 - tl) + (1 - tr)); break; + case 9: area += 0.5 * (tb + tt); break; case 6: area += 0.5 * ((1 - tb) + (1 - tt)); break; + case 5: case 10: { + const jn = join(fa, fb, fc, fd); + if (cs === 5) area += jn ? 1 - triB - triD : triA + triC; else area += jn ? 1 - triA - triC : triB + triD; + break; + } + } + } + } + return area; + } + + // Structure function S(r) = <(f(x + r) - f(x))^2> over x and y displacements, r = 1, 2, 4, ..., L/4. + function structure(f, L) { + const M = L - 1, out = []; + for (let r = 1; r <= L / 4; r *= 2) { + let s = 0; + for (let j = 0; j < L; j++) { + const row = j * L, up = ((j + r) & M) * L; + for (let i = 0; i < L; i++) { const v = f[row + i], a = f[row + ((i + r) & M)] - v, b = f[up + i] - v; s += a * a + b * b; } + } + out.push(s / (2 * L * L)); + } + return out; + } + function subsample(f, L, s) { // every s-th vertex in each direction + const n = L / s, g = new Float64Array(n * n); + for (let j = 0; j < n; j++) for (let i = 0; i < n; i++) g[j * n + i] = f[j * s * L + i * s]; + return g; + } + + // ---------- binned power-law fits ---------- + // Bins of width 1/4 octave in log2 x; per bin, the count and the sums of ln x and ln y. + const PER_OCT = 4, OFF = 64, NB = 160; + function newBins() { return { n: new Float64Array(NB), sx: new Float64Array(NB), sy: new Float64Array(NB) }; } + function binAdd(B, x, y) { + const b = Math.floor(PER_OCT * Math.log2(x)) + OFF; if (b < 0 || b >= NB) return; + B.n[b]++; B.sx[b] += Math.log(x); B.sy[b] += Math.log(y); + } + function binsCombine(list, minus) { // sum of bins, optionally minus one of them (for the jackknife) + const B = newBins(); + for (const X of list) for (let b = 0; b < NB; b++) { B.n[b] += X.n[b]; B.sx[b] += X.sx[b]; B.sy[b] += X.sy[b]; } + if (minus) for (let b = 0; b < NB; b++) { B.n[b] -= minus.n[b]; B.sx[b] -= minus.sx[b]; B.sy[b] -= minus.sy[b]; } + return B; + } + // Least squares of bin-mean ln y on bin-mean ln x over bins centred in [lo, hi] with at least minN entries. + function binFit(B, lo, hi, minN) { + const xs = [], ys = []; + for (let b = 0; b < NB; b++) { + if (B.n[b] < (minN || 1)) continue; + const centre = Math.pow(2, (b - OFF + 0.5) / PER_OCT); + if (centre < lo || centre > hi) continue; + xs.push(B.sx[b] / B.n[b]); ys.push(B.sy[b] / B.n[b]); + } + if (xs.length < 3) return { slope: NaN, n: xs.length }; + return lsq(xs, ys); + } + + // ---------- one field, and the estimates over many ---------- + // What a perimeter-area study needs from one field: loops at the mean level (P against R_g and against sqrt A, and + // sqrt A against R_g), the outlines of above-level regions only (loops with the region inside, A > 0), loops at the + // median level, clouds with their holes, the total contour length on every s-th vertex (the level set's yardstick), + // and the structure function. + function fieldRecord(f, L, { yardMax = 32 } = {}) { + const r = contours(f, L, 0, { clouds: true }), R = newBins(), A = newBins(), RA = newBins(), Ro = newBins(); + for (let k = 0; k < r.P.length; k++) { + const P = r.P[k], a = Math.abs(r.A[k]), g = r.Rg[k]; + if (!(g > 0) || !(a > 0)) continue; + binAdd(R, g, P); binAdd(A, Math.sqrt(a), P); binAdd(RA, g, Math.sqrt(a)); + if (r.A[k] > 0) binAdd(Ro, g, P); + } + const CL = newBins(); for (let k = 0; k < r.clouds.P.length; k++) binAdd(CL, Math.sqrt(r.clouds.A[k]), r.clouds.P[k]); + const m = median(f), q = contours(f, L, m), Rm = newBins(); + for (let k = 0; k < q.P.length; k++) if (q.Rg[k] > 0) binAdd(Rm, q.Rg[k], q.P[k]); + const len = [r.totalLength]; + for (let s = 2; s <= yardMax; s *= 2) len.push(s * contours(subsample(f, L, s), L / s, 0, { lengthOnly: true }).totalLength); + return { R, A, RA, Ro, CL, Rm, len, sf: structure(f, L), nLoops: r.P.length, nWrap: r.nWrap, wrapShare: r.wrapLength / r.totalLength }; + } + // Fit windows. Loops: R_g in [r0, L/16] (primary; r0 = 8 from L = 1024, else 4), and Rmin in {r0/2, r0, 2 r0}, Rmax in + // {L/32, L/16, L/8}, at least two octaves wide, for the spread. Yardstick: s in [y0, ymax] and neighbours. + function loopWindows(L) { + const r0 = L >= 1024 ? 8 : 4, list = [[r0, L / 16]]; + for (const a of [r0 / 2, r0, 2 * r0]) for (const b of [L / 32, L / 16, L / 8]) if (b / a >= 4 && !(a === r0 && b === L / 16)) list.push([a, b]); + return list; + } + function yardWindows(L, ymax) { const y0 = L >= 1024 ? 4 : 2; return [[y0, ymax], [y0 / 2, ymax], [2 * y0, ymax], [y0 / 2, ymax / 2], [y0, ymax / 2]].filter(w => w[1] / w[0] >= 2); } + // Each estimate: value in the primary window, jackknife error over fields, and half the spread over windows. + // loops: D from P ~ R_g^D; area: D from P ~ A^(D/2) (sqrt A window = c x the R_g window, c = ); + // outer: outlines of above-level regions only; median: loops at the median level; + // clouds: perimeter including holes against net area; + // yard: level set, total length l(s) ~ s^(1 - D); Heff: half the log slope of the structure function, over the + // loop window (HeffLoop) and the yardstick window (HeffYard). + function estimate(recs, L, yardMax = 32) { + const tot = {}; + const pool = (name, ex) => { if (!tot[name]) tot[name] = binsCombine(recs.map(r => r[name])); return ex == null ? tot[name] : binsCombine([tot[name]], recs[ex][name]); }; + const ratio = (ex, lo, hi) => { + const B = pool('RA', ex); let n = 0, s = 0; + for (let b = 0; b < NB; b++) { const c = Math.pow(2, (b - OFF + 0.5) / PER_OCT); if (c >= lo && c <= hi) { n += B.n[b]; s += B.sy[b] - B.sx[b]; } } + return Math.exp(s / n); + }; + const meanLog = (key, ex) => { const v = recs.filter((_, i) => i !== ex).map(r => r[key]); return v[0].map((_, k) => v.reduce((a, x) => a + Math.log(x[k]), 0) / v.length); }; + const slopeIn = (m, lo, hi) => { const x = [], y = []; m.forEach((v, k) => { const s = Math.pow(2, k); if (s >= lo && s <= hi) { x.push(Math.log(s)); y.push(v); } }); return lsq(x, y).slope; }; + const est = { + loops: (w, ex) => binFit(pool('R', ex), w[0], w[1], 2).slope, + area: (w, ex) => { const c = ratio(ex, w[0], w[1]); return binFit(pool('A', ex), c * w[0], c * w[1], 2).slope; }, + outer: (w, ex) => binFit(pool('Ro', ex), w[0], w[1], 2).slope, + median: (w, ex) => binFit(pool('Rm', ex), w[0], w[1], 2).slope, + clouds: (w, ex) => { const c = ratio(ex, w[0], w[1]); return binFit(pool('CL', ex), c * w[0], c * w[1], 2).slope; }, + yard: (w, ex) => 1 - slopeIn(meanLog('len', ex), w[0], w[1]), + Heff: (w, ex) => 0.5 * slopeIn(meanLog('sf', ex), w[0], w[1]), + }; + const n = recs.length; + // Jackknife over fields; loo keeps the leave-one-out values, so estimates from fields that share random numbers + // (the same pair index at another beta) can be differenced or fitted jointly. byWindow: the estimate in each window. + const jack = (fn, w) => { const full = fn(w), v = recs.map((_, i) => fn(w, i)), mv = v.reduce((a, x) => a + x, 0) / n; return { D: full, se: Math.sqrt((n - 1) / n * v.reduce((a, x) => a + (x - mv) * (x - mv), 0)), loo: v }; }; + const withWindows = (fn, ws) => { const byWindow = ws.map(w => fn(w)), all = byWindow.filter(Number.isFinite); return { ...jack(fn, ws[0]), win: (Math.max(...all) - Math.min(...all)) / 2, byWindow }; }; + const W = loopWindows(L), YW = yardWindows(L, yardMax); + return { + loops: withWindows(est.loops, W), area: withWindows(est.area, W), outer: withWindows(est.outer, W), median: withWindows(est.median, W), clouds: withWindows(est.clouds, W), + yard: withWindows(est.yard, YW), HeffLoop: withWindows(est.Heff, W), HeffYard: jack(est.Heff, YW[0]), c: ratio(null, W[0][0], W[0][1]), W, YW, + heffAt: k => 0.5 * slopeIn(meanLog('sf'), k * W[0][0], k * W[0][1]), // H_eff over separations k times the loop window + nLoops: recs.reduce((a, r) => a + r.nLoops, 0) / n, nWrap: recs.reduce((a, r) => a + r.nWrap, 0) / n, wrapShare: recs.reduce((a, r) => a + r.wrapShare, 0) / n, + }; + } + + const theory = beta => { const H = (beta - 2) / 2; return { H, loop: (3 - H) / 2, set: 2 - H }; }; + + const api = { + mulberry32, stream, gaussian, fft, fft2, gaussianFields, normalise, radialSpectrum, lsq, wlsq, median, select, + contours, aboveArea, structure, subsample, newBins, binAdd, binsCombine, binFit, PER_OCT, OFF, + fieldRecord, loopWindows, yardWindows, estimate, theory, + }; + if (typeof module !== 'undefined' && module.exports) module.exports = api; else root.CONTOURS = api; +})(typeof window !== 'undefined' ? window : globalThis); diff --git a/turbulence/edges.mjs b/turbulence/edges.mjs new file mode 100644 index 0000000..ba24fd8 --- /dev/null +++ b/turbulence/edges.mjs @@ -0,0 +1,102 @@ +// Cloud edges: which dimension a perimeter-area relation measures, on Gaussian random surfaces. +// Usage: node turbulence/edges.mjs (L = 2048, 32 fields per spectrum, up to 4 worker threads; about 45 s on 4 cores) +// node turbulence/edges.mjs 1024 (any power of 2 from 256; 4096 takes about 4 minutes on 4 cores) +// node turbulence/edges.mjs quick (L = 256, 12 fields per spectrum, one thread: a smoke test, about a second) +import { createRequire } from 'node:module'; +import { Worker, isMainThread, parentPort } from 'node:worker_threads'; +import os from 'node:os'; +const C = createRequire(import.meta.url)('./contours.js'); + +const SEED = 2026, BETAS = [2.4, 8 / 3, 2.8, 3, 3.4, 11 / 3], FALGARONE = [1.36, 0.02]; +const job = ({ L, beta, pair, yardMax }, ws) => C.gaussianFields(L, beta, SEED, pair, ws).map(f => C.fieldRecord(f, L, { yardMax })); + +if (!isMainThread) { + const ws = {}; + parentPort.on('message', m => (m === 'done' ? process.exit(0) : parentPort.postMessage({ id: m.id, recs: job(m, ws) }))); +} else { + const arg = process.argv[2], quick = arg === 'quick', size = quick ? 256 : Number(arg || 2048); + if (!(size >= 256) || (size & (size - 1))) throw new Error('usage: node turbulence/edges.mjs [quick | L, a power of 2 from 256]'); + const cfg = quick ? { L: 256, pairs: 6, yardMax: 8, workers: 1 } : { L: size, pairs: 16, yardMax: Math.min(32, size / 32), workers: Math.min(4, os.availableParallelism ? os.availableParallelism() : os.cpus().length) }; + if (quick) console.log('Smoke test: at L = 256 the loops are fitted over 4 to 16 cells, too few to tell (3 - H)/2 from 2 - H. The numbers below\nare not the study\'s results; run without "quick" for those.\n'); + const jobs = []; + for (const beta of BETAS) for (let pair = 0; pair < cfg.pairs; pair++) jobs.push({ L: cfg.L, beta, pair, yardMax: cfg.yardMax }); + const t0 = Date.now(), results = await run(jobs, cfg.workers); + const f3 = x => x.toFixed(3), name = b => (Math.abs(b - 8 / 3) < 1e-9 ? '8/3' : Math.abs(b - 11 / 3) < 1e-9 ? '11/3' : String(b)); + const pm = r => `${f3(r.D)} +- ${f3(r.se)} +- ${f3(r.win)}`, sig = r => Math.hypot(r.se, r.win); + console.log(`Gaussian random surfaces, power spectrum |k|^-beta, H = (beta - 2)/2: L = ${cfg.L}, ${2 * cfg.pairs} fields per spectrum ` + + `(the same random numbers for each), contours at the mean level (${((Date.now() - t0) / 1e3).toFixed(1)} s, ${cfg.workers} thread${cfg.workers > 1 ? 's' : ''})`); + console.log(`Loops fitted over R_g in [${C.loopWindows(cfg.L)[0]}] cells; errors: jackknife over fields, then half the spread over ${C.loopWindows(cfg.L).length} fit windows.\n`); + const rows = BETAS.map(beta => { const recs = results.filter((_, i) => jobs[i].beta === beta).flat(); return { beta, th: C.theory(beta), e: C.estimate(recs, cfg.L, cfg.yardMax) }; }); + + console.log('One loop: perimeter against size, P ~ R_g^D, and against area, P ~ A^(D/2); outlines of above-level regions only;'); + console.log('the median level; clouds counted with their holes (perimeter of outline and holes against net area)'); + console.log(' beta H (3-H)/2 2-H P ~ R_g^D P ~ A^(D/2) outlines only median level, P~R_g clouds with holes'); + for (const { beta, th, e } of rows) { + console.log(` ${name(beta).padEnd(5)} ${f3(th.H)} ${f3(th.loop).padStart(8)} ${f3(th.set).padStart(6)} ${pm(e.loops).padEnd(24)} ${pm(e.area).padEnd(24)} ${pm(e.outer).padEnd(24)} ${pm(e.median).padEnd(24)} ${pm(e.clouds)}`); + } + console.log('\nAll the contours at the level together: total length on every s-th grid point, l(s) ~ s^(1 - D)'); + console.log(' beta yardstick D 2-H H_eff (yardstick range) 2-H_eff'); + for (const { beta, th, e } of rows) console.log(` ${name(beta).padEnd(5)} ${pm(e.yard).padEnd(24)} ${f3(th.set)} ${f3(e.HeffYard.D).padStart(10)} ${f3(2 - e.HeffYard.D)}`); + + console.log('\nOn a finite periodic grid the roughness the loops see differs from (beta - 2)/2: higher for flat spectra, lower for steep'); + console.log('ones. H_eff is half the log slope of the structure function <(f(x + r) - f(x))^2> over the loop range. Loops against each'); + console.log('formula, z = (D - formula) / sigma:'); + console.log(' beta H_eff (3-H_eff)/2 2-H_eff z vs (3-H)/2 z vs 2-H z vs (3-H_eff)/2 z vs 2-H_eff'); + for (const { beta, th, e } of rows) { + const h = e.HeffLoop.D, z = x => ((e.loops.D - x) / sig(e.loops)).toFixed(1); + console.log(` ${name(beta).padEnd(5)} ${f3(h)} ${f3((3 - h) / 2).padStart(11)} ${f3(2 - h).padStart(8)} ${z(th.loop).padStart(13)} ${z(th.set).padStart(9)} ${z((3 - h) / 2).padStart(17)} ${z(2 - h).padStart(13)}`); + } + const wts = rows.map(r => sig(r.e.loops)), slopeOf = (H, D) => C.wlsq(H, D, wts).slope; + const fit = C.wlsq(rows.map(r => r.e.HeffLoop.D), rows.map(r => r.e.loops.D), wts); + // The six spectra share their random numbers, so their errors are correlated: a joint jackknife refits the slope + // with the same field left out of every spectrum, and the window spread refits it with every spectrum in one window. + const nF = rows[0].e.loops.loo.length, jk = Array.from({ length: nF }, (_, i) => slopeOf(rows.map(r => r.e.HeffLoop.loo[i]), rows.map(r => r.e.loops.loo[i]))); + const jm = jk.reduce((a, x) => a + x, 0) / nF, jse = Math.sqrt((nF - 1) / nF * jk.reduce((a, x) => a + (x - jm) ** 2, 0)); + const byW = rows[0].e.loops.byWindow.map((_, k) => slopeOf(rows.map(r => r.e.HeffLoop.byWindow[k]), rows.map(r => r.e.loops.byWindow[k]))).filter(Number.isFinite); + const wse = (Math.max(...byW) - Math.min(...byW)) / 2, tot = Math.hypot(jse, wse); + const off = rows.map(r => r.e.loops.D - (3 - r.e.HeffLoop.D) / 2), w = rows.map(r => 1 / sig(r.e.loops) ** 2); + const wmean = v => v.reduce((a, x, i) => a + x * w[i], 0) / w.reduce((a, x) => a + x, 0), meanOff = wmean(off); + console.log(`\nAcross the six spectra, D against H_eff: slope ${f3(fit.slope)} +- ${f3(jse)} (joint jackknife) +- ${f3(wse)} (fit windows); ` + + `Kondev-Henley -1/2, Mandelbrot -1, ${((fit.slope + 1) / tot).toFixed(1)} sigma from -1. The weighted fit alone gives +- ${f3(fit.slopeErr)}, chi^2 ` + + `${fit.chi2.toFixed(1)} for 4 degrees of freedom, but treats the spectra as independent.`); + console.log(`Loops sit ${meanOff >= 0 ? '+' : ''}${f3(meanOff)} above (3-H_eff)/2 on average (weighted), with H_eff over the loops' own range of sizes;` + + ` measured over separations 2 and 4 times larger, the offset is ${[2, 4].map(k => f3(wmean(rows.map(r => r.e.loops.D - (3 - r.e.heffAt(k)) / 2)))).join(' and ')}.`); + + // Which spectrum gives Falgarone et al.'s 1.36 +- 0.02? They fitted perimeter against area, so use the P ~ A^(D/2) + // column; it falls with beta, so interpolate linearly, and widen the band by the column's own error in quadrature. + const col = rows.map(r => ({ beta: r.beta, D: r.e.area.D, s: sig(r.e.area) })); + const at = D => { for (let i = 0; i + 1 < col.length; i++) { const a = col[i], b = col[i + 1]; if ((a.D - D) * (b.D - D) <= 0) { const t = (D - a.D) / (b.D - a.D); return { beta: a.beta + t * (b.beta - a.beta), s: a.s + t * (b.s - a.s) }; } } return { beta: NaN, s: NaN }; }; + const [D0, dD] = FALGARONE, mid = at(D0), band = Math.hypot(dD, mid.s), fb = x => (Number.isFinite(x) ? f3(x) : 'outside the spectra run'); + console.log(`\nCloud edges, D = ${D0} +- ${dD} (Falgarone, Phillips & Walker 1991, perimeter against area; +- ${dD} is a combined-fit error): if holes are` + + ` not counted, the P ~ A^(D/2) column puts it at beta = ${fb(mid.beta)}, and ${fb(at(D0 + band).beta)} to ${fb(at(D0 - band).beta)} with this run's own error (+- ${f3(mid.s)}) added.`); + { // Their stated definition (perimeter = positions externally adjacent to the contour) counts hole borders, with net area. + const colH = rows.map(r => ({ beta: r.beta, D: r.e.clouds.D, s: sig(r.e.clouds) })); + const atH = D => { for (let i = 0; i + 1 < colH.length; i++) { const a = colH[i], b = colH[i + 1]; if ((a.D - D) * (b.D - D) <= 0) { const t = (D - a.D) / (b.D - a.D); return { beta: a.beta + t * (b.beta - a.beta), s: a.s + t * (b.s - a.s) }; } } return { beta: NaN, s: NaN }; }; + const mH = atH(D0), bH = Math.hypot(dD, mH.s), c28 = rows.find(r => Math.abs(r.beta - 2.8) < 1e-9).e.clouds; + console.log(`If holes are counted, as their definition reads, the clouds-with-holes column puts it at beta = ${fb(mH.beta)}, ${fb(atH(D0 + bH).beta)} to ` + + `${fb(atH(D0 - bH).beta)}; at beta = 2.8 that column gives ${f3(c28.D)} +- ${f3(sig(c28))}.`); + } + console.log(`Without the loops' offset, (3-H)/2 alone would put 1.36 at beta = ${f3(8 - 4 * D0)}, and 2-H at ${f3(6 - 2 * D0)}.`); + const row = b => { const e = rows.find(r => Math.abs(r.beta - b) < 1e-9).e.area; return `${f3(e.D)} +- ${f3(sig(e))}`; }; + console.log(`P ~ A^(D/2): projected Kolmogorov density (beta 11/3) ${row(11 / 3)}; a planar slice of it (beta 8/3) ${row(8 / 3)}; CO maps (beta ~ 2.8) ${row(2.8)}.`); + console.log(`Per field, loops of all sizes: ${rows.map(r => `${name(r.beta)}: ${r.e.nLoops.toFixed(0)}, ${r.e.nWrap.toFixed(2)} around the torus`).join('; ')}.`); + const vals = rows.flatMap(({ e }) => ['loops', 'area', 'outer', 'median', 'clouds', 'yard'].flatMap(k => [e[k].D, e[k].se, e[k].win, ...e[k].byWindow]) + .concat([e.HeffLoop.D, e.HeffYard.D, ...e.HeffLoop.byWindow, e.heffAt(2), e.heffAt(4)])).concat([fit.slope, jse, wse, meanOff]); + if (!vals.every(Number.isFinite)) { console.error('non-finite result'); process.exit(1); } +} + +async function run(jobs, nWorkers) { + const out = new Array(jobs.length); + if (nWorkers <= 1) { const ws = {}; jobs.forEach((j, i) => { out[i] = job(j, ws); }); return out; } + await new Promise((resolve, reject) => { + let next = 0, done = 0; + const feed = w => { if (next < jobs.length) { const i = next++; w.postMessage({ id: i, ...jobs[i] }); } else w.postMessage('done'); }; + for (let k = 0; k < Math.min(nWorkers, jobs.length); k++) { + const w = new Worker(new URL(import.meta.url)); + w.on('message', m => { out[m.id] = m.recs; if (++done === jobs.length) resolve(); feed(w); }); + w.on('error', reject); + feed(w); + } + }); + return out; +} diff --git a/turbulence/model.js b/turbulence/model.js new file mode 100644 index 0000000..3f5bd50 --- /dev/null +++ b/turbulence/model.js @@ -0,0 +1,336 @@ +/* Turbulence model: a pattern that repeats over many decades of scale, and the sizes where it stops. No dependencies. + * + * 1. Kolmogorov (1941): energy fed in at a large scale L cascades to smaller scales at a rate eps per unit mass, + * the same at every scale in between, until viscosity nu takes it at eta = (nu^3/eps)^(1/4). Between the two the + * velocity difference across a distance r depends only on eps and r, so the statistics repeat at every scale. + * 2. The cascade itself, in the Sabra shell model (L'vov, Podivilov, Pomyalov, Procaccia & Vandembroucq, "Improved shell + * model of turbulence", Phys. Rev. E 58, 1811 (1998)): one complex velocity u_n per octave of wavenumber, + * du_n/dt = i (a k_{n+1} u_{n+2} u*_{n+1} + b k_n u_{n+1} u*_{n-1} - c k_{n-1} u_{n-1} u_{n-2}) - nu k_n^2 u_n + f_n + * on shells k_n = k0 lambda^n, n = 0..N-1, with u_{-2} = u_{-1} = u_N = u_{N+1} = 0. + * Invariants. Write T_m = u*_m u*_{m+1} u_{m+2} (one triad per m). The a-term of shell m, the b-term of shell m+1 + * and the c-term of shell m+2 all carry T_m (the c-term as its conjugate) with the same k_{m+1}, so for any q + * d/dt sum q^n |u_n|^2 (nonlinear part) = -2 sum_m k_{m+1} q^m (a + b q + c q^2) Im T_m. + * So E = sum |u_n|^2 is conserved if a + b + c = 0, and then H = sum (a/c)^n |u_n|^2 is too. With a = 1, + * b = c = -1/2 and lambda = 2 (the 3D-like case), a/c = -2 and H = sum (-1)^n (k_n/k0) |u_n|^2, a helicity. + * Energy flux out of shells 0..n (energy normalised as E = sum |u_n|^2, not half of it): + * Pi_n = 2 [a k_{n+1} Im T_n - c k_n Im T_{n-1}] + * With the energy flux eps and the helicity flux delta both constant, Sigma_n = obeys exactly + * 2 (a - c) k_n Sigma_n = eps - delta (c/a)^n, i.e. 3 k_n Sigma_n = eps - delta (-1/2)^n, + * the shell model's 4/5 law: a third-order correlator that scales exactly as 1/k_n. Forcing on shell 0 alone puts + * helicity in at delta = eps. The relation is derived here from the equation above; it is L'vov et al.'s eq. (16), + * S_3 = [-eps + delta (c/a)^n] / (2 k_n (a - c)) with Sigma_n = -S_3 (as read from the paper in #34's review). + * Integrator: integrating-factor RK4 with the viscous factor exact. The exponentials are precomputed by an own + * + - * / routine, so the trajectory uses only IEEE + - * / and is the same bit for bit on any conforming engine. + * Model units: k0 = 1. + * 3. Where real cascades stop, in SI units: the ion inertial length and gyroradius of the solar wind (a plasma too + * dilute for viscosity), and the sonic scale of molecular clouds, where the turbulent velocity falls to the sound + * speed: the supersonic range ends there, and a shallower subsonic cascade continues below it. + */ +(function (root) { + 'use strict'; + const C = { + c: 2.99792458e8, kB: 1.380649e-23, e: 1.602176634e-19, eps0: 8.8541878128e-12, mu0: 1.25663706212e-6, + mp: 1.67262192369e-27, mH: 1.6735575e-27, pc: 3.0856775814913673e16, AU: 1.495978707e11, + }; + + // ---------- 1. Kolmogorov scales ---------- + // Kinematic viscosity, m^2/s: air at 20 C and 1 atm; seawater at 20 C and salinity 35 g/kg. + const NU = { air: 1.51e-5, seawater: 1.05e-6 }; + // Dissipation scale, time and velocity for kinematic viscosity nu (m^2/s) and dissipation rate eps (W/kg = m^2/s^3). + const kolmogorov = (nu, eps) => ({ eta: Math.pow(nu * nu * nu / eps, 0.25), tau: Math.sqrt(nu / eps), v: Math.pow(nu * eps, 0.25) }); + // Taylor-scale Reynolds number to the ratio of the large scale to eta, with Pope's conventions (Turbulent Flows, + // 2000): L = k^(3/2)/eps, Re_L = k^(1/2) L/nu = (3/20) Re_lambda^2, L/eta = Re_L^(3/4). + const outerOverEta = reLambda => Math.pow(0.15 * reLambda * reLambda, 0.75); + + // ---------- 2. Sabra shell model ---------- + function mulberry32(seed) { // small seeded generator so runs are repeatable + return function () { seed |= 0; seed = seed + 0x6D2B79F5 | 0; let t = Math.imul(seed ^ seed >>> 15, 1 | seed); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; }; + } + // A random point on the unit circle from + - * / only: (1 - t^2, 2t)/(1 + t^2), t in [-1, 1), then a random half-plane. + function unitPhase(rng) { + const t = 2 * rng() - 1, d = 1 + t * t, s = rng() < 0.5 ? -1 : 1; + return [s * (1 - t * t) / d, s * 2 * t / d]; + } + // exp(x) from + - * / only (scaling and squaring of a Taylor series): within 2e-13 relative for |x| < 60, which + // covers every factor used here (up to |x| = 55), and about 2e-12 near -700. + function expDet(x) { + if (x > 0) return 1 / expDet(-x); + if (x < -745) return 0; + let m = 0; while (x < -0.25) { x /= 2; m++; } + let s = 1, t = 1; for (let j = 1; j <= 18; j++) { t = t * x / j; s += t; } + for (let j = 0; j < m; j++) s *= s; + return s; + } + + const SABRA = { a: 1, b: -0.5, c: -0.5, lambda: 2, k0: 1 }; + const FORCE = [{ n: 0, re: 0.005, im: 0.005 }]; // constant force on shell 0 + // Runs. CI: about 2.5 s, about 860 large-eddy times sampled (L/U = 5.6). TABLE: about 45 s. SCAN: the viscosity scan. + const RUNS = { + CI: { N: 18, nu: 1e-6, dt: 2e-3, transient: 200, blockTime: 300, nBlocks: 16, every: 10, seed: 1, fit: [3, 8] }, + TABLE: { N: 21, nu: 1e-7, dt: 5e-4, transient: 500, blockTime: 1000, nBlocks: 20, every: 10, seed: 1, fit: [3, 11] }, + SCAN: [ + { N: 16, nu: 1e-5, dt: 4e-3, transient: 300, blockTime: 300, nBlocks: 8, every: 5, seed: 1, fit: [3, 6] }, + { N: 18, nu: 1e-6, dt: 2e-3, transient: 300, blockTime: 300, nBlocks: 8, every: 5, seed: 1, fit: [3, 8] }, + { N: 21, nu: 1e-7, dt: 5e-4, transient: 300, blockTime: 300, nBlocks: 8, every: 10, seed: 1, fit: [3, 11] }, + { N: 23, nu: 1e-8, dt: 1.6e-4, transient: 300, blockTime: 150, nBlocks: 8, every: 20, seed: 1, fit: [3, 13] }, + ], + }; + const k41 = p => p / 3; + // She & Leveque, Phys. Rev. Lett. 72, 336 (1994), for Navier-Stokes turbulence. + const sheLeveque = p => p / 9 + 2 * (1 - Math.pow(2 / 3, p / 3)); + + function create({ N = 22, nu = 1e-7, a = SABRA.a, b = SABRA.b, c = SABRA.c, lambda = SABRA.lambda, k0 = SABRA.k0, force = FORCE } = {}) { + const L = N + 4, k = new Float64Array(L); // index j = n + 2; ghosts j = 0, 1, N + 2, N + 3 + let kn = k0 / (lambda * lambda); + for (let j = 0; j < L; j++) { k[j] = kn; kn *= lambda; } + const A = new Float64Array(L), B = new Float64Array(L), Cc = new Float64Array(L), fr = new Float64Array(L), fi = new Float64Array(L); + for (let j = 2; j < N + 2; j++) { A[j] = a * k[j + 1]; B[j] = b * k[j]; Cc[j] = -c * k[j - 1]; } + for (const f of force) { fr[f.n + 2] = f.re; fi[f.n + 2] = f.im; } + const arr = () => new Float64Array(L); + return { + N, nu, a, b, c, lambda, k0, force, k, A, B, Cc, fr, fi, + ur: arr(), ui: arr(), k1r: arr(), k1i: arr(), k2r: arr(), k2i: arr(), k3r: arr(), k3i: arr(), k4r: arr(), k4i: arr(), + vr: arr(), vi: arr(), E1: arr(), E2: arr(), dt: 0, t: 0, + }; + } + // Initial state: |u_n| = amp for n < 2, then x0.8 per shell (about K41's 2^(-1/3)) up to nCut, then x0.25 per shell, + // with seeded random phases. + function init(m, { seed = 1, amp = 0.05, nCut = 8 } = {}) { + const rng = mulberry32(seed); + let r = amp; + for (let n = 0; n < m.N; n++) { + const [x, y] = unitPhase(rng); + m.ur[n + 2] = r * x; m.ui[n + 2] = r * y; + if (n >= 1) r *= n < nCut ? 0.8 : 0.25; + } + m.t = 0; + } + function setDt(m, dt) { + m.dt = dt; + for (let j = 0; j < m.N + 4; j++) { const g = -m.nu * m.k[j] * m.k[j] * dt; m.E1[j] = expDet(g); m.E2[j] = expDet(g / 2); } + } + // Nonlinear term plus forcing: out = i (...) + f. + function rhs(m, ur, ui, outr, outi) { + const A = m.A, B = m.B, Cc = m.Cc, fr = m.fr, fi = m.fi, top = m.N + 2; + for (let j = 2; j < top; j++) { + const r2 = ur[j + 2], i2 = ui[j + 2], r1 = ur[j + 1], i1 = ui[j + 1], rm = ur[j - 1], im = ui[j - 1], rmm = ur[j - 2], imm = ui[j - 2]; + const x = A[j] * (r2 * r1 + i2 * i1) + B[j] * (r1 * rm + i1 * im) + Cc[j] * (rm * rmm - im * imm); + const y = A[j] * (i2 * r1 - r2 * i1) + B[j] * (i1 * rm - r1 * im) + Cc[j] * (rm * imm + im * rmm); + outr[j] = fr[j] - y; outi[j] = fi[j] + x; + } + } + // nSteps of integrating-factor RK4 (Lawson type). + function advance(m, nSteps) { + const { ur, ui, k1r, k1i, k2r, k2i, k3r, k3i, k4r, k4i, vr, vi, E1, E2 } = m, h = m.dt, h2 = h / 2, h6 = h / 6, top = m.N + 2; + for (let s = 0; s < nSteps; s++) { + rhs(m, ur, ui, k1r, k1i); + for (let j = 2; j < top; j++) { vr[j] = E2[j] * (ur[j] + h2 * k1r[j]); vi[j] = E2[j] * (ui[j] + h2 * k1i[j]); } + rhs(m, vr, vi, k2r, k2i); + for (let j = 2; j < top; j++) { vr[j] = E2[j] * ur[j] + h2 * k2r[j]; vi[j] = E2[j] * ui[j] + h2 * k2i[j]; } + rhs(m, vr, vi, k3r, k3i); + for (let j = 2; j < top; j++) { vr[j] = E1[j] * ur[j] + h * E2[j] * k3r[j]; vi[j] = E1[j] * ui[j] + h * E2[j] * k3i[j]; } + rhs(m, vr, vi, k4r, k4i); + for (let j = 2; j < top; j++) { + ur[j] = E1[j] * ur[j] + h6 * (E1[j] * k1r[j] + 2 * E2[j] * (k2r[j] + k3r[j]) + k4r[j]); + ui[j] = E1[j] * ui[j] + h6 * (E1[j] * k1i[j] + 2 * E2[j] * (k2i[j] + k3i[j]) + k4i[j]); + } + } + m.t += nSteps * h; + } + // Im T_m for the triad starting at shell m (index j = m + 2): Im(u*_m u*_{m+1} u_{m+2}). + function imT(ur, ui, j) { + const pr = ur[j] * ur[j + 1] - ui[j] * ui[j + 1], pi = ur[j] * ui[j + 1] + ui[j] * ur[j + 1]; + return pr * ui[j + 2] - pi * ur[j + 2]; + } + // E = sum |u|^2 and sum q^n |u|^2 for each q in qs. + function invariants(m, qs = [m.a / m.c]) { + let E = 0; const H = qs.map(() => 0), w = qs.map(() => 1); + for (let n = 0; n < m.N; n++) { + const e = m.ur[n + 2] * m.ur[n + 2] + m.ui[n + 2] * m.ui[n + 2]; E += e; + for (let i = 0; i < qs.length; i++) { H[i] += w[i] * e; w[i] *= qs[i]; } + } + return { E, H }; + } + + const PMAX = 8; + // Statistics in nBlocks consecutive time blocks (batch means). Per block and shell: <|u_n|^p> and + // <|Pi_n/k_n|^(p/3)>, p = 1..8; ; . Per block: injection 2 Re, dissipation 2 nu sum k^2 |u|^2. + function newStats(N, nBlocks) { + return { N, nBlocks, cnt: new Float64Array(nBlocks), S: new Float64Array(nBlocks * PMAX * N), F: new Float64Array(nBlocks * PMAX * N), + Pi: new Float64Array(nBlocks * N), Sig: new Float64Array(nBlocks * N), inj: new Float64Array(nBlocks), dis: new Float64Array(nBlocks), + E: new Float64Array(nBlocks), maxRate: 0 }; + } + function sample(m, st, blk) { + const N = m.N, ur = m.ur, ui = m.ui, k = m.k, o = blk * PMAX * N, op = blk * N; + let inj = 0, dis = 0, E = 0, rate = 0; + for (let n = 0; n < N; n++) { + const j = n + 2, e = ur[j] * ur[j] + ui[j] * ui[j], r = Math.sqrt(e); + let pw = r; for (let p = 0; p < PMAX; p++) { st.S[o + p * N + n] += pw; pw *= r; } + const tn = imT(ur, ui, j), tm = imT(ur, ui, j - 1), P = 2 * (m.a * k[j + 1] * tn - m.c * k[j] * tm); + st.Pi[op + n] += P; st.Sig[op + n] += tm; + const g = Math.cbrt(Math.abs(P) / k[j]); let pf = g; for (let p = 0; p < PMAX; p++) { st.F[o + p * N + n] += pf; pf *= g; } + inj += 2 * (m.fr[j] * ur[j] + m.fi[j] * ui[j]); dis += 2 * m.nu * k[j] * k[j] * e; E += e; + if (k[j] * r > rate) rate = k[j] * r; + } + st.inj[blk] += inj; st.dis[blk] += dis; st.E[blk] += E; st.cnt[blk]++; + if (rate > st.maxRate) st.maxRate = rate; + } + // A full run: transient, then nBlocks blocks of blockTime, sampled every `every` steps. + function simulate({ N, nu, dt, seed = 1, transient, blockTime, nBlocks, every = 10, fit, force, a, b, c, amp, nCut } = {}) { + const m = create({ N, nu, force, a, b, c }); + init(m, { seed, amp, nCut }); + setDt(m, dt); + const tSteps = Math.round(transient / dt); + advance(m, tSteps); + const st = newStats(N, nBlocks), chunks = Math.max(1, Math.round(blockTime / dt / every)); + for (let blk = 0; blk < nBlocks; blk++) for (let i = 0; i < chunks; i++) { advance(m, every); sample(m, st, blk); } + return { config: { N, nu, dt, seed, transient, blockTime: chunks * every * dt, nBlocks, every }, fit, model: m, stats: st, + steps: tSteps + nBlocks * chunks * every, finite: m.ur.every(Number.isFinite) }; + } + + function linfit(x, y) { // least squares y = icept + slope x + const n = x.length; let sx = 0, sy = 0, sxx = 0, sxy = 0; + for (let i = 0; i < n; i++) { sx += x[i]; sy += y[i]; sxx += x[i] * x[i]; sxy += x[i] * y[i]; } + const slope = (n * sxy - sx * sy) / (n * sxx - sx * sx); + return { slope, icept: (sy - slope * sx) / n }; + } + const mean = v => v.reduce((s, x) => s + x, 0) / v.length; + const sd = v => { const mu = mean(v); return Math.sqrt(v.reduce((s, x) => s + (x - mu) * (x - mu), 0) / (v.length - 1)); }; + + // Averages over all blocks, or over a subset of block indices. + function averages(res, blocks) { + const st = res.stats, N = st.N, bl = blocks || [...Array(st.nBlocks).keys()]; + let cnt = 0; for (const b of bl) cnt += st.cnt[b]; + const S = [...Array(PMAX)].map(() => new Float64Array(N)), F = [...Array(PMAX)].map(() => new Float64Array(N)), Pi = new Float64Array(N), Sig = new Float64Array(N); + let inj = 0, dis = 0, E = 0; + for (const b of bl) { + for (let p = 0; p < PMAX; p++) for (let n = 0; n < N; n++) { S[p][n] += st.S[(b * PMAX + p) * N + n]; F[p][n] += st.F[(b * PMAX + p) * N + n]; } + for (let n = 0; n < N; n++) { Pi[n] += st.Pi[b * N + n]; Sig[n] += st.Sig[b * N + n]; } + inj += st.inj[b]; dis += st.dis[b]; E += st.E[b]; + } + for (let p = 0; p < PMAX; p++) for (let n = 0; n < N; n++) { S[p][n] /= cnt; F[p][n] /= cnt; } + for (let n = 0; n < N; n++) { Pi[n] /= cnt; Sig[n] /= cnt; } + return { S, F, Pi, Sig, inj: inj / cnt, dis: dis / cnt, E: E / cnt }; + } + // Where the cascade stops: the shell (interpolated in n) where first falls to frac * eps. + function fluxCut(av, eps, frac = 0.5) { + for (let n = 1; n < av.Pi.length - 1; n++) if (av.Pi[n] >= frac * eps && av.Pi[n + 1] < frac * eps) return n + (av.Pi[n] - frac * eps) / (av.Pi[n] - av.Pi[n + 1]); + return NaN; + } + // The shell where <|u_n|^2> falls to frac of its inertial power law fitted on [n1, n2]. + function spectralCut(av, n1, n2, frac = 0.5) { + const xs = [], ys = []; for (let n = n1; n <= n2; n++) { xs.push(n); ys.push(Math.log2(av.S[1][n])); } + const f = linfit(xs, ys), d = n => Math.log2(av.S[1][n]) - (f.icept + f.slope * n), lv = Math.log2(frac); + for (let n = n2; n < av.S[1].length - 1; n++) if (d(n) >= lv && d(n + 1) < lv) return n + (d(n) - lv) / (d(n) - d(n + 1)); + return NaN; + } + // Scaling exponents over shells [n1, n2] (log2 k_n = n, so slopes against n are exponents in k): + // u: zeta_p = -slope of log2 <|u_n|^p>, the plain structure functions; + // F: -slope of log2 F_p(n), F_p = <|Pi_n/k_n|^(p/3)>, the flux-based structure functions; + // essF: extended self-similarity, the slope of log F_p against log F_3 (zeta_3 = 1 by construction). + const KEYS = ['u', 'F', 'essF']; + function exponents(av, n1, n2) { + const ns = [], f3 = []; for (let n = n1; n <= n2; n++) { ns.push(n); f3.push(Math.log2(av.F[2][n])); } + const out = { u: [], F: [], essF: [] }; + for (let p = 0; p < PMAX; p++) { + const y = ns.map(n => Math.log2(av.S[p][n])), yF = ns.map(n => Math.log2(av.F[p][n])); + out.u.push(-linfit(ns, y).slope); out.F.push(-linfit(ns, yF).slope); out.essF.push(linfit(f3, yF).slope); + } + return out; + } + // Analysis with error bars: batch means over blocks, and the spread over fit ranges [n1 + d1, n2 + d2], d1, d2 in {-1, 0, 1}. + function analyse(res, [n1, n2] = res.fit || [3, 8]) { + const m = res.model, all = averages(res), eps = all.inj, nb = res.stats.nBlocks; + const z = exponents(all, n1, n2); + const perBlock = [...Array(nb).keys()].map(b => exponents(averages(res, [b]), n1, n2)); + const be = key => [...Array(PMAX).keys()].map(p => sd(perBlock.map(e => e[key][p])) / Math.sqrt(nb)); + const ranges = []; for (const d1 of [-1, 0, 1]) for (const d2 of [-1, 0, 1]) if (n2 + d2 - (n1 + d1) >= 3) ranges.push([n1 + d1, n2 + d2]); + const rs = ranges.map(([a, b]) => exponents(all, a, b)); + const spread = key => [...Array(PMAX).keys()].map(p => { const v = rs.map(e => e[key][p]); return (Math.max(...v) - Math.min(...v)) / 2; }); + const injBlocks = [...Array(nb).keys()].map(b => res.stats.inj[b] / res.stats.cnt[b]); + const ndBlocks = [...Array(nb).keys()].map(b => fluxCut(averages(res, [b]), injBlocks[b])); + const nu = m.nu, kEta = Math.pow(eps / (nu * nu * nu), 0.25), U = Math.sqrt(all.E), L = 1 / m.k0; + const plateau = mean([1, 2, 3, 4, 5].map(n => all.Pi[n])); + return { + eps, epsErr: sd(injBlocks) / Math.sqrt(nb), dis: all.dis, E: all.E, U, L, Re: U * L / nu, kEta, nEta: Math.log2(kEta / m.k0), plateau, + nd: fluxCut(all, eps), ndErr: sd(ndBlocks) / Math.sqrt(nb), ndSpec: spectralCut(all, n1, n2), n1, n2, av: all, zeta: z, + err: Object.fromEntries(KEYS.map(key => [key, be(key)])), spread: Object.fromEntries(KEYS.map(key => [key, spread(key)])), + }; + } + // Inviscid, unforced: largest relative change of E and of each sum q^n |u|^2 over time T, checked every `every` steps, + // each normalised by sum |q|^n |u_n|^2 at t = 0 (H itself can be near 0). Stops once E has changed by more than stopAt. + function conservation({ N = 10, a = 1, b = -0.5, c = -0.5, dt = 1e-3, T = 10, seed = 7, amp = 0.05, qs, every = 100, stopAt = 1 } = {}) { + const m = create({ N, nu: 0, a, b, c, force: [] }); + init(m, { seed, amp, nCut: N }); + setDt(m, dt); + const q = qs || [a / c], i0 = invariants(m, q), norm = invariants(m, q.map(Math.abs)).H; + let dE = 0, tStop = null; const dH = q.map(() => 0); + for (let i = 0, nChunks = Math.round(T / dt / every); i < nChunks; i++) { + advance(m, every); + const iv = invariants(m, q), e = (iv.E - i0.E) / i0.E; + if (!(Math.abs(e) <= Math.abs(dE))) dE = e; + iv.H.forEach((h, j) => { const d = (h - i0.H[j]) / norm[j]; if (!(Math.abs(d) <= Math.abs(dH[j]))) dH[j] = d; }); + if (!(Math.abs(e) <= stopAt)) { tStop = m.t; break; } + } + return { dE, dH, tStop }; + } + + // ---------- 3. Plasma scales (solar wind) ---------- + // A typical solar wind at 1 AU: density, speed, field, proton temperature (slow-wind protons are cooler, a few 1e4 K). + const WIND = { n: 5e6, V: 4e5, B: 5e-9, T: 1e5 }; + const OMEGA_SUN = 2 * Math.PI / (25.38 * 86400); // sidereal rotation of the solar equator, rad/s + const inertialLength = n => C.c / Math.sqrt(n * C.e * C.e / (C.eps0 * C.mp)); // c / omega_pi, m + // Proton gyroradius v_th / Omega_p; thermal speed sqrt(kT/m) (convention 1) or sqrt(2kT/m) (convention 2). + const gyroradius = (T, B, conv = 1) => Math.sqrt(conv * C.kB * T / C.mp) / (C.e * B / C.mp); + const betaP = ({ n, T, B }) => 2 * C.mu0 * n * C.kB * T / (B * B); + // A structure of size ell swept past a spacecraft at speed V. Fluctuations elongated along the field (k_perp >> k_par) + // are seen at V sin(theta_BV) / (2 pi ell), theta_BV the angle between the field and the flow. + const taylorFrequency = (V, ell, theta = Math.PI / 2) => V * Math.sin(theta) / (2 * Math.PI * ell); + // The Parker spiral: the field winds up at tan(psi) = Omega_sun r / V. + const parkerAngle = (V, r = C.AU) => Math.atan(OMEGA_SUN * r / V); + // Proton collision time (Braginskii, Rev. Plasma Phys. 1, 205 (1965), in SI) and mean free path sqrt(2kT/m) tau. + function protonCollisions({ n, T }, lnLambda = 25) { + const tau = 3 * Math.sqrt(C.mp) * Math.pow(C.kB * T, 1.5) * Math.pow(4 * Math.PI * C.eps0, 2) / (4 * Math.sqrt(Math.PI) * n * Math.pow(C.e, 4) * lnLambda); + return { tau, mfp: Math.sqrt(2 * C.kB * T / C.mp) * tau }; + } + + // ---------- 4. Molecular clouds ---------- + const soundSpeed = (T, mu = 2.33) => Math.sqrt(C.kB * T / (mu * C.mH)); // isothermal, m/s + // Linewidth-size relations, sigma = coef (size/pc)^exp km/s, and the thermal dispersion each sigma already holds + // (km/s, in the relation's own dimensions). Larson (1981): sigma is the total 3D rms velocity of all internal + // motions, thermal ones included, with sigma_s = 0.32 km/s for 10 K molecular gas (pp. 813-814; scaled here as + // sqrt(T) at other temperatures), and L the maximum projected size. Solomon et al. (1987): sigma is the 1D intensity-weighted dispersion of the 12CO line, so it holds + // CO's own thermal width, sqrt(kT/28 m_H); the fit is (1.0 +- 0.1) S^(0.5 +- 0.05) in the size parameter S, and the + // effective radius is R = 3.4 S/sqrt(pi) = 1.92 S (p. 737), so in R it is 0.722 R^0.5 (this study's conversion; the + // paper does not print it). All as read from the papers in #34's review. + const CO_THERMAL = T => Math.sqrt(C.kB * T / (28 * C.mH)) / 1e3; + const LARSON81 = { coef: 1.10, exp: 0.38, dims: 3, thermal: T => 0.32 * Math.sqrt(T / 10), size: 'maximum projected size' }; + // Objects in Larson's Table 1 (pp. 812-813): size (pc) and sigma_3D (km/s). On p. 820 he names TMC1 and TMC2 as the + // only ones whose internal motions are subsonic (as read in #34's re-review). + const LARSON_OBJECTS = [ + { name: 'TMC2', size: 0.1, sigma: 0.35, subsonic: true }, { name: 'L1407/L1257', size: 0.12, sigma: 0.86, subsonic: false }, + { name: 'L1544/L63', size: 0.18, sigma: 0.65, subsonic: false }, { name: 'TMC1', size: 0.2, sigma: 0.40, subsonic: true }, + ]; + const RE_S = 3.4 / Math.sqrt(Math.PI); + const solomon = (coefS = 1.0, exp = 0.5) => ({ coef: coefS / Math.pow(RE_S, exp), exp, dims: 1, thermal: CO_THERMAL, size: 'effective radius' }); + const SOLOMON87 = solomon(); + // The size (pc, in the relation's own convention) where the 1D non-thermal dispersion equals the sound speed: + // sigma^2 = dims c_s^2 + thermal^2 there. + const sonicScale = (rel, T = 10) => { + const cs = soundSpeed(T) / 1e3, th = rel.thermal(T); + return Math.pow(Math.sqrt(rel.dims * cs * cs + th * th) / rel.coef, 1 / rel.exp); + }; + // The non-thermal Mach number of a dispersion sigma (km/s) read in a relation's convention: the 1D non-thermal + // dispersion over c_s. machAt: the same for the relation's own sigma at a given size (pc). + const machOf = (rel, sigma, T = 10) => { const th = rel.thermal(T); return Math.sqrt((sigma * sigma - th * th) / rel.dims) / (soundSpeed(T) / 1e3); }; + const machAt = (rel, size, T = 10) => machOf(rel, rel.coef * Math.pow(size, rel.exp), T); + + const api = { + C, NU, kolmogorov, outerOverEta, + SABRA, FORCE, RUNS, PMAX, KEYS, k41, sheLeveque, mulberry32, unitPhase, expDet, create, init, setDt, rhs, advance, imT, invariants, + newStats, sample, simulate, linfit, mean, sd, averages, fluxCut, spectralCut, exponents, analyse, conservation, + WIND, inertialLength, gyroradius, betaP, taylorFrequency, parkerAngle, protonCollisions, + soundSpeed, CO_THERMAL, LARSON81, LARSON_OBJECTS, SOLOMON87, RE_S, solomon, sonicScale, machOf, machAt, + }; + if (typeof module !== 'undefined' && module.exports) module.exports = api; else root.TURB = api; +})(typeof window !== 'undefined' ? window : globalThis); diff --git a/turbulence/seeds.mjs b/turbulence/seeds.mjs new file mode 100644 index 0000000..a42292c --- /dev/null +++ b/turbulence/seeds.mjs @@ -0,0 +1,46 @@ +// The spread over random seeds behind the README's seed figures: the test's run over 32 seeds (the 4/5 law's +// deviation and the exponents), and the viscosity scan over 8 seeds (the slope of the cascade's width against Re). +// Usage: node turbulence/seeds.mjs ci (32 seeds of the test's run, about a minute) +// node turbulence/seeds.mjs scan (8 seeds of the four-viscosity scan, about 3 minutes) +// node turbulence/seeds.mjs (both) +import { createRequire } from 'node:module'; +const T = createRequire(import.meta.url)('./model.js'); +const mode = process.argv[2] || 'both', f3 = x => x.toFixed(3), f4 = x => x.toFixed(4); +const mean = v => v.reduce((a, x) => a + x, 0) / v.length, sd = v => Math.sqrt(v.reduce((a, x) => a + (x - mean(v)) ** 2, 0) / (v.length - 1)); +const median = v => { const s = [...v].sort((a, b) => a - b), m = s.length >> 1; return s.length % 2 ? s[m] : (s[m - 1] + s[m]) / 2; }; + +if (mode === 'ci' || mode === 'both') { + const rows = []; + for (let seed = 1; seed <= 32; seed++) { + const res = T.simulate({ ...T.RUNS.CI, seed }), a = T.analyse(res), k = res.model.k, av = a.av; + rows.push({ law: Math.max(...[1, 2, 3, 4, 5].map(n => Math.abs(3 * k[n + 2] * av.Sig[n] / a.plateau - (1 - Math.pow(-0.5, n))))), + z2: a.zeta.essF[1], z4: a.zeta.essF[3], z6: a.zeta.essF[5], gap: a.nd - a.nEta }); + } + const law = rows.map(r => r.law); + console.log(`The test's run over 32 seeds: largest 4/5-law deviation in shells 1-5, median ${f4(median(law))}, maximum ${f4(Math.max(...law))}; ` + + `${law.filter(x => x > 0.005).length} of 32 above 0.005, ${law.filter(x => x > 0.01).length} above 0.01.`); + for (const k of ['z2', 'z4', 'z6']) console.log(` zeta_${k[1]}: mean ${f4(mean(rows.map(r => r[k])))}, SD ${f4(sd(rows.map(r => r[k])))}`); + console.log(` flux half-point minus k_eta: ${f3(Math.min(...rows.map(r => r.gap)))} to ${f3(Math.max(...rows.map(r => r.gap)))} shells`); +} + +if (mode === 'scan' || mode === 'both') { + const cuts = { 'flux 75%': [], 'flux 50%': [], 'flux 25%': [], 'spectrum 50%': [], 'spectrum 25%': [] }, weighted = []; + for (let seed = 1; seed <= 8; seed++) { + const xs = [], ys = {}, nd = [], ne = []; + for (const k of Object.keys(cuts)) ys[k] = []; + for (const cfg of T.RUNS.SCAN) { + const q = T.analyse(T.simulate({ ...cfg, seed })), dec = n => n * Math.log10(2); + xs.push(Math.log10(q.Re)); nd.push(dec(q.nd)); ne.push(dec(q.ndErr)); + for (const f of [0.75, 0.5, 0.25]) ys[`flux ${f * 100}%`].push(dec(T.fluxCut(q.av, q.eps, f))); + for (const f of [0.5, 0.25]) ys[`spectrum ${f * 100}%`].push(dec(T.spectralCut(q.av, q.n1, q.n2, f))); + } + for (const k of Object.keys(cuts)) cuts[k].push(T.linfit(xs, ys[k]).slope); + const w = ne.map(e => 1 / (e * e)), W = mean(w) * w.length, mx = xs.reduce((a, x, i) => a + w[i] * x, 0) / W, my = nd.reduce((a, y, i) => a + w[i] * y, 0) / W; + weighted.push(xs.reduce((a, x, i) => a + w[i] * (x - mx) * (nd[i] - my), 0) / xs.reduce((a, x, i) => a + w[i] * (x - mx) ** 2, 0)); + } + console.log('\nThe viscosity scan (every viscosity at 8 blocks) over 8 seeds: width against log10 Re, by cut'); + for (const [k, v] of Object.entries(cuts)) console.log(` ${k.padEnd(13)} ${v.map(f3).join(' ')} mean ${f3(mean(v))}, SD ${f3(sd(v))}`); + const all = Object.values(cuts).flat(); + console.log(` all cuts and seeds: ${f3(Math.min(...all))} to ${f3(Math.max(...all))}; flux half-point weighted by its errors: mean ${f3(mean(weighted))}, SD ${f3(sd(weighted))}`); + console.log(' (The spectrum cut at three quarters falls near the end of the fit range at nu = 1e-5 and is left out.)'); +} diff --git a/turbulence/table.mjs b/turbulence/table.mjs new file mode 100644 index 0000000..df9c3ed --- /dev/null +++ b/turbulence/table.mjs @@ -0,0 +1,101 @@ +// Prints where turbulence repeats itself and where it stops: the Kolmogorov scale in air and water, the Sabra cascade +// (flux, exact third-order law, scaling exponents, the dissipation cut against Reynolds number), the solar wind's ion +// scales and the molecular-cloud sonic scale. +// Usage: node turbulence/table.mjs (about a minute: the long run and the viscosity scan) +// node turbulence/table.mjs quick (a few seconds: the CI run instead, no scan) +import { createRequire } from 'node:module'; +const T = createRequire(import.meta.url)('./model.js'); +const quick = process.argv[2] === 'quick'; +const e2 = x => x.toExponential(2), f3 = x => x.toFixed(3); + +console.log('Kolmogorov scale eta = (nu^3/eps)^(1/4)\n'); +console.log('medium'.padEnd(10) + 'eps W/kg'.padStart(10) + 'eta mm'.padStart(9) + ' where'); +for (const [med, eps, where] of [['air', 0.4, 'surface layer, cold front'], ['air', 0.05, 'surface layer before the front'], ['air', 3.4e-3, '(eta = 1 mm)'], + ['air', 1e-7, 'quiet nocturnal boundary layer'], ['seawater', 1e-6, '(eta = 1 mm)'], ['seawater', 1e-7, 'upper mixed layer, by day'], ['seawater', 1e-8, 'upper mixed layer, quieter']]) { + console.log(med.padEnd(10) + e2(eps).padStart(10) + (T.kolmogorov(T.NU[med], eps).eta * 1e3).toFixed(2).padStart(9) + ' ' + where); +} + +console.log('\nHow many decades the repetition spans\n'); +for (const [name, ratio] of [['Wind tunnel, R_lambda = 6000: L/eta (Pope\'s conventions)', T.outerOverEta(6000)], + ['Atmospheric surface layer, R_lambda 1,600: L/eta', T.outerOverEta(1600)], ['the same, R_lambda 6,600', T.outerOverEta(6600)], + ['Interstellar electron density, 2e6 to 1e13 m (measured)', 1e13 / 2e6], ['the same, extended to 1e17 m (Kolmogorov assumed)', 1e17 / 2e6]]) { + console.log(`${name.padEnd(62)} ${Math.log10(ratio).toFixed(1)} decades`); +} + +const R = quick ? T.RUNS.CI : T.RUNS.TABLE; +const t0 = Date.now(), res = T.simulate(R), a = T.analyse(res), av = a.av, k = res.model.k, c = res.config; +console.log(`\nSabra shell model (${quick ? 'CI' : 'long'} run): N = ${c.N} shells, nu = ${e2(c.nu)}, dt = ${c.dt}, force 0.005(1 + i) on shell 0, ` + + `${c.nBlocks} blocks x ${c.blockTime} time units after ${c.transient} (${((Date.now() - t0) / 1e3).toFixed(1)} s)`); +console.log(`eps = ${e2(a.eps)} +- ${e2(a.epsErr)}, dissipation ${e2(a.dis)}, U = ${f3(a.U)}, large-eddy time L/U = ${(a.L / a.U).toFixed(1)}, ` + + `Re = U L/nu = ${e2(a.Re)}, k_eta = (eps/nu^3)^(1/4) at shell ${a.nEta.toFixed(2)}\n`); +console.log(' n Pi_n/eps 3 k_n Sigma_n/eps exact 1 - (-1/2)^n log2 <|u_n|^2>'); +for (let n = 0; n < c.N - 1; n++) { + console.log(String(n).padStart(3) + f3(av.Pi[n] / a.eps).padStart(11) + (3 * k[n + 2] * av.Sig[n] / a.eps).toFixed(4).padStart(20) + + (1 - Math.pow(-0.5, n)).toFixed(4).padStart(21) + Math.log2(av.S[1][n]).toFixed(2).padStart(17)); +} +console.log(`\nScaling exponents over shells ${a.n1}..${a.n2}: value +- batch-means error [half-range over fit ranges n1 +- 1, n2 +- 1]`); +console.log('u: plain <|u_n|^p>; F: flux-based <|Pi_n/k_n|^(p/3)>; essF: F against F_3 (extended self-similarity)\n'); +console.log(' p' + ['u', 'F', 'essF'].map(s => s.padStart(26)).join('') + 'K41'.padStart(8) + 'She-Leveque'.padStart(13)); +for (let p = 0; p < T.PMAX; p++) { + console.log(String(p + 1).padStart(3) + T.KEYS.map(s => `${f3(a.zeta[s][p])} +- ${a.err[s][p].toFixed(4)} [${f3(a.spread[s][p])}]`.padStart(26)).join('') + + f3(T.k41(p + 1)).padStart(8) + f3(T.sheLeveque(p + 1)).padStart(13)); +} +console.log("L'vov et al. 1998, Table I (34 shells, forcing on shells 0-1): <|u|^q> zeta_2,4,6 = 0.720 +- 0.008, 1.256 +- 0.012, 1.706 +- 0.015;"); +console.log(' flux-based, q = 2, 4, 5, 6, 7: 0.719, 1.249, 1.477, 1.691, 1.893 (as read from the paper in #34\'s review).'); +console.log('Measured in a wind tunnel at R_lambda 2000-6000: zeta_2 = 0.693 +- 0.011 (Kuchler, Bodenschatz & Bewley 2020, arXiv:2006.10993v2,'); +console.log(' Fig. 3; the abstract still gives +- 0.003); the same group\'s later analysis gives 0.698 +- 0.011 (arXiv:2106.01795v2, Fig. 5).'); +console.log(`Velocity across a size l grows as l^(zeta_2/2): ${f3(a.zeta.essF[1] / 2)} here, 1/3 in K41, ${f3(T.sheLeveque(2) / 2)} for She-Leveque, 1/2 for shocks.`); +const fl = a.av.S[3].map((s4, n) => s4 / (a.av.S[1][n] * a.av.S[1][n])), hStar = (a.zeta.essF[3] - a.zeta.essF[1]) / 2; +console.log(`Flatness <|u|^4>/<|u|^2>^2 grows toward small scales as k^(2 zeta_2 - zeta_4) = k^${f3(2 * a.zeta.essF[1] - a.zeta.essF[3])}: ` + + `${[2, 4, 6, 8, 10].map(n => `${fl[n].toFixed(2)} at shell ${n}`).join(', ')}.`); + +if (!quick) { + console.log('\nWhere the cascade stops, against Reynolds number (each run 8 blocks; the nu = 1e-7 row is the long run)\n'); + console.log(' nu N Re n_eta flux half-point n_d - n_eta spectrum half-point width, decades (3/4) log10 Re'); + const xs = [], cuts = { 'flux 75%': [], 'flux 50%': [], 'flux 25%': [], 'spectrum 50%': [], 'spectrum 25%': [] }; + const nds = [], ndErrs = []; + for (const cfg of T.RUNS.SCAN) { + const r = cfg.nu === R.nu ? { res, a } : (() => { const rr = T.simulate(cfg); return { res: rr, a: T.analyse(rr) }; })(), q = r.a; + const width = q.nd * Math.log10(2), dec = n => n * Math.log10(2); + xs.push(Math.log10(q.Re)); nds.push(width); ndErrs.push(q.ndErr * Math.log10(2)); + for (const f of [0.75, 0.5, 0.25]) cuts[`flux ${f * 100}%`].push(dec(T.fluxCut(q.av, q.eps, f))); + for (const f of [0.5, 0.25]) cuts[`spectrum ${f * 100}%`].push(dec(T.spectralCut(q.av, q.n1, q.n2, f))); + console.log(e2(cfg.nu).padStart(7) + String(cfg.N).padStart(4) + e2(q.Re).padStart(11) + q.nEta.toFixed(2).padStart(8) + `${q.nd.toFixed(2)} +- ${q.ndErr.toFixed(2)}`.padStart(18) + + (q.nd - q.nEta).toFixed(2).padStart(14) + q.ndSpec.toFixed(2).padStart(22) + width.toFixed(2).padStart(17) + (0.75 * Math.log10(q.Re)).toFixed(2).padStart(17)); + } + console.log('Width in decades (forcing shell to the cut) against log10 Re, slope by where the cut is put (where the flux or the'); + console.log('spectrum, relative to its power law, has fallen to the given fraction):'); + console.log(' ' + Object.entries(cuts).map(([name, ys]) => `${name} ${f3(T.linfit(xs, ys).slope)}`).join(', ')); + { // the flux half-points weighted by their batch-means errors + const w = ndErrs.map(e => 1 / (e * e)), W = w.reduce((a, x) => a + x, 0), mx = xs.reduce((a, x, i) => a + w[i] * x, 0) / W, my = nds.reduce((a, y, i) => a + w[i] * y, 0) / W; + const sxx = xs.reduce((a, x, i) => a + w[i] * (x - mx) ** 2, 0), b = xs.reduce((a, x, i) => a + w[i] * (x - mx) * (nds[i] - my), 0) / sxx; + console.log(` flux half-point weighted by its errors: ${f3(b)} +- ${f3(1 / Math.sqrt(sxx))} (seed 1; over 8 seeds this fit has an SD of about 0.002: node turbulence/seeds.mjs)`); + } + console.log(`K41: 3/4. Viscous loss at the cut, nu k_d^2 S_2(k_d) = eps with S_2 ~ k^-zeta_2: 1/(2 - zeta_2) = ${f3(1 / (2 - a.zeta.essF[1]))}.`); + console.log(`Multifractal cut-off Re^(1/(1 + h*)), h* = dzeta_p/dp at p = 3 ~ (zeta_4 - zeta_2)/2: ${f3(1 / (1 + hStar))}. A local Reynolds number`); + console.log(`of 1 on the rms velocity at the cut, S_2(k_d)^(1/2) = nu k_d, would give 2/(2 + zeta_2) = ${f3(2 / (2 + a.zeta.essF[1]))}, below 3/4.`); +} + +const W = T.WIND, di = T.inertialLength(W.n), col = T.protonCollisions(W), cool = T.protonCollisions({ ...W, T: 4e4 }), psi = T.parkerAngle(W.V); +console.log(`\nSolar wind at 1 AU (n = ${W.n / 1e6} cm^-3, V = ${W.V / 1e3} km/s, B = ${W.B * 1e9} nT, T = ${e2(W.T)} K, beta_p = ${T.betaP(W).toFixed(2)};` + + ` Parker-spiral angle ${(psi * 180 / Math.PI).toFixed(0)} deg)\n`); +console.log('length'.padEnd(40) + 'km'.padStart(6) + ' swept past at V: no angle, V sin(theta_BV)'); +for (const [name, ell] of [['proton inertial length c/omega_pi', di], ['proton gyroradius, v_th = sqrt(2kT/m)', T.gyroradius(W.T, W.B, 2)], + ['proton gyroradius, v_th = sqrt(kT/m)', T.gyroradius(W.T, W.B, 1)]]) { + console.log(`${name.padEnd(40)} ${(ell / 1e3).toFixed(0).padStart(5)} ${T.taylorFrequency(W.V, ell).toFixed(2).padStart(10)} Hz ${T.taylorFrequency(W.V, ell, psi).toFixed(2).padStart(10)} Hz`); +} +console.log(`proton collision time ${(col.tau / 86400).toFixed(0)} days, mean free path ${(col.mfp / T.C.AU).toFixed(2)} AU = ${e2(col.mfp / di)} inertial lengths;` + + ` at 4e4 K ${(cool.tau / 86400).toFixed(0)} days, ${(cool.mfp / T.C.AU).toFixed(2)} AU (ln Lambda = 25)`); +console.log('Observed: breaks between 0.1 and 1 Hz across conditions, at the larger ion scale (Chen et al. 2014). Sahraoui et al. (2009) find'); +console.log('0.4 Hz in one fast, high-beta interval and read it as the proton gyroscale; Leamon et al.\'s 33 intervals, as replotted in 2000,'); +console.log('run from about 0.12 to 0.87 Hz, mean near 0.38 Hz. Both estimates above lie inside the spread; the data do not choose.'); + +console.log(`\nMolecular clouds at 10 K: sound speed ${(T.soundSpeed(10) / 1e3).toFixed(3)} km/s (mu = 2.33). Where the 1D non-thermal dispersion falls to it:\n`); +console.log(`Larson (1981), sigma_3D = 1.10 L^0.38, which holds the gas's thermal sigma_s = 0.32 km/s: L = ${T.sonicScale(T.LARSON81).toFixed(3)} pc (maximum projected size);` + + ` non-thermal Mach ${[0.1, 1, 10, 100].map(l => T.machAt(T.LARSON81, l).toFixed(1)).join(', ')} at 0.1, 1, 10, 100 pc`); +console.log(` Larson's own criterion, sigma = sigma_s (p. 820): L = ${Math.pow(T.LARSON81.thermal(10) / T.LARSON81.coef, 1 / T.LARSON81.exp).toFixed(3)} pc`); +{ const ds = []; for (const c of [0.9, 1.1]) for (const e of [0.45, 0.55]) ds.push(2 * T.sonicScale(T.solomon(c, e))); + console.log(`Solomon et al. (1987), sigma_1D = 1.0 S^0.5 = ${T.SOLOMON87.coef.toFixed(3)} R^0.5, R = 1.92 S, CO's thermal width removed: R = ${T.sonicScale(T.SOLOMON87).toFixed(3)} pc` + + ` (effective radius), diameter ${(2 * T.sonicScale(T.SOLOMON87)).toFixed(3)} pc; ${Math.min(...ds).toFixed(2)} to ${Math.max(...ds).toFixed(2)} pc across within the fit's errors`); } +console.log('Larson\'s data reach about 0.1 pc; Solomon et al.\'s giant clouds are pc-sized and larger. Observed: filaments 0.10 +- 0.03 pc wide'); +console.log('(FWHM); coherent cores of about 0.1 pc. Below the sonic scale a shallower, subsonic cascade continues.');