rlx-lbm 0.2.18

Moment-encoded lattice Boltzmann (HOME-LBM) for RLX — D2Q9/D3Q27 reconstruction, collision, and periodic streaming as rlx graphs.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

// RLX — versatile ML compiler + runtime.
// Copyright (C) 2026 Eugene Hauptmann, Nataliya Kosmyna.
// SPDX-License-Identifier: MIT OR Apache-2.0

//! The moment round-trip oracle.
//!
//! A moment-encoded solver *is* its reconstruction: the populations are never
//! stored, so `reconstruct` defines the scheme. That leaves nothing to compare
//! it against — a CPU/GPU parity test compares two copies of the same formula,
//! and a regression test compares it to whatever it produced yesterday.
//!
//! The one thing that constrains it from outside is the definition of a moment:
//!
//! ```text
//! Σ f_i = ρ        Σ f_i c_α = ρ u_α        Σ f_i c_α c_β = ρ (S_αβ + cs² δ_αβ)
//! ```
//!
//! If `moments(reconstruct(m)) ≠ m`, the scheme is inconsistent no matter how
//! plausible its output looks. This matters more than it sounds: a mis-scaled
//! reconstruction does not blow up, it silently rescales the off-equilibrium
//! stress — which is the viscosity. The flow still looks like a flow. It is the
//! forward-path analogue of a finite-difference gradient check.
//!
//! # What it does not catch
//!
//! The Hermite **truncation order**. The third-order tensors are orthogonal to
//! the first three moments by construction — they live in the lattice's ghost
//! moment space — so a second-order reconstruction round-trips *exactly* while
//! changing the populations, and therefore the dynamics. Verified in
//! `round_trip_is_blind_to_truncation_order`.
//!
//! So the two oracles cover different things and neither subsumes the other:
//!
//! | oracle | constrains |
//! |---|---|
//! | this round trip | how the reconstruction is **assembled** — coefficients, scales, index layout |
//! | Taylor–Green decay vs `exp(−2νk²t)` | what the scheme **does** — effective viscosity, truncation order |
//!
//! The bug this was written against is squarely in the first column: applying
//! `1/(2cs⁴)` as 1.5 instead of 4.5 while double-counting the equilibrium terms
//! yields `S_recovered = u ⊗ u + S/3`, caught here at ~1e-3 against a 1e-12
//! threshold.

use crate::moment::{
    Moments2d, Moments3d, moments_d2q9, moments_d3q27, reconstruct_d2q9, reconstruct_d3q27,
};

/// How far a state drifts through one encode/decode cycle.
#[derive(Debug, Clone, Copy, Default)]
pub struct Residual {
    /// `|ρ' − ρ|`.
    pub rho: f64,
    /// `max |u'_α − u_α|`.
    pub u: f64,
    /// `max |S'_αβ − S_αβ|`.
    pub s: f64,
}

impl Residual {
    /// The largest component — one number to threshold on.
    pub fn max(&self) -> f64 {
        self.rho.max(self.u).max(self.s)
    }
    /// True when every component is within `tol`.
    pub fn within(&self, tol: f64) -> bool {
        self.max() <= tol
    }
}

/// Round-trip residual for a D2Q9 state.
pub fn round_trip_d2q9(m: &Moments2d) -> Residual {
    let back = moments_d2q9(&reconstruct_d2q9(m));
    Residual {
        rho: (back.rho - m.rho).abs(),
        u: (0..2)
            .map(|i| (back.u[i] - m.u[i]).abs())
            .fold(0.0, f64::max),
        s: (0..3)
            .map(|i| (back.s[i] - m.s[i]).abs())
            .fold(0.0, f64::max),
    }
}

/// Round-trip residual for a D3Q27 state.
pub fn round_trip_d3q27(m: &Moments3d) -> Residual {
    let back = moments_d3q27(&reconstruct_d3q27(m));
    Residual {
        rho: (back.rho - m.rho).abs(),
        u: (0..3)
            .map(|i| (back.u[i] - m.u[i]).abs())
            .fold(0.0, f64::max),
        s: (0..6)
            .map(|i| (back.s[i] - m.s[i]).abs())
            .fold(0.0, f64::max),
    }
}

/// Deterministic spread of physically reasonable D2Q9 states: equilibrium plus a
/// small off-equilibrium stress, at low Mach number.
pub fn sample_states_2d() -> Vec<Moments2d> {
    let mut v = Vec::new();
    for &rho in &[1.0f64, 0.85, 1.2] {
        for &(ux, uy) in &[(0.0, 0.0), (0.05, -0.03), (-0.08, 0.06), (0.1, 0.1)] {
            for &k in &[0.0f64, 1e-4, 1e-3, 5e-3] {
                v.push(Moments2d {
                    rho,
                    u: [ux, uy],
                    s: [ux * ux + k, uy * uy - 0.6 * k, ux * uy + 0.5 * k],
                });
            }
        }
    }
    v
}

/// Deterministic spread of D3Q27 states.
pub fn sample_states_3d() -> Vec<Moments3d> {
    let mut v = Vec::new();
    for &rho in &[1.0f64, 0.9, 1.15] {
        for &(ux, uy, uz) in &[
            (0.0, 0.0, 0.0),
            (0.05, -0.03, 0.02),
            (-0.07, 0.04, -0.06),
            (0.09, 0.09, -0.09),
        ] {
            for &k in &[0.0f64, 1e-4, 1e-3, 5e-3] {
                v.push(Moments3d {
                    rho,
                    u: [ux, uy, uz],
                    s: [
                        ux * ux + k,
                        uy * uy - 0.6 * k,
                        uz * uz - 0.4 * k,
                        ux * uy + 0.5 * k,
                        ux * uz - 0.2 * k,
                        uy * uz + 0.3 * k,
                    ],
                });
            }
        }
    }
    v
}