Benchmarking

* Add benchmark simulation
* Precalculate some `AtomicTransition` numbers.
* Remove all barriers.
* Sampler initialisation done by using to_vec to copy blank initialised vector.
* inner parallelism on most resource intensive systems.

examples/benchmark.rs

@@ -0,0 +1,187 @@
+//! Simulation of atoms cooled to the Doppler limit.
+
+extern crate magneto_optical_trap as lib;
+extern crate nalgebra;
+use lib::atom::{Atom, AtomicTransition, Force, Mass, Position, Velocity};
+use lib::ecs;
+use lib::initiate::NewlyCreated;
+use lib::integrator::Timestep;
+use lib::laser::cooling::CoolingLight;
+use lib::laser::force::ApplyEmissionForceOption;
+use lib::laser::gaussian::GaussianBeam;
+use lib::laser::photons_scattered::EnableScatteringFluctuations;
+use lib::magnetic::quadrupole::QuadrupoleField3D;
+use lib::output::file;
+use lib::output::file::Text;
+use nalgebra::Vector3;
+use rand::distributions::{Distribution, Normal};
+use specs::{Builder, World};
+use std::time::Instant;
+
+fn main() {
+    let now = Instant::now();
+
+    // Create the simulation world and builder for the ECS dispatcher.
+    let mut world = World::new();
+    ecs::register_components(&mut world);
+    ecs::register_resources(&mut world);
+    let mut builder = ecs::create_simulation_dispatcher_builder();
+
+    // Configure thread pool.
+    let pool = rayon::ThreadPoolBuilder::new()
+        .num_threads(12)
+        .build()
+        .unwrap();
+
+    builder.add_pool(::std::sync::Arc::new(pool));
+
+    let mut dispatcher = builder.build();
+    dispatcher.setup(&mut world.res);
+
+    // Create magnetic field.
+    world
+        .create_entity()
+        .with(QuadrupoleField3D::gauss_per_cm(18.2, Vector3::z()))
+        .with(Position {
+            pos: Vector3::new(0.0, 0.0, 0.0),
+        })
+        .build();
+
+    // Create cooling lasers.
+    let detuning = -3.0;
+    let power = 0.02;
+    let radius = 66.7e-3 / (2.0_f64.sqrt());
+    let beam_centre = Vector3::new(0.0, 0.0, 0.0);
+
+    world
+        .create_entity()
+        .with(GaussianBeam {
+            intersection: beam_centre.clone(),
+            e_radius: radius,
+            power: power,
+            direction: Vector3::new(0.0, 0.0, 1.0),
+        })
+        .with(CoolingLight::for_species(
+            AtomicTransition::rubidium(),
+            detuning,
+            -1,
+        ))
+        .build();
+    world
+        .create_entity()
+        .with(GaussianBeam {
+            intersection: beam_centre.clone(),
+            e_radius: radius,
+            power: power,
+            direction: Vector3::new(0.0, 0.0, -1.0),
+        })
+        .with(CoolingLight::for_species(
+            AtomicTransition::rubidium(),
+            detuning,
+            -1,
+        ))
+        .build();
+    world
+        .create_entity()
+        .with(GaussianBeam {
+            intersection: beam_centre.clone(),
+            e_radius: radius,
+            power: power,
+            direction: Vector3::new(-1.0, 0.0, 0.0),
+        })
+        .with(CoolingLight::for_species(
+            AtomicTransition::rubidium(),
+            detuning,
+            1,
+        ))
+        .build();
+    world
+        .create_entity()
+        .with(GaussianBeam {
+            intersection: beam_centre.clone(),
+            e_radius: radius,
+            power: power,
+            direction: Vector3::new(1.0, 0.0, 0.0),
+        })
+        .with(CoolingLight::for_species(
+            AtomicTransition::rubidium(),
+            detuning,
+            1,
+        ))
+        .build();
+    world
+        .create_entity()
+        .with(GaussianBeam {
+            intersection: beam_centre.clone(),
+            e_radius: radius,
+            power: power,
+            direction: Vector3::new(0.0, 1.0, 0.0),
+        })
+        .with(CoolingLight::for_species(
+            AtomicTransition::rubidium(),
+            detuning,
+            1,
+        ))
+        .build();
+    world
+        .create_entity()
+        .with(GaussianBeam {
+            intersection: beam_centre.clone(),
+            e_radius: radius,
+            power: power,
+            direction: Vector3::new(0.0, -1.0, 0.0),
+        })
+        .with(CoolingLight::for_species(
+            AtomicTransition::rubidium(),
+            detuning,
+            1,
+        ))
+        .build();
+
+    // Define timestep
+    world.add_resource(Timestep { delta: 1.0e-6 });
+
+    let vel_dist = Normal::new(0.0, 0.22);
+    let pos_dist = Normal::new(0.0, 1.2e-4);
+    let mut rng = rand::thread_rng();
+
+    // Add atoms
+    for _ in 0..10000 {
+        world
+            .create_entity()
+            .with(Position {
+                pos: Vector3::new(
+                    pos_dist.sample(&mut rng),
+                    pos_dist.sample(&mut rng),
+                    pos_dist.sample(&mut rng),
+                ),
+            })
+            .with(Velocity {
+                vel: Vector3::new(
+                    vel_dist.sample(&mut rng),
+                    vel_dist.sample(&mut rng),
+                    vel_dist.sample(&mut rng),
+                ),
+            })
+            .with(Force::new())
+            .with(Mass { value: 87.0 })
+            .with(AtomicTransition::rubidium())
+            .with(Atom)
+            .with(NewlyCreated)
+            .build();
+    }
+
+    // Enable fluctuation options
+    //  * Allow photon numbers to fluctuate.
+    //  * Allow random force from emission of photons.
+    world.add_resource(ApplyEmissionForceOption {});
+    world.add_resource(EnableScatteringFluctuations {});
+
+    // Run the simulation for a number of steps.
+    for _i in 0..5000 {
+        dispatcher.dispatch(&mut world.res);
+        world.maintain();
+    }
+
+    println!("Simulation completed in {} ms.", now.elapsed().as_millis());
+}

