[mod] refactor

README.rst

@@ -41,9 +41,9 @@ Several URL parameters are optionally available:
 
 * ``rate``: draw rate (the smaller the smoother but the more demanding)
 
-The default configuration is equivalent to ``length = 192``, ``nitems = 8192``, and ``rate = 0.05``, namely:
+The default configuration is equivalent to ``length = 192``, ``nitems = 8192``, and ``rate = 0.1``, namely:
 
-``https://naokihori.github.io/Collision/index.html?length=192&nitems=8192&rate=0.05``.
+``https://naokihori.github.io/Collision/index.html?length=192&nitems=8192&rate=0.1``.
 
 The particle radii are fixed to ``0.5`` for now, and the restitution coefficient between particles is set to ``0.99``.
 I assume the domain is squared-shape and the periodic and wall-bounded conditions are imposed in the horizontal and the vertical directions, respectively.
@@ -53,24 +53,39 @@ Also, the number of particles is clipped if the volume fraction exceeds ``40%``.
 Method
 ******
 
+In this project, I aim at simulating finite-sized particles following Newtonian mechanics with overlaps prohibited.
+To this end, I need to properly detect all inter-particle collisions, which inherently requires O(N_p^2) operations, where N_p is the number of particles.
+It is necessary to reduce this cost somehow to treat say millions of particles.
+
+#. Event-driven approach
+
+   One possible way to handle inter-particle interactions is to introduce repulsive forces between particles, such that they repel to each other when too close.
+   Each particle motion is given by ordinary-differential equations, which can be solved by an appropriate time marcher (e.g. Runge-Kutta method).
+   The repellent force is, however, arbitrarily chosen and to be nicely designed.
+   To avoid these problems, I adopt the so-called event-driven approach with the cost of O(N_p) for the collision detection.
+   Moreover, this is an ODE-free method, which is advantageous to eliminate the numerical errors of the time-marching schemes.
+
 #. Cell method
 
-   All inter-particle collision events are considered, which requires essentially O(N_p^2) operations, where N_p is the number of particles.
-   To drop the cost, the so-called event-driven approach combined with the cell method is used, which leads to the cost of O(N_p^2 / N_c) where N_c is the number of cells.
+   It is desired only to consider particles near-by, which is achieved by the so-called cell method splitting the whole domain into many cells.
+   This leads to the cost of O(N_p^2 / N_c) where N_c is the number of cells.
+   Since N_c can be changed arbitrarily, the cost to detect collisions results in O(1).
+   Instead, cell-particle events (particles passing through the cell boundaries) are to be considered.
+   However, the cost to update events for each step remains O(1).
 
 #. Scheduler
 
-   Scheduling requests O(N_c) operations if implemented naively.
-   To reduce the cost, a minimum binary heap with the cost of O(log N_c) for insertions and deletions is adopted.
+   Although the events are separately queued for each cell, it is necessary to fetch the latest event among all cells, which requests O(N_c) operations if implemented naively (iterate over all cells, find the latest one).
+   To reduce the cost, a minimum binary heap with the cost of O(log N_c) to be balanced is adopted.
 
 #. Local time
 
-   Updating particle positions and velocities requests O(N_p) operations.
-   This is mitigated by introducing particle-based local time, so that particles are only synchronised when drawn.
+   Updating particle positions and velocities requests O(N_p) operations, and doing this process for each step is verbose.
+   This is mitigated by introducing particle-based local time, so that particles are only updated when they are involved.
 
 ***************
 Acknowledgement
 ***************
 
-I would like to thank Prof. Stefan Luding for his stimulating lecture in a JMBC course *Particle-based modeling techniques*.
+I would like to thank Prof. Stefan Luding for a stimulating lecture in a JMBC course *Particle-based modeling techniques*.
 

src/drawer.rs

@@ -50,19 +50,18 @@ impl Drawer {
         };
         let scale: f64 = screen_size / domain_size;
         // draw circles
+        const ARCS: [f64; 2] = [0., 2. * std::f64::consts::PI];
+        context.set_fill_style(&JsValue::from_str("#8888ff"));
+        context.begin_path();
         for particle in particles.iter() {
             let particle: Ref<Particle> = particle.borrow();
-            const ARCS: [f64; 2] = [0., 2. * std::f64::consts::PI];
-            let r: f64 = particle.rad;
-            let x: f64 = particle.pos[0];
-            let y: f64 = particle.pos[1];
-            context.set_fill_style(&convert(particle.val));
-            context.begin_path();
-            context
-                .arc(scale * x, scale * y, scale * r, ARCS[0], ARCS[1])
-                .unwrap();
-            context.fill();
+            let r: f64 = scale * particle.rad;
+            let x: f64 = scale * particle.pos[0];
+            let y: f64 = scale * particle.pos[1];
+            context.move_to(x, y);
+            context.arc(x, y, r, ARCS[0], ARCS[1]).unwrap();
         }
