🧱 Mason: Decouple Momentum Equations (Fixes OCP)

.jules/mason.md

@@ -30,3 +30,7 @@
 ## 2026-10-24 - Decoupling RL Exploration Logic
 **Violation:** Dependency Inversion Principle (DIP) and Open/Closed Principle (OCP). The `TabularQAgent` hardcoded the Epsilon-Greedy exploration strategy and `rand::thread_rng()`, making it non-deterministic, hard to test, and impossible to extend with new strategies (like Boltzmann or UCB) without modifying the agent.
 **Remedy:** Strategy Pattern. Extracted `ExplorationStrategy` trait and implemented `EpsilonGreedy`. Injected the strategy into `TabularQAgent` and exposed `select_action_with_rng` for deterministic testing.
+
+## 2026-10-25 - Decoupling Fluid Momentum Equations
+**Violation:** Open/Closed Principle (OCP). The `conservation.rs` module hardcoded `navier_stokes_time_derivative` and `euler_time_derivative` as separate functions, making it difficult to extend with new physics models (e.g., Stokes flow) without modifying the module and its consumers.
+**Remedy:** Strategy Pattern. Extracted `MomentumEquation` trait and implemented `NavierStokes` and `Euler` strategies. Introduced `SpatialGradients` parameter object to unify signatures. Refactored original functions as backward-compatible wrappers.

math_explorer/src/physics/fluid_dynamics/conservation.rs

@@ -2,8 +2,9 @@
 //!
 //! Implements the core Partial Differential Equations (PDEs) of Fluid Dynamics.
 
-use super::types::{FlowState, FluidProperties};
-use nalgebra::Vector3;
+use super::error::FluidError;
+use super::types::{FlowState, FluidProperties, SpatialGradients};
+use nalgebra::{Matrix3, Vector3};
 
 /// Calculates the Material Derivative ($D/Dt$) of a scalar property.
 ///
