cellular_raza_core/backend/chili/
solvers.rs

1use cellular_raza_concepts::{CalcError, Xapy};
2use num::FromPrimitive;
3
4#[cfg(feature = "tracing")]
5use tracing::instrument;
6
7use super::{UpdateReactions, UpdateReactionsContact};
8
9/// Classical euler solver for the [Mechanics](cellular_raza_concepts::Mechanics) trait.
10///
11/// The euler solver is the most simple solver and not stable for many problems.
12/// Thus its usage is discouraged. For another general-purpose solver look at
13/// [mechanics_adams_bashforth_2].
14///
15/// The update step follows the simple equations
16/// \\begin{align}
17///     x(t_{i+1}) &= x(t_i) + \Delta t \frac{d x}{d t}(t_i)\\\\
18///     v(t_{i+1}) &= v(t_i) + \Delta t \frac{d v}{d t}(t_i)
19/// \\end{align}
20/// where $\Delta t$ is the step size and $dx/dt$ and $dv/dt$ are calculated by the
21/// [calculate_increment](cellular_raza_concepts::Mechanics::calculate_increment) method.
22#[cfg_attr(feature = "tracing", instrument(skip_all))]
23pub fn mechanics_euler<C, A, Pos, Vel, For, Float>(
24    cell: &mut C,
25    aux_storage: &mut A,
26    dt: Float,
27    rng: &mut rand_chacha::ChaCha8Rng,
28) -> Result<(), super::SimulationError>
29where
30    A: super::aux_storage::UpdateMechanics<Pos, Vel, For, 0>,
31    C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
32    C: cellular_raza_concepts::Position<Pos>,
33    C: cellular_raza_concepts::Velocity<Vel>,
34    Pos: Xapy<Float> + Clone,
35    Vel: Xapy<Float> + Clone,
36    Float: num::Float + Copy,
37{
38    let force = aux_storage.get_current_force_and_reset();
39    let velocity = cell.velocity();
40    let position = cell.pos();
41
42    let (dx, dv) = cell.calculate_increment(force)?;
43    let (dx_rand, dv_rand) = cell.get_random_contribution(rng, dt)?;
44
45    // Update values in the aux_storage
46    aux_storage.set_last_position(dx.clone());
47    aux_storage.set_last_velocity(dv.clone());
48
49    // Calculate new position and velocity of cell
50    let new_position = euler(position, dx, dt, dx_rand)?;
51    let new_velocity = euler(velocity, dv, dt, dv_rand)?;
52    cell.set_pos(&new_position);
53    cell.set_velocity(&new_velocity);
54    Ok(())
55}
56
57/// Note that the const generic for this struct is the order of the solver minus one.
58/// This is due to the fact that the AuxStorage only stores one less step than the order of the
59/// solver.
60pub(crate) struct MechanicsAdamsBashforthSolver<const N: usize>;
61
62pub(crate) trait AdamsBashforth<const N: usize> {
63    #[allow(unused)]
64    fn update<C, A, Pos, Vel, For, Float>(
65        cell: &mut C,
66        aux_storage: &mut A,
67        dt: Float,
68        rng: &mut rand_chacha::ChaCha8Rng,
69    ) -> Result<(), super::SimulationError>
70    where
71        A: super::aux_storage::UpdateMechanics<Pos, Vel, For, N>,
72        C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
73        C: cellular_raza_concepts::Position<Pos>,
74        C: cellular_raza_concepts::Velocity<Vel>,
75        Pos: Xapy<Float> + Clone,
76        Vel: Xapy<Float> + Clone,
77        Float: num::Float + FromPrimitive;
78}
79
80impl AdamsBashforth<2> for MechanicsAdamsBashforthSolver<2> {
81    #[allow(unused)]
82    fn update<C, A, Pos, Vel, For, Float>(
83        cell: &mut C,
84        aux_storage: &mut A,
85        dt: Float,
86        rng: &mut rand_chacha::ChaCha8Rng,
87    ) -> Result<(), super::SimulationError>
88    where
89        A: super::aux_storage::UpdateMechanics<Pos, Vel, For, 2>,
90        C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
91        C: cellular_raza_concepts::Position<Pos>,
92        C: cellular_raza_concepts::Velocity<Vel>,
93        Pos: Xapy<Float> + Clone,
94        Vel: Xapy<Float> + Clone,
95        Float: num::Float + FromPrimitive,
96    {
97        mechanics_adams_bashforth_3(cell, aux_storage, dt, rng)
98    }
99}
100
101impl AdamsBashforth<1> for MechanicsAdamsBashforthSolver<1> {
102    fn update<C, A, Pos, Vel, For, Float>(
103        cell: &mut C,
104        aux_storage: &mut A,
105        dt: Float,
106        rng: &mut rand_chacha::ChaCha8Rng,
107    ) -> Result<(), super::SimulationError>
108    where
109        A: super::aux_storage::UpdateMechanics<Pos, Vel, For, 1>,
110        C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
111        C: cellular_raza_concepts::Position<Pos>,
112        C: cellular_raza_concepts::Velocity<Vel>,
113        Pos: Xapy<Float> + Clone,
114        Vel: Xapy<Float> + Clone,
115        Float: num::Float + FromPrimitive,
116    {
117        mechanics_adams_bashforth_2(cell, aux_storage, dt, rng)
118    }
119}
120
121impl AdamsBashforth<0> for MechanicsAdamsBashforthSolver<0> {
122    fn update<C, A, Pos, Vel, For, Float>(
123        cell: &mut C,
124        aux_storage: &mut A,
125        dt: Float,
126        rng: &mut rand_chacha::ChaCha8Rng,
127    ) -> Result<(), super::SimulationError>
128    where
129        A: super::aux_storage::UpdateMechanics<Pos, Vel, For, 0>,
130        C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
131        C: cellular_raza_concepts::Position<Pos>,
132        C: cellular_raza_concepts::Velocity<Vel>,
133        Pos: Xapy<Float> + Clone,
134        Vel: Xapy<Float> + Clone,
135        Float: num::Float + FromPrimitive,
136    {
137        mechanics_euler(cell, aux_storage, dt, rng)
138    }
139}
140
141/// Three-step Adams-Bashforth method.
142///
143/// See also the [Wikipedia](https://en.wikipedia.org/wiki/Linear_multistep_method) article.
144/// We track previous increments of the update steps and use these in order to update the next time
145/// steps.
146///
147/// The equations for updating are given by
148/// \\begin{equation}
149///     y(t_{i+3}) = y(t_{i+2}) + \Delta t\left(\frac{23}{12}\frac{dy}{dt}(t_{i+2})  - \frac{16}{12}\frac{dy}{dt}(t_{i+1}) + \frac{5}{12}\frac{dy}{dt}(t_i)\right)
150/// \\end{equation}
151///
152/// for both the position and velocity.
153/// In the beginning of the simulation, when not enough previous increment values are known,
154/// we resort to the [mechanics_adams_bashforth_2] and [mechanics_euler] solver.
155#[cfg_attr(feature = "tracing", instrument(skip_all))]
156pub fn mechanics_adams_bashforth_3<C, A, Pos, Vel, For, Float>(
157    cell: &mut C,
158    aux_storage: &mut A,
159    dt: Float,
160    rng: &mut rand_chacha::ChaCha8Rng,
161) -> Result<(), super::SimulationError>
162where
163    A: super::aux_storage::UpdateMechanics<Pos, Vel, For, 2>,
164    C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
165    C: cellular_raza_concepts::Position<Pos>,
166    C: cellular_raza_concepts::Velocity<Vel>,
167    Pos: Xapy<Float> + Clone,
168    Vel: Xapy<Float> + Clone,
169    Float: num::Float + FromPrimitive,
170{
171    let force = aux_storage.get_current_force_and_reset();
172    let velocity = cell.velocity();
173    let position = cell.pos();
174
175    let (dx, dv) = cell.calculate_increment(force)?;
176    let (dx_rand, dv_rand) = cell.get_random_contribution(rng, dt)?;
177
178    // Update values in the aux_storage
179    aux_storage.set_last_position(dx.clone());
180    aux_storage.set_last_velocity(dv.clone());
181
182    // Calculate new position and velocity of cell
183    let n_previous_values = aux_storage.n_previous_values();
184    let mut old_pos_increments = aux_storage.previous_positions();
185    let mut old_vel_increments = aux_storage.previous_velocities();
186    let (new_position, new_velocity) = match n_previous_values {
187        2 => (
188            adams_bashforth_3(
189                position,
190                [
191                    dx,
192                    old_pos_increments.next().unwrap().clone(),
193                    old_pos_increments.next().unwrap().clone(),
194                ],
195                dt,
196                dx_rand,
197            )?,
198            adams_bashforth_3(
199                velocity,
200                [
201                    dv,
202                    old_vel_increments.next().unwrap().clone(),
203                    old_vel_increments.next().unwrap().clone(),
204                ],
205                dt,
206                dv_rand,
207            )?,
208        ),
209        1 => (
210            adams_bashforth_2(
211                position,
212                [dx, old_pos_increments.next().unwrap().clone()],
213                dt,
214                dx_rand,
215            )?,
216            adams_bashforth_2(
217                velocity,
218                [dv, old_vel_increments.next().unwrap().clone()],
219                dt,
220                dv_rand,
221            )?,
222        ),
223        _ => (
224            euler(position, dx, dt, dx_rand)?,
225            euler(velocity, dv, dt, dv_rand)?,
226        ),
227    };
228    cell.set_pos(&new_position);
229    cell.set_velocity(&new_velocity);
230    Ok(())
231}
232
233/// Two-step Adams-Bashforth method.
234///
235/// See also the [Wikipedia](https://en.wikipedia.org/wiki/Linear_multistep_method) article.
236/// We track previous increments of the update steps and use these in order to update the next time
237/// steps.
238///
239/// The equations for updating are given by
240/// \\begin{equation}
241///     y(t_{i+2}) = y(t_{i+1}) + \Delta t\left(\frac{3}{12}\frac{dy}{dt}(t_{i+1})  - \frac{1}{2}\frac{dy}{dt}(t_{i})\right)
242/// \\end{equation}
243///
244/// for both the position and velocity.
245/// In the beginning of the simulation, when not enough previous increment values are known,
246/// we resort to the [euler](mechanics_euler) solver.
247pub fn mechanics_adams_bashforth_2<C, A, Pos, Vel, For, Float>(
248    cell: &mut C,
249    aux_storage: &mut A,
250    dt: Float,
251    rng: &mut rand_chacha::ChaCha8Rng,
252) -> Result<(), super::SimulationError>
253where
254    A: super::aux_storage::UpdateMechanics<Pos, Vel, For, 1>,
255    C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
256    C: cellular_raza_concepts::Position<Pos>,
257    C: cellular_raza_concepts::Velocity<Vel>,
258    Pos: Xapy<Float> + Clone,
259    Vel: Xapy<Float> + Clone,
260    Float: num::Float + FromPrimitive,
261{
262    let force = aux_storage.get_current_force_and_reset();
263    let velocity = cell.velocity();
264    let position = cell.pos();
265
266    let (dx, dv) = cell.calculate_increment(force)?;
267    let (dx_rand, dv_rand) = cell.get_random_contribution(rng, dt)?;
268
269    // Update values in the aux_storage
270    aux_storage.set_last_position(dx.clone());
271    aux_storage.set_last_velocity(dv.clone());
272
273    // Calculate new position and velocity of cell
274    let n_previous_values = aux_storage.n_previous_values();
275    let mut old_pos_increments = aux_storage.previous_positions();
276    let mut old_vel_increments = aux_storage.previous_velocities();
277    let (new_position, new_velocity) = match n_previous_values {
278        1 => (
279            adams_bashforth_2(
280                position,
281                [dx, old_pos_increments.next().unwrap().clone()],
282                dt,
283                dx_rand,
284            )?,
285            adams_bashforth_2(
286                velocity,
287                [dv, old_vel_increments.next().unwrap().clone()],
288                dt,
289                dv_rand,
290            )?,
291        ),
292        _ => (
293            euler(position, dx, dt, dx_rand)?,
294            euler(velocity, dv, dt, dv_rand)?,
295        ),
296    };
297    cell.set_pos(&new_position);
298    cell.set_velocity(&new_velocity);
299    Ok(())
300}
301
302#[inline]
303fn euler<X, F>(x: X, dx: X, dt: F, dx_rand: X) -> Result<X, CalcError>
304where
305    X: Xapy<F>,
306    F: num::Float + Copy,
307{
308    let x_new = dx.xapy(dt, &x).xapy(F::one(), &dx_rand.xa(dt));
309    Ok(x_new)
310}
311
312#[inline]
313fn adams_bashforth_3<X, F>(x: X, dx: [X; 3], dt: F, dx_rand: X) -> Result<X, CalcError>
314where
315    X: Xapy<F>,
316    F: Copy + FromPrimitive + num::Float,
317{
318    let f0 = F::from_isize(23).unwrap() / F::from_isize(12).unwrap();
319    let f1 = -F::from_isize(16).unwrap() / F::from_isize(12).unwrap();
320    let f2 = F::from_isize(5).unwrap() / F::from_isize(12).unwrap();
321
322    let [dx0, dx1, dx2] = dx;
323
324    let x_new = dx0
325        .xapy(f0, &dx1.xapy(f1, &dx2.xa(f2)))
326        .xapy(dt, &x)
327        .xapy(F::one(), &dx_rand.xa(dt));
328    Ok(x_new)
329}
330
331#[inline]
332fn adams_bashforth_2<X, F>(x: X, dx: [X; 2], dt: F, dx_rand: X) -> Result<X, CalcError>
333where
334    X: Xapy<F>,
335    F: Copy + FromPrimitive + num::Float,
336{
337    let f0 = F::from_isize(3).unwrap() / F::from_isize(2).unwrap();
338    let f1 = -F::from_isize(1).unwrap() / F::from_isize(2).unwrap();
339
340    let [dx0, dx1] = dx;
341
342    let x_new = dx0
343        .xapy(f0, &dx1.xa(f1))
344        .xapy(dt, &x)
345        .xapy(F::one(), &dx_rand.xa(dt));
346    Ok(x_new)
347}
348
349/// Note that the const generic for this struct is the order of the solver minus one.
350/// This is due to the fact that the AuxStorage only stores one less step than the order of the
351/// solver.
352pub(crate) struct ReactionsRungeKuttaSolver<const N: usize>;
353
354pub(crate) trait RungeKutta<const N: usize> {
355    #[allow(unused)]
356    fn update<C, A, Ri, Float>(
357        cell: &mut C,
358        aux_storage: &mut A,
359        dt: Float,
360    ) -> Result<(), super::SimulationError>
361    where
362        A: UpdateReactions<Ri>,
363        C: cellular_raza_concepts::Reactions<Ri>,
364        Float: num::Float,
365        Ri: Xapy<Float>;
366}
367
368impl RungeKutta<1> for ReactionsRungeKuttaSolver<1> {
369    #[allow(unused)]
370    #[inline]
371    #[cfg_attr(feature = "tracing", instrument(skip_all))]
372    fn update<C, A, Ri, Float>(
373        cell: &mut C,
374        aux_storage: &mut A,
375        dt: Float,
376    ) -> Result<(), super::SimulationError>
377    where
378        A: UpdateReactions<Ri>,
379        C: cellular_raza_concepts::Reactions<Ri>,
380        Float: num::Float,
381        Ri: Xapy<Float>,
382    {
383        // Constants
384        let intra = cell.get_intracellular();
385
386        // Calculate the intermediate steps
387        let dintra = cell.calculate_intracellular_increment(&intra)?;
388
389        // Update the internal value of the cell
390        aux_storage.incr_conc(dintra);
391        Ok(())
392    }
393}
394
395impl RungeKutta<2> for ReactionsRungeKuttaSolver<2> {
396    #[allow(unused)]
397    #[inline]
398    #[cfg_attr(feature = "tracing", instrument(skip_all))]
399    fn update<C, A, Ri, Float>(
400        cell: &mut C,
401        aux_storage: &mut A,
402        dt: Float,
403    ) -> Result<(), super::SimulationError>
404    where
405        A: UpdateReactions<Ri>,
406        C: cellular_raza_concepts::Reactions<Ri>,
407        Float: num::Float,
408        Ri: Xapy<Float>,
409    {
410        // Constants
411        let two = Float::one() + Float::one();
412        let intra = cell.get_intracellular();
413
414        // Calculate the intermediate steps
415        let dintra1 = cell.calculate_intracellular_increment(&intra)?;
416        let dintra = cell.calculate_intracellular_increment(&dintra1.xapy(dt / two, &intra))?;
417
418        // Update the internal value of the cell
419        aux_storage.incr_conc(dintra);
420        Ok(())
421    }
422}
423
424impl RungeKutta<4> for ReactionsRungeKuttaSolver<4> {
425    #[allow(unused)]
426    #[inline]
427    #[cfg_attr(feature = "tracing", instrument(skip_all))]
428    fn update<C, A, Ri, Float>(
429        cell: &mut C,
430        aux_storage: &mut A,
431        dt: Float,
432    ) -> Result<(), super::SimulationError>
433    where
434        A: UpdateReactions<Ri>,
435        C: cellular_raza_concepts::Reactions<Ri>,
436        Float: num::Float,
437        Ri: Xapy<Float>,
438    {
439        // Constants
440        let two = Float::one() + Float::one();
441        let six = two + two + two;
442
443        let intra = cell.get_intracellular();
444
445        // Calculate the intermediate steps
446        let dintra1 = cell.calculate_intracellular_increment(&intra)?;
447        let dintra2 = cell.calculate_intracellular_increment(&dintra1.xapy(dt / two, &intra))?;
448        let dintra3 = cell.calculate_intracellular_increment(&dintra2.xapy(dt / two, &intra))?;
449        let dintra4 = cell.calculate_intracellular_increment(&dintra3.xapy(dt, &intra))?;
450        let dintra = dintra1.xapy(
451            Float::one() / six,
452            &dintra2.xapy(
453                two / six,
454                &dintra3.xapy(two / six, &dintra4.xa(Float::one() / six)),
455            ),
456        );
457
458        // Update the internal value of the cell
459        aux_storage.incr_conc(dintra);
460        Ok(())
461    }
462}
463
464/// Calculates the increment introduced by the
465/// [ReactionsContact](cellular_raza_concepts::ReactionsContact) aspect.
466#[inline]
467#[cfg_attr(feature = "tracing", instrument(skip_all))]
468pub fn reactions_contact_adams_bashforth_3rd<C, A, F, Ri>(
469    _cell: &mut C,
470    aux_storage: &mut A,
471) -> Result<(), CalcError>
472where
473    A: UpdateReactions<Ri>,
474    A: UpdateReactionsContact<Ri, 2>,
475    Ri: Xapy<F> + Clone,
476    F: FromPrimitive + num::Float,
477{
478    // let dintra0 = aux_storage
479    let dintra = aux_storage.get_conc();
480    let conc_zero = dintra.xa(F::zero());
481    let dintra_contact = <A as UpdateReactionsContact<Ri, 2>>::get_current_increment(&aux_storage);
482    aux_storage.set_last_increment(dintra_contact.clone());
483    let n_previous_values = aux_storage.n_previous_values();
484    let mut old_intracellular_increments = aux_storage.previous_increments();
485    let dintra_new = match n_previous_values {
486        2 => adams_bashforth_3(
487            dintra,
488            [
489                dintra_contact,
490                old_intracellular_increments.next().unwrap().clone(),
491                old_intracellular_increments.next().unwrap().clone(),
492            ],
493            F::one(),
494            conc_zero,
495        )?,
496        1 => adams_bashforth_2(
497            dintra,
498            [
499                dintra_contact,
500                old_intracellular_increments.next().unwrap().clone(),
501            ],
502            F::one(),
503            conc_zero,
504        )?,
505        _ => euler(dintra, dintra_contact, F::one(), conc_zero)?,
506    };
507    aux_storage.set_conc(dintra_new);
508    Ok(())
509}
510
511#[cfg(test)]
512mod test_solvers_reactions {
513    use rand::SeedableRng;
514
515    use crate::backend::chili::{UpdateReactions, local_reactions_use_increment};
516
517    use super::*;
518
519    #[test]
520    fn exponential_decay_rk4() -> Result<(), super::super::SimulationError> {
521        use cellular_raza_concepts::*;
522        struct PlainCell {
523            intracellular: f64,
524            lambda: f64,
525        }
526        impl Intracellular<f64> for PlainCell {
527            fn get_intracellular(&self) -> f64 {
528                self.intracellular
529            }
530            fn set_intracellular(&mut self, intracellular: f64) {
531                self.intracellular = intracellular;
532            }
533        }
534        impl Reactions<f64> for PlainCell {
535            fn calculate_intracellular_increment(
536                &self,
537                intracellular: &f64,
538            ) -> Result<f64, CalcError> {
539                Ok(-self.lambda * intracellular)
540            }
541        }
542        struct AuxStorage {
543            increment: f64,
544        }
545        impl UpdateReactions<f64> for AuxStorage {
546            fn get_conc(&self) -> f64 {
547                self.increment
548            }
549            fn incr_conc(&mut self, incr: f64) {
550                self.increment += incr;
551            }
552            fn set_conc(&mut self, conc: f64) {
553                self.increment = conc;
554            }
555        }
556        let y0 = 33.0;
557        let lambda = 0.2;
558        let dt = 0.1;
559        let exact_solution = |t: f64| -> f64 { y0 * (-lambda * t).exp() };
560
561        let mut cell = PlainCell {
562            intracellular: y0,
563            lambda,
564        };
565        let mut aux_storage = AuxStorage { increment: 0.0 };
566        let mut results_cr = vec![(0.0, cell.get_intracellular())];
567        let mut t = 0.0;
568        for _ in 0..100 {
569            ReactionsRungeKuttaSolver::<4>::update(&mut cell, &mut aux_storage, dt)?;
570            // This rng is just a placeholder and will not be used.
571            let mut _rng = rand_chacha::ChaCha8Rng::seed_from_u64(0);
572            local_reactions_use_increment(&mut cell, &mut aux_storage, dt, &mut _rng)?;
573            t += dt;
574            results_cr.push((t, cell.get_intracellular()));
575        }
576        for (t, res) in results_cr {
577            let res_exact = exact_solution(t);
578            assert!((res - res_exact).abs() < 1e-6);
579        }
580        Ok(())
581    }
582}
583
584#[cfg(test)]
585mod test_solvers {
586    use super::*;
587
588    #[test]
589    fn euler_exp_decay() {
590        let y0 = 27.0;
591        let lambda = 0.3;
592        let dt = 0.001;
593        let rhs = |y: f64| -> f64 { -lambda * y };
594        let exact_solution = |t: f64| -> f64 { y0 * (-lambda * t).exp() };
595
596        // The following truncation error is taken from wikipedia:
597        // https://en.wikipedia.org/wiki/Euler_method#Global_truncation_error
598        // Lipschitz constant of RHS
599        let lipschitz_constant = lambda;
600        // Upper bound on second derivative of y
601        let second_derivative_bound = y0 * lambda.powi(2);
602
603        let expected_global_truncation_error = |t: f64| -> f64 {
604            dt * second_derivative_bound / (2.0 * lipschitz_constant)
605                * ((lipschitz_constant * t).exp() - 1.0)
606        };
607
608        let mut y = y0;
609        let mut t = 0.0;
610        for _ in 0..100 {
611            let dy = rhs(y);
612            y = euler(y, dy, dt, 0.0).unwrap();
613            t += dt;
614            let e = expected_global_truncation_error(t);
615            assert!((y - exact_solution(t)).abs() < e);
616        }
617    }
618
619    #[test]
620    fn adams_bashforth_2_harmonic_oscillator() {
621        #[derive(Clone, Copy, Debug)]
622        struct Vec2(f32, f32);
623        impl<'a> core::ops::Add<&'a Vec2> for Vec2 {
624            type Output = Vec2;
625            fn add(self, rhs: &'a Vec2) -> Self::Output {
626                Vec2(self.0 + rhs.0, self.1 + rhs.1)
627            }
628        }
629        impl core::ops::Add<Vec2> for Vec2 {
630            type Output = Self;
631            fn add(self, rhs: Vec2) -> Self::Output {
632                Vec2(self.0 + rhs.0, self.1 + rhs.1)
633            }
634        }
635        impl<'a> core::ops::Mul<f32> for &'a Vec2 {
636            type Output = Vec2;
637            fn mul(self, rhs: f32) -> Self::Output {
638                Vec2(self.0 * rhs, self.1 * rhs)
639            }
640        }
641        impl core::ops::Neg for Vec2 {
642            type Output = Vec2;
643            fn neg(self) -> Self::Output {
644                Self(-self.0, -self.1)
645            }
646        }
647        impl num::Zero for Vec2 {
648            fn is_zero(&self) -> bool {
649                self.0.is_zero() && self.1.is_zero()
650            }
651            fn set_zero(&mut self) {
652                self.0 = 0.0;
653                self.1 = 0.0;
654            }
655            fn zero() -> Self {
656                Vec2(0., 0.)
657            }
658        }
659
660        // Define parameters and initial values
661        let y0 = Vec2(2.8347, 0.0);
662        let omega: f32 = 0.319;
663        let dt: f32 = 0.045;
664
665        // Write down rhs and exact solution
666        let rhs = |y: Vec2| -> Vec2 { Vec2(y.1, -omega.powi(2) * y.0) };
667        let exact_solution =
668            |t: f32| -> Vec2 { Vec2(y0.0 * (omega * t).cos(), -y0.0 * omega * (omega * t).sin()) };
669        // This is taken from this math.stackexchange post:
670        // https://math.stackexchange.com/questions/1326502/determine-the-local-truncation-error-of-the-following-method
671        // Third order derivatives
672        let third_derivative_bound = Vec2(y0.0 * omega.powi(3), y0.0 * omega.powi(3));
673        let lipschitz_constant = Vec2(1.0, omega);
674        let local_truncation_error = &third_derivative_bound * (5f32 / 12.0 * dt.powi(2));
675        // See this wikipedia article:
676        // https://en.wikipedia.org/wiki/Truncation_error_(numerical_integration)#Relationship_between_local_and_global_truncation_errors
677        let global_truncation_error = |t: f32| -> Vec2 {
678            Vec2(
679                ((lipschitz_constant.0 * t).exp() - 1.0) * local_truncation_error.0
680                    / dt
681                    / lipschitz_constant.0,
682                ((lipschitz_constant.1 * t).exp() - 1.0) * local_truncation_error.1
683                    / dt
684                    / lipschitz_constant.1,
685            )
686        };
687
688        // Numerically integrate equation
689        let mut t = 0.0;
690        let mut y = y0.clone();
691        let mut dy_storage = [
692            rhs(exact_solution(t - dt)),
693            rhs(exact_solution(t - 2.0 * dt)),
694        ];
695
696        for i in 0..100 {
697            let dy = rhs(y);
698            dy_storage[1] = dy_storage[0];
699            dy_storage[0] = dy;
700            y = adams_bashforth_2(y, dy_storage, dt, Vec2(0.0, 0.0)).unwrap();
701            t += dt;
702            let e = global_truncation_error(t);
703            let d1 = (y.0 - exact_solution(t).0).abs();
704            let d2 = (y.1 - exact_solution(t).1).abs();
705            if i > 0 {
706                assert!(d1 < e.0);
707                assert!(d2 < e.1);
708            }
709        }
710    }
711}
cellular_raza_core/backend/chili/solvers.rs

cellular_raza_core/backend/chili/
solvers.rs