+        context.fill();
     }
 
     pub fn update_canvas_size(&self) {
@@ -73,50 +72,3 @@ impl Drawer {
         canvas.set_height(h as u32);
     }
 }
-
-fn convert(val: f64) -> JsValue {
-    let rgbcoefs: [[f64; 3]; 5] = [
-        [
-            2.672303238499781e-01,
-            5.015408860973969e-03,
-            3.290548382054911e-01,
-        ],
-        [
-            8.867281107764821e-01,
-            1.415434679048477e+00,
-            6.427369217396137e-01,
-        ],
-        [
-            -6.777660845884058e+00,
-            -8.089902514371242e-01,
-            2.998258532949060e+00,
-        ],
-        [
-            1.102198635856048e+01,
-            7.296293729490473e-01,
-            -9.057970794130403e+00,
-        ],
-        [
-            -4.404685706758277e+00,
-            -4.355228476501643e-01,
-            5.230151793650696e+00,
-        ],
-    ];
-    // fit polynomial
-    let mut rgb: [f64; 3] = [0., 0., 0.];
-    for (n, rgbcoef) in rgbcoefs.iter().enumerate() {
-        for m in 0..3 {
-            rgb[m] += rgbcoef[m] * val.powi(n as i32);
-        }
-    }
-    // truncate
-    for m in 0..3 {
-        rgb[m] = if rgb[m] < 0. { 0. } else { rgb[m] };
-        rgb[m] = if 1. < rgb[m] { 1. } else { rgb[m] };
-    }
-    let r: u8 = (255. * rgb[0]) as u8;
-    let g: u8 = (255. * rgb[1]) as u8;
-    let b: u8 = (255. * rgb[2]) as u8;
-    let string = format!("#{r:02x}{g:02x}{b:02x}");
-    JsValue::from_str(&string)
-}

src/lib.rs