@@ -30,7 +31,7 @@ pub fn material_derivative_scalar(
 pub fn material_derivative_vector(
     local_change: Vector3<f64>,
     velocity: Vector3<f64>,
-    gradient_tensor: &nalgebra::Matrix3<f64>,
+    gradient_tensor: &Matrix3<f64>,
 ) -> Vector3<f64> {
     // (\mathbf{u} \cdot \nabla) \mathbf{A} corresponds to Jacobian * velocity vector
     // J = [ dA_x/dx  dA_x/dy  dA_x/dz ]
@@ -49,6 +50,80 @@ pub fn continuity_divergence(velocity_divergence: f64) -> f64 {
     velocity_divergence
 }
 
+// --- Strategy Pattern for Momentum Equation ---
+
+/// Defines a strategy for calculating the time evolution of velocity.
+pub trait MomentumEquation {
+    /// Computes $\frac{\partial \mathbf{u}}{\partial t}$ based on the conservation of momentum.
+    fn compute_time_derivative(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError>;
+}
+
+/// Navier-Stokes Strategy (Viscous Flow).
+#[derive(Debug, Clone, Copy, Default)]
+pub struct NavierStokes;
+
+impl MomentumEquation for NavierStokes {
+    fn compute_time_derivative(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError> {
+        let laplacian = gradients
+            .laplacian_velocity
+            .ok_or(FluidError::MissingLaplacian)?;
+
+        let nu = properties.kinematic_viscosity();
+        let convection = convective_acceleration(&state.velocity, &gradients.velocity_gradient);
+        let pressure = pressure_acceleration(&gradients.pressure_gradient, properties.density);
+
+        // Viscous term: nu * del^2 u
+        let viscous_term = laplacian * nu;
+
+        Ok(convection + pressure + viscous_term + body_force)
+    }
+}
+
+/// Euler Strategy (Inviscid Flow).
+#[derive(Debug, Clone, Copy, Default)]
+pub struct Euler;
+
+impl MomentumEquation for Euler {
+    fn compute_time_derivative(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError> {
+        let convection = convective_acceleration(&state.velocity, &gradients.velocity_gradient);
+        let pressure = pressure_acceleration(&gradients.pressure_gradient, properties.density);
+
+        Ok(convection + pressure + body_force)
+    }
+}
+
+// --- Helper Functions (DRY) ---
+
+fn convective_acceleration(velocity: &Vector3<f64>, gradient: &Matrix3<f64>) -> Vector3<f64> {
+    // -(u . del) u
+    -(gradient * velocity)
+}
+
+fn pressure_acceleration(pressure_gradient: &Vector3<f64>, density: f64) -> Vector3<f64> {
+    // -(1/rho) grad p
+    -pressure_gradient / density
+}
+
+// --- Legacy Functions (Wrappers) ---
+
 /// Computes the time evolution of velocity based on the Navier-Stokes Equations.
 ///
 /// $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{g}$$
@@ -63,47 +138,50 @@ pub fn continuity_divergence(velocity_divergence: f64) -> f64 {
 /// * `body_force`: External forces (e.g., gravity $\mathbf{g}$). Note: Input is acceleration vector (force per unit mass), or force vector if divided by rho manually.
 ///   Standard form usually takes body force density $\mathbf{f}$. If $\mathbf{f} = \rho \mathbf{g}$, then term is $\mathbf{g}$.
 ///   Here we assume `body_force` is $\mathbf{g}$ (acceleration).
+#[deprecated(note = "Use NavierStokes strategy instead")]
 pub fn navier_stokes_time_derivative(
     properties: &FluidProperties,
     state: &FlowState,
-    velocity_gradient: &nalgebra::Matrix3<f64>,
+    velocity_gradient: &Matrix3<f64>,
     pressure_gradient: Vector3<f64>,
     laplacian_velocity: Vector3<f64>,
     body_force_accel: Vector3<f64>,
 ) -> Vector3<f64> {
-    let nu = properties.kinematic_viscosity();
-    let rho = properties.density;
-
-    // Convective term: -(u . del) u
-    let convection = -(velocity_gradient * state.velocity);
-
-    // Pressure term: -(1/rho) grad p
-    let pressure_term = -pressure_gradient / rho;
-
-    // Viscous term: nu * del^2 u
-    let viscous_term = laplacian_velocity * nu;
+    let gradients = SpatialGradients {
+        velocity_gradient: *velocity_gradient,
+        pressure_gradient,
+        laplacian_velocity: Some(laplacian_velocity),
+    };
 
-    // Sum
-    convection + pressure_term + viscous_term + body_force_accel
+    NavierStokes
+        .compute_time_derivative(properties, state, &gradients, body_force_accel)
+        .expect("Laplacian guaranteed by function signature")
 }
 
 /// Computes the time evolution of velocity based on the Euler Equations (Inviscid).
 ///
 /// $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \mathbf{g}$$
 ///
 /// Assumes $\mu = 0$.
+#[deprecated(note = "Use Euler strategy instead")]
 pub fn euler_time_derivative(
     rho: f64,
     state: &FlowState,
-    velocity_gradient: &nalgebra::Matrix3<f64>,
+    velocity_gradient: &Matrix3<f64>,
     pressure_gradient: Vector3<f64>,
     body_force_accel: Vector3<f64>,
 ) -> Vector3<f64> {
-    // Convective term: -(u . del) u
-    let convection = -(velocity_gradient * state.velocity);
+    // Euler strategy requires FluidProperties, but ignores viscosity.
+    // We create a dummy properties with the given density.
+    let properties = FluidProperties::new(rho, 0.0);
 
-    // Pressure term: -(1/rho) grad p
-    let pressure_term = -pressure_gradient / rho;
+    let gradients = SpatialGradients {
+        velocity_gradient: *velocity_gradient,
+        pressure_gradient,
+        laplacian_velocity: None,
+    };
 
-    convection + pressure_term + body_force_accel
+    Euler
+        .compute_time_derivative(&properties, state, &gradients, body_force_accel)
+        .expect("Euler does not require laplacian")
 }

math_explorer/src/physics/fluid_dynamics/error.rs

@@ -0,0 +1,26 @@
+//! Errors for Fluid Dynamics.
+
+use std::fmt;
+
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub enum FluidError {
+    /// The Laplacian of velocity is required but was not provided.
+    MissingLaplacian,
+    /// Generic invalid configuration.
+    InvalidConfiguration(String),
+}
+
+impl fmt::Display for FluidError {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            FluidError::MissingLaplacian => {
+                write!(f, "Laplacian of velocity is required for this operation")
+            }
+            FluidError::InvalidConfiguration(msg) => {
+                write!(f, "Invalid fluid configuration: {}", msg)
+            }
+        }
+    }
+}
+
+impl std::error::Error for FluidError {}

math_explorer/src/physics/fluid_dynamics/mod.rs

@@ -1,5 +1,6 @@
 pub mod analysis;
 pub mod conservation;
+pub mod error;
 pub mod regimes;
 pub mod turbulence;
 pub mod types;

math_explorer/src/physics/fluid_dynamics/tests.rs

@@ -3,10 +3,12 @@ mod tests {
     use crate::physics::fluid_dynamics::{
         analysis::{bernoulli_constant, reynolds_number, shear_stress},
         conservation::{
-            continuity_divergence, material_derivative_scalar, navier_stokes_time_derivative,
+            Euler, MomentumEquation, NavierStokes, continuity_divergence,
+            material_derivative_scalar, navier_stokes_time_derivative,
         },
+        error::FluidError,
         regimes::{FlatPlateClassifier, FlowClassifier, FlowRegime, PipeFlowClassifier},
-        types::{FlowState, FluidProperties},
+        types::{FlowState, FluidProperties, SpatialGradients},
     };
     use nalgebra::{Matrix3, Vector3};
 
@@ -90,6 +92,7 @@ mod tests {
     }
 
     #[test]
+    #[allow(deprecated)] // Testing deprecated wrapper
     fn test_navier_stokes_simple_couette() {
         let props = FluidProperties::new(1.0, 1.0); // rho=1, mu=1 -> nu=1
         let state = FlowState::new(Vector3::new(1.0, 0.0, 0.0), 0.0);
@@ -113,4 +116,64 @@ mod tests {
         assert_eq!(continuity_divergence(0.0), 0.0);
         assert_eq!(continuity_divergence(0.5), 0.5);
     }
+
+    #[test]
+    fn test_navier_stokes_strategy() {
+        let props = FluidProperties::new(1.0, 1.0);
+        let state = FlowState::new(Vector3::new(1.0, 0.0, 0.0), 0.0);
+        let g = Vector3::zeros();
+
+        let gradients = SpatialGradients {
+            velocity_gradient: Matrix3::zeros(),
+            pressure_gradient: Vector3::new(2.0, 0.0, 0.0),
+            laplacian_velocity: Some(Vector3::zeros()),
+        };
+
+        let strategy = NavierStokes;
+        let accel = strategy
+            .compute_time_derivative(&props, &state, &gradients, g)
+            .unwrap();
+
+        // Pressure term: -1/rho * grad_p = -2.0
+        assert!((accel.x - (-2.0)).abs() < 1e-9);
+    }
+
+    #[test]
+    fn test_navier_stokes_missing_laplacian() {
+        let props = FluidProperties::new(1.0, 1.0);
+        let state = FlowState::new(Vector3::zeros(), 0.0);
+        let g = Vector3::zeros();
+
+        let gradients = SpatialGradients {
+            velocity_gradient: Matrix3::zeros(),
+            pressure_gradient: Vector3::zeros(),
+            laplacian_velocity: None, // Missing!
+        };
+
+        let strategy = NavierStokes;
+        let result = strategy.compute_time_derivative(&props, &state, &gradients, g);
+
+        assert_eq!(result, Err(FluidError::MissingLaplacian));
+    }
+
+    #[test]
+    fn test_euler_strategy() {
+        let props = FluidProperties::new(1.0, 0.0);
+        let state = FlowState::new(Vector3::new(1.0, 0.0, 0.0), 0.0);
+        let g = Vector3::zeros();
+
+        // Laplacian is None, but Euler shouldn't care
+        let gradients = SpatialGradients {
+            velocity_gradient: Matrix3::zeros(),
+            pressure_gradient: Vector3::new(2.0, 0.0, 0.0),
+            laplacian_velocity: None,
+        };
+
+        let strategy = Euler;
+        let accel = strategy
+            .compute_time_derivative(&props, &state, &gradients, g)
+            .unwrap();
+
+        assert!((accel.x - (-2.0)).abs() < 1e-9);
+    }
 }

math_explorer/src/physics/fluid_dynamics/types.rs

@@ -56,3 +56,33 @@ impl FlowState {
         Self { velocity, pressure }
     }
 }
+
+/// Encapsulates spatial gradients of flow properties.
+///
+/// This struct groups related derivatives to avoid primitive obsession and
+/// long argument lists in momentum equations.
+#[derive(Debug, Clone, Copy)]
+pub struct SpatialGradients {
+    /// Gradient of velocity ($\nabla \mathbf{u}$).
+    pub velocity_gradient: nalgebra::Matrix3<f64>,
+    /// Gradient of pressure ($\nabla p$).
+    pub pressure_gradient: Vector3<f64>,
+    /// Laplacian of velocity ($\nabla^2 \mathbf{u}$).
+    /// Optional because it is not required for Euler (inviscid) flow.
+    pub laplacian_velocity: Option<Vector3<f64>>,
+}
+
+impl SpatialGradients {
+    /// Creates a new `SpatialGradients` struct.
+    pub fn new(
+        velocity_gradient: nalgebra::Matrix3<f64>,
+        pressure_gradient: Vector3<f64>,
+        laplacian_velocity: Option<Vector3<f64>>,
+    ) -> Self {
+        Self {
+            velocity_gradient,
+            pressure_gradient,
+            laplacian_velocity,
+        }
+    }
+}