feat[physics] :: improve numerical stability and refactor simulation loop#9
Conversation
Add explicit Result-based error handling to physics APIs, implement numerical stability with input validation and safe math operations, split physics simulation into clear phases, and reorganize codebase with dedicated numerics/ and constraints/ modules. Key improvements: - Convert public APIs to return Result for proper error propagation - Add validation and clamping to prevent numerical instability - Separate integration, collision detection, and resolution phases - Create numerics/ module for constants and utility functions - Create constraints/ module for physics constraint system - Update pre-commit hooks to v6.0.0 and fix CI configuration - Expand test coverage for edge cases and error conditions - Update examples and tests for new Result-based API This makes the physics engine production-ready with robust error handling, numerical stability, and maintainable module boundaries.
Summary of ChangesHello @bniladridas, I'm Gemini Code Assist1! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed! This pull request significantly enhances the robustness and numerical stability of the physics engine by introducing dedicated modules for managing physics constraints and numerical utilities. It refactors the core simulation loop into distinct, well-defined phases, and implements comprehensive error handling and validation for critical operations. These changes are designed to prevent common simulation issues, leading to a more reliable and predictable physics system. Highlights
Ignored Files
Using Gemini Code AssistThe full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips. Invoking Gemini You can request assistance from Gemini at any point by creating a comment using either
Customization To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a Limitations & Feedback Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here. You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension. Footnotes
|
There was a problem hiding this comment.
Code Review
This pull request significantly enhances the physics engine's robustness by introducing input validation, numerical stability improvements, and better code structure. Key changes include adding a numerics module for constants and utilities, a constraints module, and refactoring the PhysicsEngine to use a phased simulation loop and return Result types for failable operations. The changes are very positive, but there are several areas for improvement, mainly around code duplication and the use of magic numbers. My review focuses on eliminating duplicated code by using the new numerics module, replacing hardcoded values with named constants, and fixing a logic issue in the new constraints implementation.
| // Temporary inline implementations until module issue is resolved | ||
| const MIN_DISTANCE_EPSILON: f32 = 1e-6; | ||
| const MAX_VELOCITY: f32 = 1000.0; | ||
|
|
||
| fn clamp<T: PartialOrd>(value: T, min: T, max: T) -> T { | ||
| if value < min { | ||
| min | ||
| } else if value > max { | ||
| max | ||
| } else { | ||
| value | ||
| } | ||
| } | ||
|
|
||
| fn safe_divide(numerator: f32, denominator: f32, default: f32) -> f32 { | ||
| if denominator.abs() < MIN_DISTANCE_EPSILON { | ||
| default | ||
| } else { | ||
| numerator / denominator | ||
| } | ||
| } | ||
|
|
||
| fn clamp_velocity(velocity: &mut Vector2<f32>) { | ||
| velocity.x = clamp(velocity.x, -MAX_VELOCITY, MAX_VELOCITY); | ||
| velocity.y = clamp(velocity.y, -MAX_VELOCITY, MAX_VELOCITY); | ||
| } | ||
|
|
||
| fn validate_vector(vec: &Vector2<f32>) -> Result<(), String> { | ||
| if !vec.x.is_finite() || !vec.y.is_finite() { | ||
| Err(format!("Vector contains non-finite values: {:?}", vec)) | ||
| } else { | ||
| Ok(()) | ||
| } | ||
| } | ||
|
|
||
| fn approx_eq(a: f32, b: f32, epsilon: f32) -> bool { | ||
| (a - b).abs() < epsilon | ||
| } |
There was a problem hiding this comment.
These helper functions and constants are duplicates of the ones defined in the new numerics module (engine/numerics/utilities.rs and engine/numerics/constants.rs). This code duplication should be removed to ensure a single source of truth. The functions and constants from the numerics module should be imported and used throughout this file instead.
| if position.x <= min.x || position.x >= max.x { | ||
| velocity.x *= 1.0 - self.damping; | ||
| } | ||
| if position.y <= min.y || position.y >= max.y { | ||
| velocity.y *= 1.0 - self.damping; | ||
| } |
There was a problem hiding this comment.
The current logic for applying damping at boundaries is incorrect. It dampens the velocity regardless of the object's direction of movement. This means an object moving away from a boundary will still be slowed down. The damping should only be applied if the object is moving towards the boundary.
Also, note that this implementation only dampens the velocity; it doesn't cause a bounce (velocity reversal). This is inconsistent with the wall collision handling in physics_engine.rs, which both reverses and dampens velocity. If a bounce is intended, the velocity component should be negated.
| if position.x <= min.x || position.x >= max.x { | |
| velocity.x *= 1.0 - self.damping; | |
| } | |
| if position.y <= min.y || position.y >= max.y { | |
| velocity.y *= 1.0 - self.damping; | |
| } | |
| if (position.x <= min.x && velocity.x < 0.0) || (position.x >= max.x && velocity.x > 0.0) { | |
| velocity.x *= 1.0 - self.damping; | |
| } | |
| if (position.y <= min.y && velocity.y < 0.0) || (position.y >= max.y && velocity.y > 0.0) { | |
| velocity.y *= 1.0 - self.damping; | |
| } |
| pub fn new(position: Vector2<f32>, velocity: Vector2<f32>, mass: f32) -> Result<Self, String> { | ||
| // Temporary inline validation | ||
| if mass < 1e-3 { | ||
| return Err(format!("Mass {} is below minimum threshold {}", mass, 1e-3)); | ||
| } | ||
| if mass > 1e6 { | ||
| return Err(format!("Mass {} exceeds maximum threshold {}", mass, 1e6)); | ||
| } | ||
| if !mass.is_finite() { | ||
| return Err("Mass must be finite".to_string()); | ||
| } | ||
| if !position.x.is_finite() || !position.y.is_finite() { | ||
| return Err(format!( | ||
| "Position contains non-finite values: {:?}", | ||
| position | ||
| )); | ||
| } | ||
| if !velocity.x.is_finite() || !velocity.y.is_finite() { | ||
| return Err(format!( | ||
| "Velocity contains non-finite values: {:?}", | ||
| velocity | ||
| )); | ||
| } | ||
|
|
||
| Ok(Object { | ||
| position, | ||
| velocity, | ||
| mass, | ||
| radius: 0.5, | ||
| } | ||
| }) | ||
| } |
There was a problem hiding this comment.
The validation logic within Object::new is duplicated. It should be refactored to use the validate_mass and validate_vector functions from the numerics::utilities module to avoid code duplication.
Additionally, it uses hardcoded magic numbers for mass validation (e.g., 1e-3 on line 53, 1e6 on line 56). These should be replaced with the MIN_MASS and MAX_MASS constants from the numerics::constants module.
| _ => { | ||
| // Other constraint types not yet implemented | ||
| } |
There was a problem hiding this comment.
The apply method currently does nothing for Distance and Axis constraints. This can lead to silent failures where constraints are expected to be applied but are not. It's better to use unimplemented!() to make it explicit that these are not yet supported. This will cause a panic during development if these constraint types are used, which is desirable for unimplemented features.
| _ => { | |
| // Other constraint types not yet implemented | |
| } | |
| _ => { | |
| unimplemented!("This constraint type is not yet implemented"); | |
| } |
| use crate::numerics::constants::*; | ||
| /// Numerical utilities and helper functions | ||
| /// Provides safe mathematical operations and validation functions | ||
| use nalgebra::Vector2; |
There was a problem hiding this comment.
Module-level documentation should use //! instead of ///. /// documents the following item (in this case, a use statement), whereas //! documents the enclosing module. These doc comments should also be at the very top of the file, before any use statements.
| use crate::numerics::constants::*; | |
| /// Numerical utilities and helper functions | |
| /// Provides safe mathematical operations and validation functions | |
| use nalgebra::Vector2; | |
| //! Numerical utilities and helper functions | |
| //! Provides safe mathematical operations and validation functions | |
| use crate::numerics::constants::*; | |
| use nalgebra::Vector2; |
|
|
||
| Ok(PhysicsEngine { | ||
| objects: vec![default_object], // default object | ||
| gravity: Vector2::new(0.0, -9.81), // downward acceleration |
There was a problem hiding this comment.
The gravity value is hardcoded here. The numerics::constants module now provides an EARTH_GRAVITY constant. Using the constant here will improve consistency and maintainability.
| gravity: Vector2::new(0.0, -9.81), // downward acceleration | |
| gravity: Vector2::new(0.0, -crate::numerics::EARTH_GRAVITY), // downward acceleration |
| // Elastic collision impulse calculation | ||
| let m1 = self.objects[i].mass; | ||
| let m2 = self.objects[j].mass; | ||
| let restitution = 1.0; // perfectly elastic |
There was a problem hiding this comment.
| let m1 = self.objects[i].mass; | ||
| let m2 = self.objects[j].mass; | ||
| let restitution = 1.0; // perfectly elastic | ||
| let impulse_scalar = -(1.0 + restitution) * velocity_along_normal * (m1 * m2) / (m1 + m2); |
There was a problem hiding this comment.
This formula for impulse_scalar is mathematically equivalent to the previous one (-(1.0 + restitution) * velocity_along_normal / (1.0 / m1 + 1.0 / m2)), but the previous form is generally preferred in physics engines. It handles cases with infinite mass (e.g., for static objects where inverse mass is 0) more gracefully without requiring special checks. While not an issue with the current implementation due to mass validation, using the inverse mass formulation is a more robust and standard practice.
| let impulse_scalar = -(1.0 + restitution) * velocity_along_normal * (m1 * m2) / (m1 + m2); | |
| let impulse_scalar = -(1.0 + restitution) * velocity_along_normal / (1.0 / m1 + 1.0 / m2); |
| if approx_eq(obj.position.x, self.bounds.0.x, MIN_DISTANCE_EPSILON) | ||
| || approx_eq(obj.position.x, self.bounds.1.x, MIN_DISTANCE_EPSILON) | ||
| { | ||
| obj.velocity.x = -obj.velocity.x * 0.8; |
| obj.position.x = clamp(obj.position.x, -1e6, 1e6); | ||
| obj.position.y = clamp(obj.position.y, -1e6, 1e6); |
Summary
This PR improves the robustness and numerical stability of the physics engine by restructuring core simulation logic, introducing dedicated numerical utilities, and adding explicit error handling throughout the API.
Key Changes
Numerical stability & error handling
PhysicsEngine::new,add_object,simulate) now returnResultto prevent invalid engine states.Modularized physics logic
constraintsmodule for defining object constraints (e.g. fixed position, boundaries).numericsmodule to centralize constants and math helpers.Refactored simulation loop
Tooling
pre-commit-hookstov6.0.0.