sva_samples/physics/
willemsen_bilbao_serafin.rs

1// Concern: one bowed-string call site's finite-difference state | Non-concern: argument evaluation, builtin dispatch | IO: (params, sr) -> a site; () -> a sample
2
3//! Elasto-plastic friction: Dupont/Hayward/Armstrong/Altpeter, IEEE TAC 47(5) (2002). First
4//! bowed use: Serafin/Avanzini/Rocchesso, SMAC-03. Corrected eqs. (7)-(9)/FD coupling:
5//! Willemsen/Bilbao/Serafin, DAFx-19 pp. 40-46. Reuses `chaigne_askenfelt`'s FD string.
6
7use crate::physics::Solver;
8
9use crate::physics::bound::Bound::*;
10use crate::physics::bound::all;
11use crate::physics::chaigne_askenfelt::{StringGrid, stencil_update, stiff_string_grid};
12
13#[derive(Clone, Debug, PartialEq)]
14pub struct WillemsenBilbaoSerafinParams {
15    pub f0: f64,
16    pub b: f64,
17    pub bow_pos: f64,
18    pub bow_vel: f64,
19    pub bow_force: f64,
20    pub mu_s: f64,
21    pub mu_c: f64,
22    pub stribeck_vel: f64,
23    pub bristle_stiffness: f64,
24    pub bristle_damping: f64,
25    pub viscous_friction: f64,
26    pub damp_dc: f64,
27    pub damp_freq: f64,
28}
29
30impl WillemsenBilbaoSerafinParams {
31    /// DAFx-19 Table 1's reference values, at the fundamental asked for.
32    pub fn at(f0: f64) -> WillemsenBilbaoSerafinParams {
33        WillemsenBilbaoSerafinParams {
34            f0,
35            b: 1e-4,
36            bow_pos: 0.25,
37            bow_vel: 0.1,
38            bow_force: 10.0,
39            mu_s: 0.8,
40            mu_c: 0.3,
41            stribeck_vel: 0.1,
42            bristle_stiffness: 1e4,
43            bristle_damping: 0.1,
44            viscous_friction: 0.4,
45            damp_dc: 1.0,
46            damp_freq: 5e-3,
47        }
48    }
49}
50
51impl WillemsenBilbaoSerafinParams {
52    pub fn valid(&self) -> bool {
53        all(&[
54            (self.f0, Positive),
55            (self.b, NonNegative),
56            (self.bow_pos, OpenUnit),
57            (self.bow_vel, Finite),
58            (self.bow_force, Positive),
59            (self.mu_s, Positive),
60            (self.mu_c, Positive),
61            (self.stribeck_vel, Positive),
62            (self.bristle_stiffness, Positive),
63            (self.bristle_damping, NonNegative),
64            (self.viscous_friction, NonNegative),
65            (self.damp_dc, NonNegative),
66            (self.damp_freq, NonNegative),
67        ])
68    }
69}
70
71const WIRE_DENSITY_KG_M3: f64 = 7850.0;
72const WIRE_RADIUS_M: f64 = 5e-4;
73/// Fixed per Table 1; `c = 2 f0 L` scales instead (`chaigne_askenfelt` inverts that).
74const WIRE_LENGTH_M: f64 = 1.0;
75const NR_MAX_ITERATIONS: usize = 50;
76const NR_TOLERANCE: f64 = 1e-7;
77
78/// `c` comes back with the grid so `tension` shares the derivation.
79fn build_grid(params: &WillemsenBilbaoSerafinParams, sr: f64) -> (StringGrid, f64) {
80    let rho = std::f64::consts::PI * WIRE_RADIUS_M * WIRE_RADIUS_M * WIRE_DENSITY_KG_M3;
81    let c = 2.0 * params.f0 * WIRE_LENGTH_M;
82    let kappa = c * WIRE_LENGTH_M * params.b.sqrt() / std::f64::consts::PI;
83    let grid = stiff_string_grid(
84        rho,
85        c,
86        WIRE_LENGTH_M,
87        kappa,
88        params.damp_dc,
89        params.damp_freq,
90        sr,
91    );
92    (grid, c)
93}
94
95fn sgn(x: f64) -> f64 {
96    if x > 0.0 {
97        1.0
98    } else if x < 0.0 {
99        -1.0
100    } else {
101        0.0
102    }
103}
104
105/// Eq. (7): `abs()` around `z_ss`, missing in the pre-DAFx-19 literature.
106fn steady_state(v: f64, s0: f64, f_c: f64, f_s: f64, stribeck_vel: f64) -> f64 {
107    sgn(v) / s0 * (f_c + (f_s - f_c) * (-(v / stribeck_vel).powi(2)).exp())
108}
109
110/// Eq. (8)-(9): `sgn(z)` on the sine offset; `|z|=z_ba`/`|z_ss|` resolved, not undefined.
111fn adhesion(v: f64, z: f64, z_ba: f64, zss: f64) -> f64 {
112    if sgn(v) != sgn(z) {
113        return 0.0;
114    }
115    let (az, azss) = (z.abs(), zss.abs());
116    if az <= z_ba {
117        0.0
118    } else if az >= azss {
119        1.0
120    } else {
121        let span = (azss - z_ba).max(1e-300);
122        let s = sgn(z);
123        0.5 * (1.0 + s * (std::f64::consts::PI * (z - s * 0.5 * (azss + z_ba)) / span).sin())
124    }
125}
126
127/// Eq. (6)/(15). `v == 0` avoids dividing by `z_ss(0) = 0`.
128fn bristle_rate(v: f64, z: f64, s0: f64, f_c: f64, f_s: f64, stribeck_vel: f64, z_ba: f64) -> f64 {
129    if v == 0.0 {
130        return 0.0;
131    }
132    let zss = steady_state(v, s0, f_c, f_s, stribeck_vel);
133    let a = adhesion(v, z, z_ba, zss);
134    v * (1.0 - a * z / zss)
135}
136
137/// Eq. (4), noise term `s3*w` dropped (disclosed v1 gap).
138fn friction_force(v: f64, z: f64, r: f64, s0: f64, s1: f64, s2: f64) -> f64 {
139    s0 * z + s1 * r + s2 * v
140}
141
142struct Coupling {
143    coeff: f64,
144    b_known: f64,
145    s0: f64,
146    s1: f64,
147    s2: f64,
148    f_c: f64,
149    f_s: f64,
150    stribeck_vel: f64,
151    z_ba: f64,
152    z_prev: f64,
153    r_prev: f64,
154    dt: f64,
155}
156
157impl Coupling {
158    /// `g1` re-derives eq. (17)-(19): theirs assume centred damping, the FD string's is
159    /// backward. `g2` keeps their eq. (20)-(21).
160    fn residual(&self, v: f64, z: f64) -> (f64, f64) {
161        let r = bristle_rate(
162            v,
163            z,
164            self.s0,
165            self.f_c,
166            self.f_s,
167            self.stribeck_vel,
168            self.z_ba,
169        );
170        let f = friction_force(v, z, r, self.s0, self.s1, self.s2);
171        let g1 = v + self.coeff * f - self.b_known;
172        let a = 2.0 * (z - self.z_prev) / self.dt - self.r_prev;
173        (g1, r - a)
174    }
175
176    /// A central-difference Jacobian, not Algorithm 1's analytic one.
177    fn newton_step(&self, v: f64, z: f64) -> Option<(f64, f64)> {
178        let (g1, g2) = self.residual(v, z);
179        let hv = v.abs().max(1e-3) * 1e-6;
180        let hz = z.abs().max(self.z_ba).max(1e-9) * 1e-6;
181        let (g1_vp, g2_vp) = self.residual(v + hv, z);
182        let (g1_vm, g2_vm) = self.residual(v - hv, z);
183        let (g1_zp, g2_zp) = self.residual(v, z + hz);
184        let (g1_zm, g2_zm) = self.residual(v, z - hz);
185        let dg1_dv = (g1_vp - g1_vm) / (2.0 * hv);
186        let dg2_dv = (g2_vp - g2_vm) / (2.0 * hv);
187        let dg1_dz = (g1_zp - g1_zm) / (2.0 * hz);
188        let dg2_dz = (g2_zp - g2_zm) / (2.0 * hz);
189        let det = dg1_dv * dg2_dz - dg1_dz * dg2_dv;
190        if !det.is_finite() || det.abs() < 1e-300 {
191            return None;
192        }
193        let dv = (g1 * dg2_dz - g2 * dg1_dz) / det;
194        let dz = (dg1_dv * g2 - dg2_dv * g1) / det;
195        let (new_v, new_z) = (v - dv, z - dz);
196        if !new_v.is_finite() || !new_z.is_finite() {
197            return None;
198        }
199        Some((new_v, new_z))
200    }
201}
202
203pub struct WillemsenBilbaoSerafinSite {
204    strings: Vec<StringGrid>,
205    contact_index: usize,
206    tension: f64,
207    dt: f64,
208    bow_vel: f64,
209    f_c: f64,
210    f_s: f64,
211    stribeck_vel: f64,
212    s0: f64,
213    s1: f64,
214    s2: f64,
215    z_ba: f64,
216    z: f64,
217    z_prev: f64,
218    /// `r^{n-1}`, eq. (20)-(21)'s extra trapezoidal state.
219    r_prev: f64,
220    last_v: f64,
221    last_f: f64,
222}
223
224impl WillemsenBilbaoSerafinSite {
225    pub fn new(params: &WillemsenBilbaoSerafinParams, sr: f64) -> WillemsenBilbaoSerafinSite {
226        let (grid, c) = build_grid(params, sr);
227        let contact_index = (params.bow_pos * grid.n as f64)
228            .round()
229            .clamp(1.0, (grid.n - 1) as f64) as usize;
230        let tension = c * c * grid.rho;
231        let f_c = params.mu_c * params.bow_force;
232        let f_s = params.mu_s * params.bow_force;
233        // Table 1: z_ba = 0.7*f_C/s0 — off the kinetic (mu_c) force, not the static one.
234        let z_ba = 0.7 * f_c / params.bristle_stiffness;
235        WillemsenBilbaoSerafinSite {
236            strings: vec![grid],
237            contact_index,
238            tension,
239            dt: 1.0 / sr,
240            bow_vel: params.bow_vel,
241            f_c,
242            f_s,
243            stribeck_vel: params.stribeck_vel,
244            s0: params.bristle_stiffness,
245            s1: params.bristle_damping,
246            s2: params.viscous_friction,
247            z_ba,
248            z: 0.0,
249            z_prev: 0.0,
250            r_prev: 0.0,
251            last_v: 0.0,
252            last_f: 0.0,
253        }
254    }
255}
256
257impl Solver for WillemsenBilbaoSerafinSite {
258    fn step(&mut self) -> f64 {
259        let dt = self.dt;
260        let l = self.contact_index;
261        let coupling_prev = (self.z, self.r_prev);
262
263        let grid = &mut self.strings[0];
264        let n = grid.n;
265        let mut free_next_l = 0.0;
266        for i in 1..n {
267            let next = stencil_update(grid, i, 0.0, 0.0);
268            if i == l {
269                free_next_l = next;
270            }
271            grid.y_next[i] = next;
272        }
273        grid.y_next[0] = 0.0;
274        grid.y_next[n] = 0.0;
275
276        let coupling = Coupling {
277            coeff: dt / (2.0 * grid.rho * grid.dx),
278            b_known: (free_next_l - grid.y_prev[l]) / (2.0 * dt) - self.bow_vel,
279            s0: self.s0,
280            s1: self.s1,
281            s2: self.s2,
282            f_c: self.f_c,
283            f_s: self.f_s,
284            stribeck_vel: self.stribeck_vel,
285            z_ba: self.z_ba,
286            z_prev: coupling_prev.0,
287            r_prev: coupling_prev.1,
288            dt,
289        };
290
291        let (mut v, mut z) = (self.last_v, self.z);
292        let mut converged = false;
293        for _ in 0..NR_MAX_ITERATIONS {
294            let Some((new_v, new_z)) = coupling.newton_step(v, z) else {
295                break;
296            };
297            let step_norm = ((new_v - v).powi(2) + (new_z - z).powi(2)).sqrt();
298            v = new_v;
299            z = new_z;
300            if step_norm < NR_TOLERANCE {
301                converged = true;
302                break;
303            }
304        }
305        debug_assert!(converged, "the friction solve never settled");
306
307        let r = bristle_rate(
308            v,
309            z,
310            self.s0,
311            self.f_c,
312            self.f_s,
313            self.stribeck_vel,
314            self.z_ba,
315        );
316        let f = friction_force(v, z, r, self.s0, self.s1, self.s2);
317        let grid = &mut self.strings[0];
318        grid.y_next[l] = free_next_l - (dt * dt / (grid.rho * grid.dx)) * f;
319
320        self.z_prev = self.z;
321        self.z = z;
322        self.r_prev = r;
323        self.last_v = v;
324        self.last_f = f;
325
326        let sample = self.tension * (grid.y_now[n] - grid.y_now[n - 1]) / grid.dx;
327
328        std::mem::swap(&mut grid.y_prev, &mut grid.y_now);
329        std::mem::swap(&mut grid.y_now, &mut grid.y_next);
330
331        sample
332    }
333}
sva_samples/physics/willemsen_bilbao_serafin.rs

sva_samples/physics/
willemsen_bilbao_serafin.rs
