Space

Simulation Loop and Data Flow

The simulation uses a multi-threaded architecture to ensure a responsive user interface. The core physics calculations run in a separate thread, communicating with the main UI thread via Qt’s signal and slot mechanism. This prevents the UI from freezing while the physics engine is busy.

The diagram below illustrates the data flow between the main components:

Step-by-step breakdown:

1. The PhysicsWorker calls Space.compute_dynamics() to run one full step of the simulation (integration, force calculation, collision handling).

2. Once the step is complete, the new positions and orientations are returned to the PhysicsWorker.

3. The worker emits the updated_body_data signal, passing the new data as NumPy arrays. This is a thread-safe way to communicate with the UI thread.

4. The Display class’s _update_frame slot, running in the UI thread, receives this signal. It updates the visual representation of the bodies and then uses QTimer.singleShot(0) to schedule the next physics step. Using a zero-delay timer posts an event to the Qt event loop, ensuring the UI remains responsive and updates are not queued faster than they can be processed.

5. The Qt event loop processes the timer event, which calls the perform_single_step slot on the PhysicsWorker, starting the cycle over again.

This architecture effectively decouples the rendering rate from the simulation rate, leading to a smooth user experience.

Adaptive Time-Stepping for Stability and Performance

The simulation employs an adaptive time-stepping algorithm to ensure both stability and efficiency. A fixed, large time step could cause bodies to gain too much velocity in a single frame or even pass through each other (an effect known as “tunnelling”). A fixed, small time step would be safe but computationally wasteful when bodies are far apart and interacting weakly.

The adaptive algorithm, implemented in the _update_adaptive_dt method, calculates an optimal dt for each frame based on the state of the system in the previous frame. This is achieved by considering two primary constraints:

1. Acceleration-Based Constraint (Positional Accuracy)

   This constraint is designed to maintain positional accuracy. It prevents any single body from moving too far in one step due to high acceleration. The time step is calculated based on the maximum acceleration observed in the system, using the kinematic equation for displacement under constant acceleration: \(d = \frac{1}{2} a t^2\).

   The goal is to limit the maximum distance traveled to a small target value, TARGET_DX. The time step is therefore derived as:

   \[dt_a^2 = \frac{2 \times \text{TARGET_DX}}{\left\| a_{max} \right\|}\]

2. Velocity-Based Constraint (Tunnelling Prevention)

   This constraint prevents tunnelling. It is a form of the Courant-Friedrichs-Lewy (CFL) condition, which ensures that a body does not travel more than a fraction of its own size in a single time step. The time step is calculated based on the fastest-moving body in the system:

   \[dt_v^2 = \left( 0.5 \times \frac{\text{radius}_{\text{fastest}}}{\left\| v_{max} \right\|} \right)^2\]

Combining and Smoothing

The final target time step for the frame is the minimum of the two calculated values, capped by a global maximum MAX_DT:

\[\text{target_dt} = \sqrt{\min(\text{MAX_DT}^2, dt_a^2, dt_v^2)}\]

To prevent rapid, jerky oscillations in the simulation speed, the change in dt is smoothed using an exponential moving average (a low-pass filter). This ensures a more stable and visually pleasing simulation.

Collision Detection and Response

The simulation uses a sophisticated, multi-stage process to handle collisions efficiently and robustly.

1. Collision Detection

Detecting collisions in a system with many bodies can be computationally expensive. A naive approach would check every body against every other body (an \(O(n^2)\) operation). To avoid this, a two-stage process is used in _build_collision_graph:

• Broad Phase: The goal is to quickly eliminate pairs of bodies that are definitely not colliding. This is achieved using a k-d tree (scipy.spatial.cKDTree), a data structure that allows for very fast spatial queries. By querying the tree for all pairs of bodies that are within a certain distance of each other, we generate a much smaller list of “candidate pairs”.

• Narrow Phase: For the candidate pairs identified in the broad phase, a precise, vectorized intersection test is performed. This check calculates the exact distance between the two bodies in each pair and compares it to the sum of their radii.

The result of this process is a collision graph, where each body is a node and an edge exists between any two bodies that are currently colliding.

2. Collision Response

