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 `stiff_string`'s FD string.
6
7use crate::error::SampleError;
8use crate::physics::{Solver, Varies};
9
10/// The bow drives the string: its speed and its force may move every sample.
11pub const VARYING: &[(&str, Varies)] = &[
12    ("bow_vel", Varies::PerSample),
13    ("bow_force", Varies::PerSample),
14];
15
16use crate::physics::bound::Bound::*;
17use crate::physics::bound::all;
18use crate::physics::stiff_string::{
19    StringGrid, Wire, dispersive_grid, grid_tension, point_weights, read_at, spread, stencil_update,
20};
21
22#[derive(Clone, Debug, PartialEq)]
23pub struct WillemsenBilbaoSerafinParams {
24    pub f0: f64,
25    pub b: f64,
26    pub bow_pos: f64,
27    pub bow_vel: f64,
28    pub bow_force: f64,
29    pub mu_s: f64,
30    pub mu_c: f64,
31    pub stribeck_vel: f64,
32    pub bristle_stiffness: f64,
33    pub bristle_damping: f64,
34    pub viscous_friction: f64,
35    pub damp_dc: f64,
36    pub damp_freq: f64,
37}
38
39impl WillemsenBilbaoSerafinParams {
40    /// DAFx-19 Table 1's reference values, at the fundamental asked for.
41    pub fn at(f0: f64) -> WillemsenBilbaoSerafinParams {
42        WillemsenBilbaoSerafinParams {
43            f0,
44            b: 1e-4,
45            bow_pos: 0.25,
46            bow_vel: 0.1,
47            bow_force: 10.0,
48            mu_s: 0.8,
49            mu_c: 0.3,
50            stribeck_vel: 0.1,
51            bristle_stiffness: 1e4,
52            bristle_damping: 0.1,
53            viscous_friction: 0.4,
54            damp_dc: 1.0,
55            damp_freq: 5e-3,
56        }
57    }
58}
59
60impl WillemsenBilbaoSerafinParams {
61    pub fn valid(&self) -> bool {
62        all(&[
63            (self.f0, Positive),
64            (self.b, NonNegative),
65            (self.bow_pos, OpenUnit),
66            (self.bow_vel, Finite),
67            (self.bow_force, Positive),
68            (self.mu_s, Positive),
69            (self.mu_c, Positive),
70            (self.stribeck_vel, Positive),
71            (self.bristle_stiffness, Positive),
72            (self.bristle_damping, NonNegative),
73            (self.viscous_friction, NonNegative),
74            (self.damp_dc, NonNegative),
75            (self.damp_freq, NonNegative),
76        ])
77    }
78}
79
80const WIRE_DENSITY_KG_M3: f64 = 7850.0;
81const WIRE_RADIUS_M: f64 = 5e-4;
82/// Fixed per Table 1; `c = 2 f0 L` scales instead (`chaigne_askenfelt` inverts that).
83const WIRE_LENGTH_M: f64 = 1.0;
84const NR_MAX_ITERATIONS: usize = 50;
85const NR_TOLERANCE: f64 = 1e-7;
86
87/// Table 1's wire at `c = 2 f0 L`.
88fn build_grid(params: &WillemsenBilbaoSerafinParams, sr: f64) -> Option<StringGrid> {
89    let rho = std::f64::consts::PI * WIRE_RADIUS_M * WIRE_RADIUS_M * WIRE_DENSITY_KG_M3;
90    let wire = Wire {
91        rho,
92        c: 2.0 * params.f0 * WIRE_LENGTH_M,
93        length: WIRE_LENGTH_M,
94    };
95    dispersive_grid(
96        wire,
97        params.f0,
98        params.b,
99        params.damp_dc,
100        params.damp_freq,
101        sr,
102    )
103}
104
105fn sgn(x: f64) -> f64 {
106    if x > 0.0 {
107        1.0
108    } else if x < 0.0 {
109        -1.0
110    } else {
111        0.0
112    }
113}
114
115/// Eq. (7): `abs()` around `z_ss`, missing in the pre-DAFx-19 literature.
116fn steady_state(v: f64, s0: f64, f_c: f64, f_s: f64, stribeck_vel: f64) -> f64 {
117    sgn(v) / s0 * (f_c + (f_s - f_c) * (-(v / stribeck_vel).powi(2)).exp())
118}
119
120/// Eq. (8)-(9): `sgn(z)` on the sine offset; `|z|=z_ba`/`|z_ss|` resolved, not undefined.
121fn adhesion(v: f64, z: f64, z_ba: f64, zss: f64) -> f64 {
122    if sgn(v) != sgn(z) {
123        return 0.0;
124    }
125    let (az, azss) = (z.abs(), zss.abs());
126    if az <= z_ba {
127        0.0
128    } else if az >= azss {
129        1.0
130    } else {
131        let span = (azss - z_ba).max(1e-300);
132        let s = sgn(z);
133        0.5 * (1.0 + s * (std::f64::consts::PI * (z - s * 0.5 * (azss + z_ba)) / span).sin())
134    }
135}
136
137/// Eq. (6)/(15). `v == 0` avoids dividing by `z_ss(0) = 0`.
138fn bristle_rate(v: f64, z: f64, s0: f64, f_c: f64, f_s: f64, stribeck_vel: f64, z_ba: f64) -> f64 {
139    if v == 0.0 {
140        return 0.0;
141    }
142    let zss = steady_state(v, s0, f_c, f_s, stribeck_vel);
143    let a = adhesion(v, z, z_ba, zss);
144    v * (1.0 - a * z / zss)
145}
146
147/// Eq. (4), noise term `s3*w` dropped (disclosed v1 gap).
148fn friction_force(v: f64, z: f64, r: f64, s0: f64, s1: f64, s2: f64) -> f64 {
149    s0 * z + s1 * r + s2 * v
150}
151
152#[derive(Clone)]
153struct Coupling {
154    coeff: f64,
155    b_known: f64,
156    s0: f64,
157    s1: f64,
158    s2: f64,
159    f_c: f64,
160    f_s: f64,
161    stribeck_vel: f64,
162    z_ba: f64,
163    z_prev: f64,
164    r_prev: f64,
165    dt: f64,
166}
167
168impl Coupling {
169    /// `g1` re-derives eq. (17)-(19): theirs assume centred damping, the FD string's is
170    /// backward. `g2` keeps their eq. (20)-(21).
171    fn residual(&self, v: f64, z: f64) -> (f64, f64) {
172        let r = bristle_rate(
173            v,
174            z,
175            self.s0,
176            self.f_c,
177            self.f_s,
178            self.stribeck_vel,
179            self.z_ba,
180        );
181        let f = friction_force(v, z, r, self.s0, self.s1, self.s2);
182        let g1 = v + self.coeff * f - self.b_known;
183        let a = 2.0 * (z - self.z_prev) / self.dt - self.r_prev;
184        (g1, r - a)
185    }
186
187    /// A central-difference Jacobian, not Algorithm 1's analytic one.
188    fn newton_step(&self, v: f64, z: f64) -> Option<(f64, f64)> {
189        let (g1, g2) = self.residual(v, z);
190        let hv = v.abs().max(1e-3) * 1e-6;
191        let hz = z.abs().max(self.z_ba).max(1e-9) * 1e-6;
192        let (g1_vp, g2_vp) = self.residual(v + hv, z);
193        let (g1_vm, g2_vm) = self.residual(v - hv, z);
194        let (g1_zp, g2_zp) = self.residual(v, z + hz);
195        let (g1_zm, g2_zm) = self.residual(v, z - hz);
196        let dg1_dv = (g1_vp - g1_vm) / (2.0 * hv);
197        let dg2_dv = (g2_vp - g2_vm) / (2.0 * hv);
198        let dg1_dz = (g1_zp - g1_zm) / (2.0 * hz);
199        let dg2_dz = (g2_zp - g2_zm) / (2.0 * hz);
200        let det = dg1_dv * dg2_dz - dg1_dz * dg2_dv;
201        if !det.is_finite() || det.abs() < 1e-300 {
202            return None;
203        }
204        let dv = (g1 * dg2_dz - g2 * dg1_dz) / det;
205        let dz = (dg1_dv * g2 - dg2_dv * g1) / det;
206        let (new_v, new_z) = (v - dv, z - dz);
207        if !new_v.is_finite() || !new_z.is_finite() {
208            return None;
209        }
210        Some((new_v, new_z))
211    }
212}
213
214#[derive(Clone)]
215pub struct WillemsenBilbaoSerafinSite {
216    strings: Vec<StringGrid>,
217    /// Read and spread alike.
218    bow: Vec<f64>,
219    bow_self: f64,
220    tension: f64,
221    dt: f64,
222    bow_vel: f64,
223    /// The `(bow_vel, bow_force)` the site was opened with, and the force `f_c`, `f_s` and
224    /// `z_ba` are derived for.
225    own: (f64, f64),
226    force: f64,
227    mu: (f64, f64),
228    f_c: f64,
229    f_s: f64,
230    stribeck_vel: f64,
231    s0: f64,
232    s1: f64,
233    s2: f64,
234    z_ba: f64,
235    z: f64,
236    z_prev: f64,
237    /// `r^{n-1}`, eq. (20)-(21)'s extra trapezoidal state.
238    r_prev: f64,
239    last_v: f64,
240    last_f: f64,
241    /// Samples stepped, which a refusal names.
242    steps: usize,
243}
244
245impl WillemsenBilbaoSerafinSite {
246    pub fn new(params: &WillemsenBilbaoSerafinParams, sr: f64) -> Result<Self, SampleError> {
247        let grid = build_grid(params, sr).ok_or(SampleError::StringPastRate {
248            model: "willemsen_bilbao_serafin",
249        })?;
250        let bow = point_weights(grid.n, params.bow_pos);
251        let bow_self = read_at(&bow, &bow);
252        let tension = grid_tension(&grid, 1.0 / sr);
253        let f_c = params.mu_c * params.bow_force;
254        let f_s = params.mu_s * params.bow_force;
255        // Table 1: z_ba = 0.7*f_C/s0 — off the kinetic (mu_c) force, not the static one.
256        let z_ba = 0.7 * f_c / params.bristle_stiffness;
257        Ok(WillemsenBilbaoSerafinSite {
258            strings: vec![grid],
259            bow,
260            bow_self,
261            tension,
262            dt: 1.0 / sr,
263            bow_vel: params.bow_vel,
264            own: (params.bow_vel, params.bow_force),
265            force: params.bow_force,
266            mu: (params.mu_c, params.mu_s),
267            f_c,
268            f_s,
269            stribeck_vel: params.stribeck_vel,
270            s0: params.bristle_stiffness,
271            s1: params.bristle_damping,
272            s2: params.viscous_friction,
273            z_ba,
274            z: 0.0,
275            z_prev: 0.0,
276            r_prev: 0.0,
277            last_v: 0.0,
278            last_f: 0.0,
279            steps: 0,
280        })
281    }
282}
283
284impl Solver for WillemsenBilbaoSerafinSite {
285    fn bytes(&self) -> usize {
286        let grids: usize = self.strings.iter().map(StringGrid::bytes).sum();
287        size_of::<Self>() + grids + super::floats(&self.bow)
288    }
289
290    fn step(&mut self, args: &[f64]) -> Result<f64, SampleError> {
291        let vel = args.first().copied().unwrap_or(self.own.0);
292        let force = args.get(1).copied().unwrap_or(self.own.1);
293        for ((name, _), value, holds) in [
294            (VARYING[0], vel, vel.is_finite()),
295            (VARYING[1], force, force > 0.0 && force.is_finite()),
296        ] {
297            if !holds {
298                return Err(SampleError::ArgumentOutOfRange {
299                    model: "willemsen_bilbao_serafin",
300                    name,
301                    bits: value.to_bits(),
302                    sample: self.steps as u64,
303                });
304            }
305        }
306        self.bow_vel = vel;
307        if force.to_bits() != self.force.to_bits() {
308            self.force = force;
309            self.f_c = self.mu.0 * force;
310            self.f_s = self.mu.1 * force;
311            self.z_ba = 0.7 * self.f_c / self.s0;
312        }
313        let dt = self.dt;
314        let coupling_prev = (self.z, self.r_prev);
315
316        let grid = &mut self.strings[0];
317        let n = grid.n;
318        for i in 1..n {
319            grid.y_next[i] = stencil_update(grid, i, 0.0, 0.0);
320        }
321        grid.y_next[0] = 0.0;
322        grid.y_next[n] = 0.0;
323        let free_next = read_at(&self.bow, &grid.y_next);
324        let bow_prev = read_at(&self.bow, &grid.y_prev);
325
326        let coupling = Coupling {
327            coeff: dt * self.bow_self / (2.0 * grid.rho * grid.dx),
328            b_known: (free_next - bow_prev) / (2.0 * dt) - self.bow_vel,
329            s0: self.s0,
330            s1: self.s1,
331            s2: self.s2,
332            f_c: self.f_c,
333            f_s: self.f_s,
334            stribeck_vel: self.stribeck_vel,
335            z_ba: self.z_ba,
336            z_prev: coupling_prev.0,
337            r_prev: coupling_prev.1,
338            dt,
339        };
340
341        let (mut v, mut z) = (self.last_v, self.z);
342        let mut converged = false;
343        for _ in 0..NR_MAX_ITERATIONS {
344            let Some((new_v, new_z)) = coupling.newton_step(v, z) else {
345                break;
346            };
347            let step_norm = ((new_v - v).powi(2) + (new_z - z).powi(2)).sqrt();
348            v = new_v;
349            z = new_z;
350            if step_norm < NR_TOLERANCE {
351                converged = true;
352                break;
353            }
354        }
355        if !converged {
356            return Err(SampleError::ContactUnsettled {
357                model: "willemsen_bilbao_serafin",
358                sample: self.steps,
359            });
360        }
361        self.steps += 1;
362
363        let r = bristle_rate(
364            v,
365            z,
366            self.s0,
367            self.f_c,
368            self.f_s,
369            self.stribeck_vel,
370            self.z_ba,
371        );
372        let f = friction_force(v, z, r, self.s0, self.s1, self.s2);
373        let grid = &mut self.strings[0];
374        spread(grid, &self.bow, -f, dt);
375
376        self.z_prev = self.z;
377        self.z = z;
378        self.r_prev = r;
379        self.last_v = v;
380        self.last_f = f;
381
382        let sample = self.tension * (grid.y_now[n] - grid.y_now[n - 1]) / grid.dx;
383
384        std::mem::swap(&mut grid.y_prev, &mut grid.y_now);
385        std::mem::swap(&mut grid.y_now, &mut grid.y_next);
386
387        Ok(sample)
388    }
389}
sva_samples/physics/willemsen_bilbao_serafin.rs

sva_samples/physics/
willemsen_bilbao_serafin.rs
