sva_samples/physics/
chaigne_doutaut.rs

1// Concern: one hammer/free-free-bar call site's finite-difference state | Non-concern: argument evaluation, builtin dispatch | IO: (params, sr) -> a site; () -> a sample
2
3//! Chaigne & Doutaut, JASA 101(1) 539-557 (1997): `u_tt = -kappa^2 u_xxxx`, free-free BCs,
4//! uniform cross-section, driven by `hammer.rs`.
5
6use crate::error::SampleError;
7use crate::physics::Solver;
8
9use crate::physics::bound::Bound::*;
10use crate::physics::bound::all;
11use crate::physics::hammer::Hammer;
12
13#[derive(Clone, Debug, PartialEq)]
14pub struct ChaigneDoutautParams {
15    pub f0: f64,
16    pub strike_pos: f64,
17    pub vel: f64,
18    pub hammer_mass: f64,
19    pub hammer_k: f64,
20    pub hammer_p: f64,
21    pub damp_dc: f64,
22    pub damp_freq: f64,
23}
24
25impl ChaigneDoutautParams {
26    /// `chaigne_askenfelt`'s hammer defaults, at the fundamental asked for.
27    pub fn at(f0: f64) -> ChaigneDoutautParams {
28        ChaigneDoutautParams {
29            f0,
30            strike_pos: 0.15,
31            vel: 3.2,
32            hammer_mass: 2.9e-3,
33            hammer_k: 2.6646e8,
34            hammer_p: 2.5,
35            damp_dc: 0.6,
36            damp_freq: 1.6e-4,
37        }
38    }
39}
40
41impl ChaigneDoutautParams {
42    pub fn valid(&self) -> bool {
43        all(&[
44            (self.f0, Positive),
45            (self.strike_pos, OpenUnit),
46            (self.vel, Positive),
47            (self.hammer_mass, Positive),
48            (self.hammer_k, Positive),
49            (self.hammer_p, Positive),
50            (self.damp_dc, NonNegative),
51            (self.damp_freq, NonNegative),
52        ])
53    }
54}
55
56/// Aluminum, tuned to land a C4 bar near ~45cm.
57const BAR_YOUNGS_MODULUS_PA: f64 = 7.0e10;
58const BAR_DENSITY_KG_M3: f64 = 2700.0;
59const BAR_THICKNESS_M: f64 = 0.01;
60/// Cancels out of `kappa`; sets the contact-injection mass scale.
61const BAR_WIDTH_M: f64 = 1.0;
62/// Free-free root of `cos(beta L)cosh(beta L) = 1`.
63const BETA1_L: f64 = 4.730040744862704;
64/// A free-bar anti-node for nearly every low partial, off the strike.
65const PICKUP_POS: f64 = 0.93;
66
67#[derive(Clone)]
68pub(crate) struct BarGrid {
69    pub(crate) u_now: Vec<f64>,
70    pub(crate) u_prev: Vec<f64>,
71    pub(crate) u_next: Vec<f64>,
72    pub(crate) n: usize,
73    pub(crate) dx: f64,
74    pub(crate) rho: f64,
75    #[allow(dead_code)]
76    pub(crate) kappa: f64,
77    pub(crate) stiff_sq: f64,
78    pub(crate) damp_a: f64,
79    pub(crate) damp_b: f64,
80    pub(crate) bar_substeps: usize,
81}
82
83/// Free-end ghosts: `u_xx=0` gives `u_{-1}=2u_0-u_1`, then `u_xxx=0` gives
84/// `u_{-2}=4u_0-4u_1+u_2`.
85pub(crate) fn free_ghost(u: &[f64], n: usize, idx: isize) -> f64 {
86    if idx >= 0 && idx as usize <= n {
87        return u[idx as usize];
88    }
89    if idx < 0 {
90        match idx {
91            -1 => 2.0 * u[0] - u[1],
92            -2 => 4.0 * u[0] - 4.0 * u[1] + u[2],
93            _ => unreachable!(),
94        }
95    } else {
96        match idx as usize - n {
97            1 => 2.0 * u[n] - u[n - 1],
98            2 => 4.0 * u[n] - 4.0 * u[n - 1] + u[n - 2],
99            _ => unreachable!(),
100        }
101    }
102}
103
104/// von Neumann: `mu = kappa*dt/dx^2 <= 0.5`, at a 90% margin.
105pub(crate) fn bar_grid(
106    rho: f64,
107    kappa: f64,
108    length: f64,
109    damp_dc: f64,
110    damp_freq: f64,
111    sr: f64,
112) -> BarGrid {
113    let dt = 1.0 / sr;
114    let mu_max = 0.5;
115    let mu_target = 0.9 * mu_max;
116    let dx_target = (kappa * dt / mu_target).sqrt();
117    let n = ((length / dx_target).round() as usize).max(4);
118    let dx = length / n as f64;
119    // The n>=4 floor can outrun `dx_target` for a short (high-`f0`) bar; sub-step then.
120    let mu_at_dt = kappa * dt / (dx * dx);
121    let bar_substeps = (mu_at_dt / mu_target).ceil().max(1.0) as usize;
122    let dt_sub = dt / bar_substeps as f64;
123    let stiff_sq = kappa * kappa * dt_sub * dt_sub / dx.powi(4);
124
125    BarGrid {
126        u_now: vec![0.0; n + 1],
127        u_prev: vec![0.0; n + 1],
128        u_next: vec![0.0; n + 1],
129        n,
130        dx,
131        rho,
132        kappa,
133        stiff_sq,
134        damp_a: 2.0 * damp_dc * dt_sub,
135        damp_b: 2.0 * damp_freq * dt_sub / (dx * dx),
136        bar_substeps,
137    }
138}
139
140/// `f0 = kappa*(beta1 L)^2/(2 pi L^2)`, solved for `L`.
141fn build_grid(params: &ChaigneDoutautParams, sr: f64) -> BarGrid {
142    let kappa =
143        (BAR_YOUNGS_MODULUS_PA / BAR_DENSITY_KG_M3).sqrt() * BAR_THICKNESS_M / 12.0f64.sqrt();
144    let length = (kappa * BETA1_L * BETA1_L / (2.0 * std::f64::consts::PI * params.f0)).sqrt();
145    let rho = BAR_DENSITY_KG_M3 * BAR_WIDTH_M * BAR_THICKNESS_M;
146    bar_grid(rho, kappa, length, params.damp_dc, params.damp_freq, sr)
147}
148
149#[derive(Clone)]
150pub struct ChaigneDoutautSite {
151    bar: BarGrid,
152    hammer: Hammer,
153    detached: bool,
154    contact_index: usize,
155    pickup_index: usize,
156    dt: f64,
157}
158
159impl ChaigneDoutautSite {
160    pub fn new(params: &ChaigneDoutautParams, sr: f64) -> ChaigneDoutautSite {
161        let bar = build_grid(params, sr);
162        let contact_index = (params.strike_pos * bar.n as f64)
163            .round()
164            .clamp(1.0, (bar.n - 1) as f64) as usize;
165        let pickup_index = (PICKUP_POS * bar.n as f64)
166            .round()
167            .clamp(1.0, (bar.n - 1) as f64) as usize;
168        let dt = 1.0 / sr;
169        ChaigneDoutautSite {
170            bar,
171            hammer: Hammer::new(
172                params.hammer_mass,
173                params.hammer_k,
174                params.hammer_p,
175                params.vel,
176                dt,
177            ),
178            detached: false,
179            contact_index,
180            pickup_index,
181            dt,
182        }
183    }
184}
185
186impl ChaigneDoutautSite {
187    fn advance(&mut self) -> f64 {
188        let bar = &mut self.bar;
189        let n = bar.n;
190        let u_h = bar.u_now[self.contact_index];
191
192        let mut forces = [0.0f64; 1];
193        self.hammer.substeps(
194            self.dt,
195            &[u_h],
196            std::slice::from_mut(&mut self.detached),
197            &mut forces,
198        );
199        let force = forces[0];
200
201        let sample = bar.u_now[self.pickup_index];
202
203        let dt_sub = self.dt / bar.bar_substeps as f64;
204        let injection = (dt_sub * dt_sub / (bar.rho * bar.dx)) * force;
205        for _ in 0..bar.bar_substeps {
206            for i in 0..=n {
207                let ii = i as isize;
208                let lap_now = free_ghost(&bar.u_now, n, ii + 1) - 2.0 * bar.u_now[i]
209                    + free_ghost(&bar.u_now, n, ii - 1);
210                let biharm = free_ghost(&bar.u_now, n, ii + 2)
211                    - 4.0 * free_ghost(&bar.u_now, n, ii + 1)
212                    + 6.0 * bar.u_now[i]
213                    - 4.0 * free_ghost(&bar.u_now, n, ii - 1)
214                    + free_ghost(&bar.u_now, n, ii - 2);
215                let lap_prev = free_ghost(&bar.u_prev, n, ii + 1) - 2.0 * bar.u_prev[i]
216                    + free_ghost(&bar.u_prev, n, ii - 1);
217                let mut next = 2.0 * bar.u_now[i]
218                    - bar.u_prev[i]
219                    - bar.stiff_sq * biharm
220                    - bar.damp_a * (bar.u_now[i] - bar.u_prev[i])
221                    + bar.damp_b * (lap_now - lap_prev);
222                if i == self.contact_index {
223                    next += injection;
224                }
225                bar.u_next[i] = next;
226            }
227            std::mem::swap(&mut bar.u_prev, &mut bar.u_now);
228            std::mem::swap(&mut bar.u_now, &mut bar.u_next);
229        }
230
231        sample
232    }
233}
234
235impl Solver for ChaigneDoutautSite {
236    fn bytes(&self) -> usize {
237        let bar = &self.bar;
238        size_of::<Self>()
239            + super::floats(&bar.u_now)
240            + super::floats(&bar.u_prev)
241            + super::floats(&bar.u_next)
242    }
243
244    fn step(&mut self, _args: &[f64]) -> Result<f64, SampleError> {
245        Ok(self.advance())
246    }
247}
sva_samples/physics/chaigne_doutaut.rs

sva_samples/physics/
chaigne_doutaut.rs