Once collisions are detected, they must be resolved. A key challenge is that resolving one collision (e.g., A-B) can affect the velocity of bodies involved in other collisions (e.g., A-C). Handling these dependencies incorrectly can lead to instability. The simulation solves this with a graph-based approach in _resolve_all_collisions:

• Graph Colouring: The collision graph is “coloured” using a greedy algorithm (_colour_graph). This assigns a colour (an integer) to each body such that no two adjacent bodies in the graph share the same colour. This partitions the colliding bodies into independent sets. All bodies of the same colour can have their collisions resolved simultaneously without causing conflicts.

• Staged Resolution: The simulation iterates through each colour, resolving all collisions associated with that colour in a single stage. This ensures that the resolution of one collision does not interfere with another in the same stage.

• Interpenetration Resolution: Due to discrete time steps, bodies may slightly overlap before a collision is detected. The _resolve_interpenetration method applies a small positional correction to push them apart. This is crucial for stability, as it prevents bodies from sinking into each other and causing unrealistic forces.

• Impulse-Based Physics: The change in velocity is calculated using an impulse-based model, which is well-suited for rigid body dynamics. The total impulse is composed of two parts:

  1. Restitution (Normal) Impulse: This impulse acts along the collision normal (the line connecting the centres of the two bodies). It is responsible for the “bounciness” of the collision and is governed by the coeff_restitution.

  2. Friction Impulse: This impulse acts tangentially to the collision surface and opposes the relative sliding motion. The model handles both static and kinetic friction by clamping the friction impulse to the limit defined by Coulomb’s law (\(J_t \le \mu_s J_n\)).

This multi-stage, impulse-based approach ensures that collisions are handled in a way that is both physically realistic and numerically stable.

Graph Colouring for Collision Resolution

The simulation uses a graph colouring algorithm to resolve collisions safely and efficiently. In each time step, the simulation identifies pairs of bodies that are colliding. To handle these collisions without introducing artificial biases or instabilities, the collisions are processed in a specific order determined by graph colouring.

1. Collision Graph Construction

The first step is to construct a collision graph, where each node represents a body, and an edge connects two nodes if the corresponding bodies are in contact.

2. Graph Colouring Algorithm

A graph colouring algorithm assigns a “colour” to each node (body) in the collision graph such that no two adjacent nodes (colliding bodies) have the same colour. This partitioning ensures that all bodies of the same colour can be processed simultaneously without conflicts. The simulation uses a greedy colouring algorithm:

1. Iterate through each body in the simulation.

2. For each body, identify the colours used by its colliding neighbours.

3. Assign the smallest available colour (non-negative integer) that is not used by any of its neighbours.

3. Collision Resolution Stages

The collision resolution process is divided into stages, one for each colour in the graph. In each stage, all collisions between bodies of the same colour are resolved. Because no two bodies of the same colour are in contact, these collisions can be processed independently and, in principle, in parallel.

Benefits of Graph Colouring

• Parallelism: Collisions within the same colour group can be processed independently, allowing for parallel execution and improved performance.

• Stability: By resolving collisions in stages based on graph colouring, the simulation avoids conflicts and artificial biases that can arise from processing collisions in an arbitrary order.

• Efficiency: The greedy colouring algorithm is relatively simple and efficient, adding minimal overhead to the collision resolution process.

The number of colours required to colour the collision graph depends on the complexity of the collision scenario. In practice, the number of colours is typically small, allowing for efficient collision resolution.

Collision Response Model

Once collisions are detected and partitioned using graph colouring, they must be resolved. The simulation uses an impulse-based model, which is standard for rigid body physics, as it correctly handles the instantaneous changes in velocity and angular velocity that occur during a collision.

The core idea is to calculate a total impulse vector, \(J\), and apply it to the colliding bodies to change their momentum. This impulse is composed of two orthogonal components: the normal impulse and the tangential (friction) impulse.

1. Normal Impulse (Restitution and Positional Correction)

The normal impulse acts along the collision normal (the line connecting the centres of the two bodies). It is responsible for two things:

• Restitution (Bounciness): It reverses the relative velocity along the normal, scaled by the coeff_restitution. A value of 1.0 would be a perfectly elastic collision, while 0.0 would be perfectly inelastic.