examples/doppler_limit.rs

@@ -27,14 +27,6 @@ fn main() {
     ecs::register_resources(&mut world);
     let mut builder = ecs::create_simulation_dispatcher_builder();
 
-    // Configure thread pool.
-    let pool = rayon::ThreadPoolBuilder::new()
-        .num_threads(6)
-        .build()
-        .unwrap();
-
-    builder.add_pool(::std::sync::Arc::new(pool));
-
     // Configure simulation output.
     builder = builder.with(
         file::new::<Velocity, Text>("vel.txt".to_string(), 10),
@@ -153,7 +145,7 @@ fn main() {
     let mut rng = rand::thread_rng();
 
     // Add atoms
-    for _ in 0..10000 {
+    for _ in 0..1000 {
         world
             .create_entity()
             .with(Position {

src/atom.rs

@@ -1,6 +1,6 @@
 //! Common atom components and systems.
 
-use crate::constant::{BOHRMAG, C};
+use crate::constant::{BOHRMAG, C, HBAR, PI};
 use crate::output::file::BinaryConversion;
 use nalgebra::Vector3;
 use serde::{Deserialize, Serialize};
@@ -132,12 +132,21 @@ pub struct AtomicTransition {
 	pub linewidth: f64,
 	/// Saturation intensity, in units of W/m^2.
 	pub saturation_intensity: f64,
+	/// Precalculate prefactor used in the determination of rate coefficients.
+	pub rate_prefactor: f64,
 }
 
 impl Component for AtomicTransition {
 	type Storage = VecStorage<Self>;
 }
 impl AtomicTransition {
+	pub fn calculate(mut self) -> Self {
+		self.rate_prefactor =
+			3. / 4. * C * C / HBAR * (self.frequency).powf(-3.) * self.linewidth * self.linewidth
+				/ 2. / PI;
+		self
+	}
+
 	/// Creates an `AtomicTransition` component populated with parameters for Rubidium.
 	/// The parameters are taken from Daniel Steck's Data sheet on Rubidium-87.
 	pub fn rubidium() -> Self {
@@ -148,7 +157,9 @@ impl AtomicTransition {
 			frequency: C / 780.0e-9,
 			linewidth: 6.065e6,          // [Steck, Rubidium87]
 			saturation_intensity: 16.69, // [Steck, Rubidium 87, D2 cycling transition]
+			rate_prefactor: 0.0,         // set in calculate
 		}
+		.calculate()
 	}
 
 	/// Creates an `AtomicTransition` component populated with parameters for Strontium.
@@ -161,7 +172,9 @@ impl AtomicTransition {
 			frequency: 650_759_219_088_937.,
 			linewidth: 32e6,             // [Nosske2017]
 			saturation_intensity: 430.0, // [Nosske2017, 43mW/cm^2]
+			rate_prefactor: 0.0,         // set in calculate
 		}
+		.calculate()
 	}
 
 	/// Creates an `AtomicTransition` component populated with parameters for red Strontium transition.
@@ -174,7 +187,9 @@ impl AtomicTransition {
 			frequency: 434_829_121_311_000., // NIST, doi:10.1063/1.344917
 			linewidth: 7_400.,               // [Schreck2013]
 			saturation_intensity: 0.0295,    // [SChreck2013, 3 µW/cm^2]
+			rate_prefactor: 0.0,             // set in calculate
 		}
+		.calculate()
 	}
 
 	/// Creates an `AtomicTransition` component populated with parameters for Erbium.
@@ -186,7 +201,9 @@ impl AtomicTransition {
 			frequency: 5.142e14,
 			linewidth: 190e3,
 			saturation_intensity: 0.13,
+			rate_prefactor: 0.0, // set in calculate
 		}
+		.calculate()
 	}
 	/// Creates an `AtomicTransition` component populated with parameters for Erbium 401 .
 	pub fn erbium_401() -> Self {
@@ -197,7 +214,9 @@ impl AtomicTransition {
 			frequency: 7.476e14,
 			linewidth: 30e6,
 			saturation_intensity: 56.0,
+			rate_prefactor: 0.0, // set in calculate
 		}
+		.calculate()
 	}
 
 	pub fn gamma(&self) -> f64 {

src/ecs.rs

@@ -42,13 +42,13 @@ pub fn create_simulation_dispatcher_builder() -> DispatcherBuilder<'static, 'sta
 	builder = builder.with(LargerEarlyTimestepOptimizationSystem, "opt", &[]);
 	builder = builder.with(ClearForceSystem, "clear", &[]);
 	builder = builder.with(DeflagNewAtomsSystem, "deflag", &[]);
-	builder.add_barrier();
+	//builder.add_barrier();
 	builder = magnetic::add_systems_to_dispatch(builder, &[]);
-	builder.add_barrier();
+	//builder.add_barrier();
 	builder = laser::add_systems_to_dispatch(builder, &[]);
-	builder.add_barrier();
+	//builder.add_barrier();
 	builder = atom_sources::add_systems_to_dispatch(builder, &[]);
-	builder.add_barrier();
+	//builder.add_barrier();
 	builder = builder.with(ApplyGravitationalForceSystem, "add_gravity", &["clear"]);
 	builder = builder.with(
 		EulerIntegrationSystem,

src/laser/doppler.rs

@@ -96,14 +96,17 @@ impl<'a> System<'a> for InitialiseDopplerShiftSamplersSystem {
         WriteStorage<'a, DopplerShiftSamplers>,
     );
     fn run(&mut self, (cooling, cooling_index, mut samplers): Self::SystemData) {
+        use rayon::prelude::*;
+        use specs::ParJoin;
+
         let mut content = Vec::new();
         for (_, _) in (&cooling, &cooling_index).join() {
             content.push(DopplerShiftSampler::default());
         }
 
-        for mut sampler in (&mut samplers).join() {
-            sampler.contents = content.clone();
-        }
+        (&mut samplers).par_join().for_each(|mut sampler| {
+            sampler.contents = content.to_vec();
+        });
     }
 }