Merge pull request #581 from dimforge/solver-pick

feat: rework solver parameters to make it easy to recover the old behaviors

src/dynamics/integration_parameters.rs

@@ -44,17 +44,65 @@ pub struct IntegrationParameters {
     pub max_penetration_correction: Real,
     /// The maximal distance separating two objects that will generate predictive contacts (default: `0.002`).
     pub prediction_distance: Real,
-    /// Number of iterations performed to solve friction constraints at solver iteration (default: `2`).
-    pub num_friction_iteration_per_solver_iteration: usize,
     /// The number of solver iterations run by the constraints solver for calculating forces (default: `4`).
     pub num_solver_iterations: NonZeroUsize,
+    /// Number of addition friction resolution iteration run during the last solver sub-step (default: `4`).
+    pub num_additional_friction_iterations: usize,
+    /// Number of internal Project Gauss Seidel (PGS) iterations run at each solver iteration (default: `1`).
+    pub num_internal_pgs_iterations: usize,
     /// Minimum number of dynamic bodies in each active island (default: `128`).
     pub min_island_size: usize,
     /// Maximum number of substeps performed by the  solver (default: `1`).
     pub max_ccd_substeps: usize,
 }
 
 impl IntegrationParameters {
+    /// Configures the integration parameters to match the old PGS solver
+    /// from Rapier version <= 0.17.
+    ///
+    /// This solver was slightly faster than the new one but resulted
+    /// in less stable joints and worse convergence rates.
+    ///
+    /// This should only be used for comparison purpose or if you are
+    /// experiencing problems with the new solver.
+    ///
+    /// NOTE: this does not affect any [`RigidBody::additional_solver_iterations`] that will
+    ///       still create solver iterations based on the new "small-steps" PGS solver.
+    /// NOTE: this resets [`Self::erp`], [`Self::damping_ratio`], [`Self::joint_erp`],
+    ///       [`Self::joint_damping_ratio`] to their former default values.
+    pub fn switch_to_standard_pgs_solver(&mut self) {
+        self.num_internal_pgs_iterations =
+            self.num_solver_iterations.get() * self.num_internal_pgs_iterations;
+        self.num_solver_iterations = NonZeroUsize::new(1).unwrap();
+        self.erp = 0.8;
+        self.damping_ratio = 0.25;
+        self.joint_erp = 1.0;
+        self.joint_damping_ratio = 1.0;
+    }
+
+    /// Configures the integration parameters to match the new "small-steps" PGS solver
+    /// from Rapier version >= 0.18.
+    ///
+    /// The "small-steps" PGS solver is the default one given by [`Self::default()`] so
+    /// calling this function is generally not needed unless
+    /// [`Self::switch_to_standard_pgs_solver()`] was called.
+    ///
+    /// This solver results in more stable joints and significantly better convergence
+    /// rates but is slightly slower in its default settings.
+    ///
+    /// NOTE: this resets [`Self::erp`], [`Self::damping_ratio`], [`Self::joint_erp`],
+    ///       [`Self::joint_damping_ratio`] to their default values.
+    pub fn switch_to_small_steps_pgs_solver(&mut self) {
+        self.num_solver_iterations = NonZeroUsize::new(self.num_internal_pgs_iterations).unwrap();
+        self.num_internal_pgs_iterations = 1;
+
+        let default = Self::default();
+        self.erp = default.erp;
+        self.damping_ratio = default.damping_ratio;
+        self.joint_erp = default.joint_erp;
+        self.joint_damping_ratio = default.joint_damping_ratio;
+    }
+
     /// The inverse of the time-stepping length, i.e. the steps per seconds (Hz).
     ///
     /// This is zero if `self.dt` is zero.
@@ -154,7 +202,8 @@ impl Default for IntegrationParameters {
             allowed_linear_error: 0.001,
             max_penetration_correction: Real::MAX,
             prediction_distance: 0.002,
-            num_friction_iteration_per_solver_iteration: 2,
+            num_internal_pgs_iterations: 1,
+            num_additional_friction_iterations: 4,
             num_solver_iterations: NonZeroUsize::new(4).unwrap(),
             // TODO: what is the optimal value for min_island_size?
             // It should not be too big so that we don't end up with

src/dynamics/joint/rope_joint.rs

@@ -118,8 +118,11 @@ impl RopeJoint {
 
     /// The maximum distance allowed between the attached objects.
     #[must_use]
-    pub fn max_distance(&self) -> Option<Real> {
-        self.data.limits(JointAxis::X).map(|l| l.max)
+    pub fn max_distance(&self) -> Real {
+        self.data
+            .limits(JointAxis::X)
+            .map(|l| l.max)
+            .unwrap_or(Real::MAX)
     }
 
     /// Sets the maximum allowed distance between the attached objects.
@@ -146,8 +149,6 @@ pub struct RopeJointBuilder(pub RopeJoint);
 
 impl RopeJointBuilder {
     /// Creates a new builder for rope joints.
-    ///
-    /// This axis is expressed in the local-space of both rigid-bodies.
     pub fn new(max_dist: Real) -> Self {
         Self(RopeJoint::new(max_dist))
     }

src/dynamics/joint/spring_joint.rs

@@ -7,7 +7,7 @@ use crate::math::{Point, Real};
 #[repr(transparent)]
 /// A spring-damper joint, applies a force proportional to the distance between two objects.
 ///
-/// The spring is integrated implicitly, implying that an even undamped spring will still be subject to some
+/// The spring is integrated implicitly, implying that even an undamped spring will be subject to some
 /// amount of numerical damping (so it will eventually come to a rest). More solver iterations, or smaller
 /// timesteps, will lower the effect of numerical damping, providing a more realistic result.
 pub struct SpringJoint {

src/dynamics/rigid_body.rs

@@ -89,7 +89,7 @@ impl RigidBody {
     /// and every rigid-body interacting directly or indirectly with it (through joints
     /// or contacts). This implies a performance hit.
     ///
-    /// The default value is 0, meaning [`IntegrationParameters::num_solver_iterations`] will
+    /// The default value is 0, meaning exactly [`IntegrationParameters::num_solver_iterations`] will
     /// be used as number of solver iterations for this body.
     pub fn set_additional_solver_iterations(&mut self, additional_iterations: usize) {
         self.additional_solver_iterations = additional_iterations;

src/dynamics/solver/parallel_velocity_solver.rs

@@ -105,7 +105,7 @@ impl ParallelVelocitySolver {
          */
         {
             for i in 0..params.num_velocity_iterations_per_small_step {
-                let solve_friction = params.num_friction_iteration_per_solver_iteration + i
+                let solve_friction = params.num_additional_friction_iterations + i
                     >= params.num_velocity_iterations_per_small_step;
                 // Solve joints.
                 solve!(
@@ -156,11 +156,10 @@ impl ParallelVelocitySolver {
             }
 
             // Solve the remaining friction iterations.
-            let remaining_friction_iterations = if params
-                .num_friction_iteration_per_solver_iteration
+            let remaining_friction_iterations = if params.num_additional_friction_iterations
                 > params.num_velocity_iterations_per_small_step
             {
-                params.num_friction_iteration_per_solver_iteration
+                params.num_additional_friction_iterations
                     - params.num_velocity_iterations_per_small_step
             } else {
                 0

src/dynamics/solver/velocity_solver.rs

@@ -152,14 +152,14 @@ impl VelocitySolver {
     pub fn solve_constraints(
         &mut self,
         params: &IntegrationParameters,
-        num_solver_iterations: usize,
+        num_substeps: usize,
         bodies: &mut RigidBodySet,
         multibodies: &mut MultibodyJointSet,
         contact_constraints: &mut SolverConstraintsSet<ContactConstraintTypes>,
         joint_constraints: &mut SolverConstraintsSet<JointConstraintTypes>,
     ) {
-        for small_step_id in 0..num_solver_iterations {
-            let is_last_small_step = small_step_id == num_solver_iterations - 1;
+        for substep_id in 0..num_substeps {
+            let is_last_substep = substep_id == num_substeps - 1;
 
             for (solver_vels, incr) in self
                 .solver_vels
@@ -176,28 +176,27 @@ impl VelocitySolver {
              * Update & solve constraints.
              */
             joint_constraints.update(&params, multibodies, &self.solver_bodies);
-            contact_constraints.update(&params, small_step_id, multibodies, &self.solver_bodies);
+            contact_constraints.update(&params, substep_id, multibodies, &self.solver_bodies);
 
-            joint_constraints.solve(&mut self.solver_vels, &mut self.generic_solver_vels);
-            contact_constraints
-                .solve_restitution(&mut self.solver_vels, &mut self.generic_solver_vels);
-
-            let num_friction_iterations = if is_last_small_step {
-                params.num_friction_iteration_per_solver_iteration * num_solver_iterations
-                    - (num_solver_iterations - 1)
-            } else {
-                1
-            };
-
-            for _ in 0..num_friction_iterations {
+            for _ in 0..params.num_internal_pgs_iterations {
+                joint_constraints.solve(&mut self.solver_vels, &mut self.generic_solver_vels);
+                contact_constraints
+                    .solve_restitution(&mut self.solver_vels, &mut self.generic_solver_vels);
                 contact_constraints
                     .solve_friction(&mut self.solver_vels, &mut self.generic_solver_vels);
             }
 
+            if is_last_substep {
+                for _ in 0..params.num_additional_friction_iterations {
+                    contact_constraints
+                        .solve_friction(&mut self.solver_vels, &mut self.generic_solver_vels);
+                }
+            }
+
             /*
              * Integrate positions.
              */
-            self.integrate_positions(&params, is_last_small_step, bodies, multibodies);
+            self.integrate_positions(&params, is_last_substep, bodies, multibodies);
 
             /*
              * Resolution without bias.
@@ -211,7 +210,7 @@ impl VelocitySolver {
     pub fn integrate_positions(
         &mut self,
         params: &IntegrationParameters,
-        is_last_small_step: bool,
+        is_last_substep: bool,
         bodies: &mut RigidBodySet,
         multibodies: &mut MultibodyJointSet,
     ) {
@@ -241,9 +240,9 @@ impl VelocitySolver {
             multibody.velocities.copy_from(&solver_vels);
             multibody.integrate(params.dt);
             // PERF: we could have a mode where it doesn’t write back to the `bodies` yet.
-            multibody.forward_kinematics(bodies, !is_last_small_step);
+            multibody.forward_kinematics(bodies, !is_last_substep);
 
-            if !is_last_small_step {
+            if !is_last_substep {
                 // These are very expensive and not needed if we don’t
                 // have to run another step.
                 multibody.update_dynamics(params.dt, bodies);
@@ -260,7 +259,7 @@ impl VelocitySolver {
     pub fn writeback_bodies(
         &mut self,
         params: &IntegrationParameters,
-        num_solver_iterations: usize,
+        num_substeps: usize,
         islands: &IslandManager,
         island_id: usize,
         bodies: &mut RigidBodySet,
@@ -302,9 +301,9 @@ impl VelocitySolver {
                 //       solver integrating at intermediate sub-steps. Should we just switch
                 //       to interpolation?
                 rb.integrated_vels.linvel =
-                    solver_body.integrated_vels.linvel / num_solver_iterations as Real;
+                    solver_body.integrated_vels.linvel / num_substeps as Real;
                 rb.integrated_vels.angvel =
-                    solver_body.integrated_vels.angvel / num_solver_iterations as Real;
+                    solver_body.integrated_vels.angvel / num_substeps as Real;
                 rb.vels = new_vels;
                 rb.pos.next_position = solver_body.position;
             }

src_testbed/testbed.rs

@@ -101,6 +101,12 @@ bitflags! {
     }
 }
 
+#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+pub enum RapierSolverType {
+    SmallStepsPgs,
+    StandardPgs,
+}
+
 #[derive(Resource)]
 pub struct TestbedState {
     pub running: RunMode,
@@ -121,6 +127,7 @@ pub struct TestbedState {
     pub example_names: Vec<&'static str>,
     pub selected_example: usize,
     pub selected_backend: usize,
+    pub solver_type: RapierSolverType,
     pub physx_use_two_friction_directions: bool,
     pub snapshot: Option<PhysicsSnapshot>,
     nsteps: usize,
@@ -204,6 +211,7 @@ impl TestbedApp {
             example_names: Vec::new(),
             selected_example: 0,
             selected_backend: RAPIER_BACKEND,
+            solver_type: RapierSolverType::SmallStepsPgs,
             physx_use_two_friction_directions: true,
             nsteps: 1,
             camera_locked: false,
@@ -1189,6 +1197,14 @@ fn update_testbed(
 
             if state.selected_example != prev_example {
                 harness.physics.integration_parameters = IntegrationParameters::default();
+
+                match state.solver_type {
+                    RapierSolverType::SmallStepsPgs => {} // It’s already the default.
+                    RapierSolverType::StandardPgs => harness
+                        .physics
+                        .integration_parameters
+                        .switch_to_standard_pgs_solver(),
+                }
             }
 
             let selected_example = state.selected_example;

src_testbed/ui.rs

@@ -5,8 +5,8 @@ use std::num::NonZeroUsize;
 use crate::debug_render::DebugRenderPipelineResource;
 use crate::harness::Harness;
 use crate::testbed::{
-    RunMode, TestbedActionFlags, TestbedState, TestbedStateFlags, PHYSX_BACKEND_PATCH_FRICTION,
-    PHYSX_BACKEND_TWO_FRICTION_DIR,
+    RapierSolverType, RunMode, TestbedActionFlags, TestbedState, TestbedStateFlags,
+    PHYSX_BACKEND_PATCH_FRICTION, PHYSX_BACKEND_TWO_FRICTION_DIR,
 };
 
 use crate::PhysicsState;
@@ -107,17 +107,52 @@ pub fn update_ui(
             integration_parameters.num_solver_iterations =
                 NonZeroUsize::new(num_iterations).unwrap();
         } else {
+            let mut changed = false;
+            egui::ComboBox::from_label("solver type")
+                .width(150.0)
+                .selected_text(format!("{:?}", state.solver_type))
+                .show_ui(ui, |ui| {
+                    let solver_types = [
+                        RapierSolverType::SmallStepsPgs,
+                        RapierSolverType::StandardPgs,
+                    ];
+                    for sty in solver_types {
+                        changed = ui
+                            .selectable_value(&mut state.solver_type, sty, format!("{sty:?}"))
+                            .changed()
+                            || changed;
+                    }
+                });
+
+            if changed {
+                match state.solver_type {
+                    RapierSolverType::SmallStepsPgs => {
+                        integration_parameters.switch_to_small_steps_pgs_solver()
+                    }
+                    RapierSolverType::StandardPgs => {
+                        integration_parameters.switch_to_standard_pgs_solver()
+                    }
+                }
+            }
+
             let mut num_iterations = integration_parameters.num_solver_iterations.get();
             ui.add(Slider::new(&mut num_iterations, 1..=40).text("num solver iters."));
             integration_parameters.num_solver_iterations =
                 NonZeroUsize::new(num_iterations).unwrap();
 
             ui.add(
                 Slider::new(
-                    &mut integration_parameters.num_friction_iteration_per_solver_iteration,
+                    &mut integration_parameters.num_internal_pgs_iterations,
+                    1..=40,
+                )
+                .text("num internal PGS iters."),
+            );
+            ui.add(
+                Slider::new(
+                    &mut integration_parameters.num_additional_friction_iterations,
                     1..=40,
                 )
-                .text("frict. iters. per solver iters."),
+                .text("num additional frict. iters."),
             );
         }