• Positional Correction (Baumgarte Stabilization): Due to discrete time steps, bodies can slightly overlap (interpenetrate). If not corrected, this leads to instability. The impulse calculation includes a bias term that creates a small corrective push to resolve this overlap over the next frame. This is a robust and stable method for preventing bodies from sinking into each other.

The final normal impulse, \(J_n\), combines both the restitution and the positional bias, ensuring that bodies both bounce off each other and are pushed apart correctly.

2. Friction Impulse

The friction impulse, \(J_t\), acts tangentially to the collision surface and opposes the relative sliding motion between the two bodies. The simulation implements a Coulomb friction model:

• Static Friction: If the tangential sliding velocity is low, the friction impulse is calculated to be exactly what is needed to make the relative tangential velocity zero, provided this impulse is less than the static friction limit (\(J_t \le \mu_s J_n\)).

• Kinetic Friction: If the impulse required to stop the sliding exceeds the static friction limit, it means the bodies are actively sliding. In this case, the friction impulse is clamped to the kinetic friction limit (\(J_t = \mu_k J_n\)), acting in the direction opposite to the sliding.

Applying the Impulse

The total impulse, \(J = J_n + J_t\), is applied to each body. The apply_impulse method in the BodyProxy class correctly updates both the linear and angular velocity of the body based on the impulse-momentum theorem:

• Linear Velocity Change: \(\Delta v = J / m\)

• Angular Velocity Change: \(\Delta \omega = I^{-1} (r \times J)\), where \(r\) is the vector from the body’s centre of mass to the point of contact.

This staged, impulse-based approach ensures that collisions are handled in a way that is both physically realistic and numerically stable.

Initial Position Generation

Generating the initial positions for all bodies presents a challenge: they must not overlap. A simple approach of placing bodies one by one and checking for collisions becomes extremely inefficient as the number of bodies or the density increases, as finding a valid spot for the last few bodies can require many attempts.

The simulation solves this with a fast, vectorized, iterative filtering approach in _generate_non_overlapping_positions.

The Algorithm

1. Generate Excess Candidates: Instead of generating exactly n positions, the algorithm starts by creating a larger number of candidate positions (e.g., n * 1.2) randomly distributed in the space.

2. Iterative Filtering and Replacement: The algorithm then enters a loop that iteratively refines this set of candidates:

   1. A k-d tree is built from the current candidate positions.

   2. A broad-phase check quickly identifies all pairs of candidates that are potentially overlapping.

   3. A narrow-phase check confirms which of these candidates are actually colliding.

   4. All colliding candidates are removed from the set.

   5. The set is replenished with an equal number of new, randomly generated candidates to maintain the size of the pool.

3. Termination: This loop continues until an iteration occurs where no overlaps are found, or a maximum number of iterations is reached. The method then returns the first n positions from the clean set.

This “generate and filter” strategy is significantly more performant than a one-by-one placement strategy because it leverages highly optimized, vectorized NumPy and SciPy operations to process all bodies at once.

Simulation Warm-up

A crucial step for ensuring a stable start to the simulation is the “warm-up” phase, implemented in the _warm_up_simulation method.

The Problem

The simulation begins with all bodies having zero initial acceleration. The adaptive time-stepping algorithm calculates the time step for the current frame based on the accelerations from the previous frame. On the very first frame, it would see zero acceleration and incorrectly choose the maximum allowed dt. If bodies start close together, this large initial time step could cause them to gain enormous velocities in a single step, leading to an unstable “explosion” that violates energy conservation.

The Solution

To prevent this, the simulation runs a few “warm-up” iterations before the main loop begins. During this phase: 1. A very small, safe time step is used to kickstart the process. 2. The core dynamics calculation (compute_dynamics) is called several times. 3. This populates the acceleration arrays with realistic values based on the initial positions. 4. The adaptive time-stepper then converges to a safe, stable dt value.

This pre-calculation ensures that the first visible frame of the simulation begins with a sensible time step, preventing numerical instability from the start.

Energy Conservation and Stability Check

A fundamental principle of a closed physical system is the conservation of energy. While perfect conservation is impossible in a discrete numerical simulation, the total energy should remain nearly constant. A significant, non-physical increase in energy is a clear sign of numerical instability.