@@ -14,9 +14,9 @@ pub struct Collision {
 
 #[wasm_bindgen::prelude::wasm_bindgen]
 impl Collision {
-    pub fn new(length: f64, nitems: usize, rate: f64) -> Collision {
+    pub fn new(length: f64, nitems: usize, rate: f64, seed: f64) -> Collision {
         let lengths: [f64; NDIMS] = [length, length];
-        let simulator = Simulator::new(rate, lengths, nitems);
+        let simulator = Simulator::new(rate, lengths, nitems, seed);
         let drawer = Drawer::new();
         Collision { simulator, drawer }
     }

src/main.rs

@@ -5,19 +5,20 @@ mod simulator;
 use crate::simulator::{Simulator, NDIMS};
 
 fn main() {
+    const SEED: f64 = 0.;
     let mut time: f64 = 0.;
-    let time_max: f64 = 100.;
+    let time_max: f64 = 200.;
     let sync_rate: f64 = 1.;
-    let lengths: [f64; NDIMS] = [32., 32.];
-    let nparticles: usize = 256;
-    let mut simulator = Simulator::new(sync_rate, lengths, nparticles);
+    let lengths: [f64; NDIMS] = [16., 16.];
+    let nparticles: usize = 128;
+    let mut simulator = Simulator::new(sync_rate, lengths, nparticles, SEED);
     let ncells: &[usize; NDIMS] = simulator.get_ncells();
     println!("number of cells: ({}, {})", ncells[0], ncells[1]);
     loop {
         simulator.integrate();
         time += sync_rate;
         println!("time: {:8.2e}", time);
-        if time_max < time {
+        if time_max <= time {
             break;
         }
     }

src/simulator.rs

@@ -28,11 +28,11 @@ pub struct Simulator {
 }
 
 impl Simulator {
-    pub fn new(sync_rate: f64, lengths: [f64; NDIMS], nparticles: usize) -> Simulator {
+    pub fn new(sync_rate: f64, lengths: [f64; NDIMS], nparticles: usize, seed: f64) -> Simulator {
         let time: f64 = 0.;
         let (ncells, cells): ([usize; NDIMS], Vec<Rc<RefCell<Cell>>>) = cell::init_cells(&lengths);
         let particles: Vec<Rc<RefCell<Particle>>> =
-            particle::init_particles(&lengths, &ncells, &cells, nparticles, time);
+            particle::init_particles(&lengths, &ncells, &cells, nparticles, time, seed);
         let mut scheduler = Scheduler::new(&cells);
         event::init_events(&lengths, &cells, &mut scheduler);
         Simulator {

src/simulator/cell.rs

@@ -9,7 +9,7 @@ use crate::simulator::NDIMS;
 /// Reference cell size.
 ///
 /// N.B.: This is a "typical" size of a cell and is not necessarily the exact size.
-///       For safety give 4 times larger than the radius of particles.
+///       For safety give more than 4 times larger than the radius of particles.
 const CELL_SIZE: f64 = 3.;
 
 /// Used to take into account the periodicity.
@@ -23,7 +23,7 @@ pub struct Cell {
     pub index: usize,
     pub bounds: [Extrema<f64>; NDIMS],
     pub particles: Rc<RefCell<Vec<Rc<RefCell<Particle>>>>>,
-    pub events: Rc<RefCell<Vec<Rc<RefCell<Event>>>>>,
+    pub events: Rc<RefCell<Vec<Event>>>,
     pub positions: [CellPosition; NDIMS],
     pub neighbours: [Extrema<usize>; NDIMS],
 }
@@ -171,7 +171,7 @@ pub fn init_cells(lengths: &[f64; NDIMS]) -> ([usize; NDIMS], Vec<Rc<RefCell<Cel
             },
         ];
         let particles = Rc::new(RefCell::new(Vec::<Rc<RefCell<Particle>>>::new()));
-        let events = Rc::new(RefCell::new(Vec::<Rc<RefCell<Event>>>::new()));
+        let events = Rc::new(RefCell::new(Vec::<Event>::new()));
         let positions: [CellPosition; NDIMS] = [
             if 0 == i {
                 CellPosition::NegativeEdge

src/simulator/debug.rs

@@ -6,38 +6,42 @@ use crate::simulator::particle::Particle;
 
 /// Checks if the particle knows which cells it belongs, and vice versa.
 #[allow(dead_code)]
-pub fn check_recognition(p: &Rc<RefCell<Particle>>, cell: &Ref<Cell>) {
+pub fn check_recognition(p: &Rc<RefCell<Particle>>, c0: &Rc<RefCell<Cell>>) {
     // particle -> cell check
-    let mut is_included: bool = false;
-    let qs: &Rc<RefCell<Vec<Rc<RefCell<Particle>>>>> = &cell.particles;
-    let qs: Ref<Vec<Rc<RefCell<Particle>>>> = qs.borrow();
-    for q in qs.iter() {
-        if Rc::ptr_eq(p, q) {
-            is_included = true;
-            break;
+    {
+        let mut is_included: bool = false;
+        let qs: &Rc<RefCell<Vec<Rc<RefCell<Particle>>>>> = &c0.borrow().particles;
+        let qs: Ref<Vec<Rc<RefCell<Particle>>>> = qs.borrow();
+        for q in qs.iter() {
+            if Rc::ptr_eq(p, q) {
+                is_included = true;
+                break;
+            }
+        }
+        if !is_included {
+            panic!(
+                "cell {} does not recognise it contains the particle {}",
+                c0.borrow().index,
+                p.borrow().index
+            );
         }
-    }
-    if !is_included {
-        panic!(
-            "cell {} does not recognise it contains the particle {}",
-            cell.index,
-            p.borrow().index
-        );
     }
     // cell -> particle check
-    let mut is_included: bool = false;
-    let cell_indices: &Vec<usize> = &p.borrow().cell_indices;
-    for &cell_index in cell_indices.iter() {
-        if cell.index == cell_index {
-            is_included = true;
-            break;
+    {
+        let mut is_included: bool = false;
+        let cells: &Vec<Rc<RefCell<Cell>>> = &p.borrow().cells;
+        for c1 in cells.iter() {
+            if Rc::ptr_eq(c0, c1) {
+                is_included = true;
+                break;
+            }
+        }
+        if !is_included {
+            panic!(
+                "particle {} does not recognise it belongs to the cell {}",
+                p.borrow().index,
+                c0.borrow().index
+            );
         }
-    }
-    if !is_included {
-        panic!(
-            "particle {} does not recognise it belongs to the cell {}",
-            p.borrow().index,
-            cell.index
-        );
     }
 }