moldyn_core/simulation/
sum.rs

1//! TODO document
2
3use crate::SimulationArgs;
4use crate::{Force, LennardJonesForce};
5use crate::{Particle, simulation::Simulation};
6use std::sync::Arc;
7
8/// The [DirectSum] simulation method is the most intuitive way to process
9/// a molecular dynamics simulation. It bases the computation on the
10/// [Direct Sum](https://en.wikipedia.org/wiki/Direct_sum) method.
11///
12/// **Newton Pair Optimization**: The only optimization [DirectSum] performs
13/// is avoiding computing the same pair of particles twice.
14pub struct DirectSum {
15    // TODO explain in slides why Arc works and Box does not
16    force: Arc<dyn Force>,
17    particles: Vec<Particle>,
18    args: SimulationArgs,
19}
20
21impl Simulation for DirectSum {
22    fn system_name(&self) -> &str {
23        "direct-sum"
24    }
25
26    fn particles(&self) -> &[Particle] {
27        &self.particles
28    }
29
30    fn particles_mut(&mut self) -> &mut [Particle] {
31        &mut self.particles
32    }
33
34    // index-based approach because two mutable iterators were problematic
35    fn for_each_particle_pairs_mut(&mut self, f: &mut dyn FnMut(&mut Particle, &mut Particle)) {
36        let count = self.particle_count();
37
38        for i in 0..count {
39            // newtons third law: skip same pairs
40            for j in (i + 1)..count {
41                // TODO maybe move this line one above, check efficiency of split_at_mut
42                // https://doc.rust-lang.org/std/vec/struct.Vec.html#method.split_at_mut
43                let (left, right) = self.particles.split_at_mut(j);
44
45                // conceptually:
46                //
47                // [p1,  p2,  p3,  p4,  p5]
48                //        i         j
49                // [p1,  p2,  p3],[p4,  p5]
50                //       ^^        ^^ avoid borrow issue with split_at_mut
51
52                f(&mut left[i], &mut right[0]);
53            }
54        }
55    }
56
57    fn particle_count(&self) -> usize {
58        self.particles.len()
59    }
60
61    fn add_particles(&mut self, particles: Vec<Particle>) {
62        self.particles.extend(particles);
63    }
64
65    fn get_force(&self) -> Arc<dyn Force> {
66        self.force.clone()
67    }
68
69    fn set_force(&mut self, force: Arc<dyn Force>) {
70        self.force = force;
71    }
72
73    fn args(&self) -> SimulationArgs {
74        self.args.clone()
75    }
76
77    fn set_args(&mut self, args: SimulationArgs) {
78        self.args = args;
79    }
80}
81
82impl Default for DirectSum {
83    fn default() -> Self {
84        Self {
85            force: Arc::new(LennardJonesForce::default()),
86            particles: Vec::new(),
87            args: SimulationArgs::default(),
88        }
89    }
90}
91
92// https://doc.rust-lang.org/nightly/unstable-book/library-features/test.html
93// build.rs is specifically to allow this type of benchmarks
94#[cfg(all(test, nightly))]
95mod benchmark {
96    use crate::{
97        CustomForce, DirectSum, LennardJonesForce, NewtonForce, Particle, Simulation,
98        SimulationArgs, Vec3,
99    };
100    use meval::Expr;
101    use std::sync::Arc;
102    use test::Bencher;
103
104    #[bench]
105    fn halleys_comet(b: &mut Bencher) {
106        let particles = vec![
107            // sun
108            Particle::from_data(Vec3::zero(), Vec3::zero(), 1.0),
109            // earth
110            Particle::from_data(Vec3::new(0.0, 1.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), 3.0e-6),
111            // jupiter
112            Particle::from_data(
113                Vec3::new(0.0, 5.36, 0.0),
114                Vec3::new(-0.425, 0.0, 0.0),
115                9.55e-4,
116            ),
117            // halleys comet
118            Particle::from_data(
119                Vec3::new(34.75, 0.0, 0.0),
120                Vec3::new(0.0, 0.0296, 0.0),
121                1.0e-14,
122            ),
123        ];
124
125        let mut simulation = DirectSum::default();
126        simulation.set_force(Arc::new(NewtonForce::default()));
127        simulation.add_particles(particles);
128
129        b.iter(|| {
130            simulation.step(0.01);
131        });
132    }
133
134    #[bench]
135    fn halleys_comet_lennard_jones(b: &mut Bencher) {
136        let particles = vec![
137            // sun
138            Particle::from_data(Vec3::zero(), Vec3::zero(), 1.0),
139            // earth
140            Particle::from_data(Vec3::new(0.0, 1.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), 3.0e-6),
141            // jupiter
142            Particle::from_data(
143                Vec3::new(0.0, 5.36, 0.0),
144                Vec3::new(-0.425, 0.0, 0.0),
145                9.55e-4,
146            ),
147            // halleys comet
148            Particle::from_data(
149                Vec3::new(34.75, 0.0, 0.0),
150                Vec3::new(0.0, 0.0296, 0.0),
151                1.0e-14,
152            ),
153        ];
154
155        let mut simulation = DirectSum::default();
156        simulation.set_force(Arc::new(LennardJonesForce::default()));
157        simulation.add_particles(particles);
158
159        b.iter(|| {
160            simulation.step(0.01);
161        });
162    }
163
164    #[bench]
165    fn halleys_comet_custom_newton_force(b: &mut Bencher) {
166        let particles = vec![
167            // sun
168            Particle::from_data(Vec3::zero(), Vec3::zero(), 1.0),
169            // earth
170            Particle::from_data(Vec3::new(0.0, 1.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), 3.0e-6),
171            // jupiter
172            Particle::from_data(
173                Vec3::new(0.0, 5.36, 0.0),
174                Vec3::new(-0.425, 0.0, 0.0),
175                9.55e-4,
176            ),
177            // halleys comet
178            Particle::from_data(
179                Vec3::new(34.75, 0.0, 0.0),
180                Vec3::new(0.0, 0.0296, 0.0),
181                1.0e-14,
182            ),
183        ];
184
185        let mut simulation = DirectSum::default();
186        let force_expr = CustomForce::from_expr("M / r");
187        let force = simulation.set_force(Arc::new(NewtonForce::default()));
188        simulation.add_particles(particles);
189
190        b.iter(|| {
191            simulation.step(0.01);
192        });
193    }
194
195    #[bench]
196    fn halleys_comet_custom_lennard_jones(b: &mut Bencher) {
197        let particles = vec![
198            // sun
199            Particle::from_data(Vec3::zero(), Vec3::zero(), 1.0),
200            // earth
201            Particle::from_data(Vec3::new(0.0, 1.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), 3.0e-6),
202            // jupiter
203            Particle::from_data(
204                Vec3::new(0.0, 5.36, 0.0),
205                Vec3::new(-0.425, 0.0, 0.0),
206                9.55e-4,
207            ),
208            // halleys comet
209            Particle::from_data(
210                Vec3::new(34.75, 0.0, 0.0),
211                Vec3::new(0.0, 0.0296, 0.0),
212                1.0e-14,
213            ),
214        ];
215
216        // let frac = self.sigma / distance;
217        //     let frac6 = frac.powi(6);
218        //     let frac12 = frac6.powi(2);
219
220        //     4.0 * self.epsilon * (frac12 - frac6)
221
222        let mut simulation = DirectSum::default();
223        let force_expr = CustomForce::from_expr("20.0 * ((1 / r)^12 - (1 / r)^6)");
224        let force = simulation.set_force(Arc::new(NewtonForce::default()));
225        simulation.add_particles(particles);
226
227        b.iter(|| {
228            simulation.step(0.01);
229        });
230    }
231
232    #[bench]
233    fn halleys_comet_custom_heavy(b: &mut Bencher) {
234        let particles = vec![
235            // sun
236            Particle::from_data(Vec3::zero(), Vec3::zero(), 1.0),
237            // earth
238            Particle::from_data(Vec3::new(0.0, 1.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), 3.0e-6),
239            // jupiter
240            Particle::from_data(
241                Vec3::new(0.0, 5.36, 0.0),
242                Vec3::new(-0.425, 0.0, 0.0),
243                9.55e-4,
244            ),
245            // halleys comet
246            Particle::from_data(
247                Vec3::new(34.75, 0.0, 0.0),
248                Vec3::new(0.0, 0.0296, 0.0),
249                1.0e-14,
250            ),
251        ];
252
253        let mut simulation = DirectSum::default();
254        let force_expr = CustomForce::from_expr("(1 / r)^12 - (1 / r)^11 + (1 / r)^10 - (1 / r)^9 + (1 / r)^6 - (1 / r)^5");
255        let force = simulation.set_force(Arc::new(NewtonForce::default()));
256        simulation.add_particles(particles);
257
258        b.iter(|| {
259            simulation.step(0.01);
260        });
261    }
262
263    #[bench]
264    fn ten_bodies(b: &mut Bencher) {
265        let mut particles = vec![];
266
267        for x in 0..10 {
268            particles.push(Particle::at(x as f64, 0.0, 0.0).with_mass(1.0));
269        }
270
271        let mut simulation = DirectSum::default();
272        simulation.set_force(Arc::new(NewtonForce::default()));
273        simulation.add_particles(particles);
274
275        b.iter(|| {
276            simulation.step(0.01);
277        });
278    }
279
280    #[bench]
281    fn ten_bodies_lennard_jones(b: &mut Bencher) {
282        let mut particles = vec![];
283
284        for x in 0..10 {
285            particles.push(Particle::at(x as f64, 0.0, 0.0).with_mass(1.0));
286        }
287
288        let mut simulation = DirectSum::default();
289        simulation.set_force(Arc::new(LennardJonesForce::default()));
290        simulation.add_particles(particles);
291
292        b.iter(|| {
293            simulation.step(0.01);
294        });
295    }
296
297    #[bench]
298    fn hundred_bodies(b: &mut Bencher) {
299        let mut particles = vec![];
300
301        for x in 0..10 {
302            for y in 0..10 {
303                particles.push(Particle::at(x as f64, y as f64, 0.0).with_mass(1.0));
304            }
305        }
306
307        let mut simulation = DirectSum::default();
308        simulation.set_force(Arc::new(NewtonForce::default()));
309        simulation.add_particles(particles);
310
311        b.iter(|| {
312            simulation.step(0.01);
313        });
314    }
315
316    #[bench]
317    fn hundred_bodies_lennard_jones(b: &mut Bencher) {
318        let mut particles = vec![];
319
320        for x in 0..10 {
321            for y in 0..10 {
322                particles.push(Particle::at(x as f64, y as f64, 0.0).with_mass(1.0));
323            }
324        }
325
326        let mut simulation = DirectSum::default();
327        simulation.set_force(Arc::new(LennardJonesForce::default()));
328        simulation.add_particles(particles);
329
330        b.iter(|| {
331            simulation.step(0.01);
332        });
333    }
334
335    #[bench]
336    fn thousand_bodies(b: &mut Bencher) {
337        let mut particles = vec![];
338
339        for x in 0..10 {
340            for y in 0..10 {
341                for z in 0..10 {
342                    particles.push(Particle::at(x as f64, y as f64, 0.0).with_mass(1.0));
343                }
344            }
345        }
346
347        let mut simulation = DirectSum::default();
348        simulation.set_force(Arc::new(NewtonForce::default()));
349        simulation.add_particles(particles);
350
351        b.iter(|| {
352            simulation.step(0.01);
353        });
354    }
355
356    #[bench]
357    fn thousand_bodies_lennard_jones(b: &mut Bencher) {
358        let mut particles = vec![];
359
360        for x in 0..10 {
361            for y in 0..10 {
362                for z in 0..10 {
363                    particles.push(Particle::at(x as f64, y as f64, 0.0).with_mass(1.0));
364                }
365            }
366        }
367
368        let mut simulation = DirectSum::default();
369        simulation.set_force(Arc::new(LennardJonesForce::default()));
370        simulation.add_particles(particles);
371
372        b.iter(|| {
373            simulation.step(0.01);
374        });
375    }
376}
moldyn_core/simulation/sum.rs

moldyn_core/simulation/
sum.rs