The simulation includes a built-in diagnostic check, implemented in the _log_system_energy method, to monitor the system’s total energy over time. This check runs periodically (controlled by the --log-freq command-line argument) and calculates the three components of mechanical energy:

• Translational Kinetic Energy: \(\sum \frac{1}{2} m_i v_i^2\)

• Rotational Kinetic Energy: \(\sum \frac{1}{2} I_i \omega_i^2\)

• Gravitational Potential Energy: \(\sum_{i<j} -G \frac{m_i m_j}{r_{ij}}\)

The total energy is logged to the console if diagnostics are enabled.

Diagnosing Instability

If the total energy increases by more than a small tolerance between checks, a warning is logged. An energy increase is a strong indicator that the collision response is adding energy to the system, often due to complex, multi-body pile-ups.

To aid in debugging, the warning message automatically includes the state (velocities and angular velocities) of the bodies that were colliding in the previous frame. Since the collision response is the most likely source of instability, this information allows a developer to immediately inspect the conditions that led to the energy gain and diagnose the issue.

Module Overview

The core physics engine of the simulation.

This module provides the Space class, which manages the simulation loop, including gravity calculations, collision detection, and collision response for all bodies.

module:

space

author:

Le Bars, Yoann

class src.space.CollisionContext(b1: BodyProxy, b2: BodyProxy, n_vec: ndarray, r1_vec: ndarray, r2_vec: ndarray, v_rel: ndarray, overlap: float)

Bases: object

A data container to hold all relevant information for a single collision event. This simplifies passing collision data between functions.

Variables:

• b1 (BodyProxy) – First body.

• b2 (BodyProxy) – Second body.

• n_vec (np.ndarray) – The collision normal vector.

• r1_vec (np.ndarray) – The relative position vector from body 1.

• r2_vec (np.ndarray) – The relative position vector from body 2.

• v_rel (np.ndarray) – The relative velocity vector.

• overlap (float) – The amount of overlap.

b1: BodyProxy

b2: BodyProxy

classmethod from_bodies(b1: BodyProxy, b2: BodyProxy) → CollisionContext | None

Factory method to create a CollisionContext if two bodies are colliding.

Parameters:

• b1 (BodyProxy) – First body.

• b2 (BodyProxy) – Second body.

Returns:

A CollisionContext object if a collision occurs, otherwise None.

Return type:

CollisionContext | None

n_vec: ndarray

overlap: float

r1_vec: ndarray

r2_vec: ndarray

v_rel: ndarray

class src.space.Space(config: SimulationConfig, initial_dt: float)

Bases: object

Class describing a space in which move several bodies.

Variables:

• _dt (float) – Model time step (in s).

• _previous_dt (float) – Previous model time step (in s).

• _config (SimulationConfig) – The simulation configuration object.

• _bodies (Bodies) – The main SoA container.

• _profiling (_ProfilingData) – A container for performance profiling timers.

• _collision_graph (list[list[int]]) – Graph of collisions.

• _colours (list[int]) – Colours used for spatialised collision handling.

• _pool (Pool) – A pool of worker processes for parallel computation.

• _num_processes (int) – The number of worker processes in the pool.

• _last_total_energy (float) – The total energy from the previous logged step.

• _current_iteration (int) – Helper for detailed logging.

• _logged_collision_details (bool) – Helper for detailed logging.

property bodies: list[BodyProxy]

Provides an AoS-like list of body proxies for the display controller.

Returns:

List of body proxies.

Return type:

list[BodyProxy]

cleanup() → None

Clean up resources, such as closing the multiprocessing pool.

compute_dynamics(iteration: int) → None

Computes one step of the simulation, including dynamics and collisions.

This method follows the same order of operations as the C++ implementation to ensure consistent behavior and valid benchmark comparisons.

Parameters:

iteration (int) – Current iteration.

get_body_data_for_display() → tuple[ndarray, ndarray, ndarray, ndarray]

Provides current simulation data required for display as NumPy arrays. This is optimised for efficient data transfer to the UI thread.

Returns:

A tuple containing (positions, quaternions, torques, angular_accelerations).

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]

print_profiling_report() → None

Prints a formatted report of the time spent in different stages of the simulation if diagnostics are enabled.