⚙️ Core: Momentum Equation Strategy Pattern

.jules/systems_core.md

@@ -58,3 +58,8 @@
 **Violation:** The `gillespie_step_time` function in `epidemiology::stochastic` had hardcoded `rand::thread_rng()` dependency (Side Effect) and hardcoded 2-reaction logic, violating Dependency Inversion and Open/Closed Principle.
 **Refactor:** Applied the **Strategy Pattern**. Defined `StochasticSystem` trait for reaction networks and `GillespieSolver` struct with injected RNG (`R: Rng`).
 **Trade-off:** Increased complexity (Generics, Trait definition) compared to a single function, but enabled deterministic testing (seeded RNG) and support for any reaction network (Composability).
+
+## 2027-02-15 - Strategy Pattern for Momentum Equation
+**Violation:** The fluid dynamics conservation module used separate functions (`navier_stokes_time_derivative`, `euler_time_derivative`) with repetitive logic and "Tuple Soup" signatures, violating SRP, DRY, and OCP.
+**Refactor:** Applied **Strategy Pattern**. Extracted `MomentumEquation` trait and `SpatialGradients` parameter object. Implemented `NavierStokes` and `Euler` strategies.
+**Trade-off:** Introduced `FluidError` handling (breaking change wrapped in `Result`), but enabled unified solver interfaces and extensibility for new flow regimes (e.g., Stokes flow).

math_explorer/src/physics/fluid_dynamics/conservation.rs

@@ -2,7 +2,8 @@
 //!
 //! Implements the core Partial Differential Equations (PDEs) of Fluid Dynamics.
 
-use super::types::{FlowState, FluidProperties};
+use super::error::FluidError;
+use super::types::{FlowState, FluidProperties, SpatialGradients};
 use nalgebra::Vector3;
 
 /// Calculates the Material Derivative ($D/Dt$) of a scalar property.
@@ -49,61 +50,112 @@ pub fn continuity_divergence(velocity_divergence: f64) -> f64 {
     velocity_divergence
 }
 
-/// Computes the time evolution of velocity based on the Navier-Stokes Equations.
+/// Strategy for calculating the time evolution of velocity (acceleration).
+pub trait MomentumEquation {
+    /// Calculates $\frac{\partial \mathbf{u}}{\partial t}$.
+    fn calculate_acceleration(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError>;
+}
+
+/// Navier-Stokes equations for viscous flow.
 ///
 /// $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{g}$$
+#[derive(Debug, Clone, Copy, Default)]
+pub struct NavierStokes;
+
+impl MomentumEquation for NavierStokes {
+    fn calculate_acceleration(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError> {
+        let nu = properties.kinematic_viscosity();
+        let rho = properties.density;
+
+        let laplacian = gradients
+            .laplacian_velocity
+            .ok_or(FluidError::MissingLaplacian)?;
+
+        // Convective term: -(u . del) u
+        let convection = -(gradients.velocity_gradient * state.velocity);
+
+        // Pressure term: -(1/rho) grad p
+        let pressure_term = -gradients.pressure_gradient / rho;
+
+        // Viscous term: nu * del^2 u
+        let viscous_term = laplacian * nu;
+
+        // Sum
+        Ok(convection + pressure_term + viscous_term + body_force_accel)
+    }
+}
+
+/// Euler equations for inviscid flow.
+///
+/// $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \mathbf{g}$$
 ///
-/// Returns the local acceleration $\frac{\partial \mathbf{u}}{\partial t}$.
+/// Assumes $\mu = 0$.
+#[derive(Debug, Clone, Copy, Default)]
+pub struct Euler;
+
+impl MomentumEquation for Euler {
+    fn calculate_acceleration(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError> {
+        let rho = properties.density;
+
+        // Convective term: -(u . del) u
+        let convection = -(gradients.velocity_gradient * state.velocity);
+
+        // Pressure term: -(1/rho) grad p
+        let pressure_term = -gradients.pressure_gradient / rho;
+
+        Ok(convection + pressure_term + body_force_accel)
+    }
+}
+
+/// Computes the time evolution of velocity based on the Navier-Stokes Equations.
 ///
-/// * `properties`: Fluid properties ($\rho, \mu$).
-/// * `state`: Current flow state ($\mathbf{u}, p$).
-/// * `velocity_gradient`: Jacobian of velocity ($\nabla \mathbf{u}$).
-/// * `pressure_gradient`: Gradient of pressure ($\nabla p$).
-/// * `laplacian_velocity`: Laplacian of velocity ($\nabla^2 \mathbf{u}$).
-/// * `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).
+/// **Refactored:** Now a wrapper around the `NavierStokes` strategy.
 pub fn navier_stokes_time_derivative(
     properties: &FluidProperties,
     state: &FlowState,
     velocity_gradient: &nalgebra::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;
-
-    // Sum
-    convection + pressure_term + viscous_term + body_force_accel
+) -> Result<Vector3<f64>, FluidError> {
+    let gradients = SpatialGradients::new(
+        *velocity_gradient,
+        pressure_gradient,
+        Some(laplacian_velocity),
+    );
+    NavierStokes.calculate_acceleration(properties, state, &gradients, body_force_accel)
 }
 
 /// 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$.
+/// **Refactored:** Now a wrapper around the `Euler` strategy.
 pub fn euler_time_derivative(
     rho: f64,
     state: &FlowState,
     velocity_gradient: &nalgebra::Matrix3<f64>,
     pressure_gradient: Vector3<f64>,
     body_force_accel: Vector3<f64>,
-) -> Vector3<f64> {
-    // Convective term: -(u . del) u
-    let convection = -(velocity_gradient * state.velocity);
-
-    // Pressure term: -(1/rho) grad p
-    let pressure_term = -pressure_gradient / rho;
-
-    convection + pressure_term + body_force_accel
+) -> Result<Vector3<f64>, FluidError> {
+    let gradients = SpatialGradients::new(*velocity_gradient, pressure_gradient, None);
+    // Create temporary properties with correct density and 0 viscosity.
+    let props = FluidProperties::new(rho, 0.0);
+    Euler.calculate_acceleration(&props, state, &gradients, body_force_accel)
 }

math_explorer/src/physics/fluid_dynamics/error.rs

@@ -0,0 +1,26 @@
+//! Error types for Fluid Dynamics.
+
+use std::error::Error;
+use std::fmt;
+
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub enum FluidError {
+    /// The Laplacian of velocity is required but was not provided.
+    MissingLaplacian,
+    /// Invalid configuration for the flow model.
+    InvalidConfiguration(String),
+    /// Other errors.
+    Other(String),
+}
+
+impl fmt::Display for FluidError {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            FluidError::MissingLaplacian => write!(f, "Missing Laplacian of velocity"),
+            FluidError::InvalidConfiguration(msg) => write!(f, "Invalid configuration: {}", msg),
+            FluidError::Other(msg) => write!(f, "Fluid dynamics error: {}", msg),
+        }
+    }
+}
+
+impl 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, euler_time_derivative,
+            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};
 
@@ -99,12 +101,14 @@ mod tests {
         let lap_vel = Vector3::zeros();
         let g = Vector3::zeros();
 
-        let accel = navier_stokes_time_derivative(&props, &state, &vel_grad, p_grad, lap_vel, g);
+        let accel = navier_stokes_time_derivative(&props, &state, &vel_grad, p_grad, lap_vel, g)
+            .expect("Valid NS calculation");
         assert_eq!(accel, Vector3::zeros());
 
         let p_grad_x = Vector3::new(2.0, 0.0, 0.0);
         let accel_p =
-            navier_stokes_time_derivative(&props, &state, &vel_grad, p_grad_x, lap_vel, g);
+            navier_stokes_time_derivative(&props, &state, &vel_grad, p_grad_x, lap_vel, g)
+                .expect("Valid NS calculation");
         assert!((accel_p.x - (-2.0)).abs() < 1e-9);
     }
 
@@ -113,4 +117,30 @@ mod tests {
         assert_eq!(continuity_divergence(0.0), 0.0);
         assert_eq!(continuity_divergence(0.5), 0.5);
     }
+
+    #[test]
+    fn test_euler_equation() {
+        // Inviscid flow, so viscosity shouldn't matter (Euler ignores it)
+        let rho = 2.0;
+        let state = FlowState::new(Vector3::zeros(), 0.0);
+        let vel_grad = Matrix3::zeros();
+        let p_grad = Vector3::new(4.0, 0.0, 0.0);
+        let g = Vector3::zeros();
+
+        // dp/dx = 4, rho = 2 -> accel = -4/2 = -2
+        let accel = euler_time_derivative(rho, &state, &vel_grad, p_grad, g)
+            .expect("Valid Euler calculation");
+        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 gradients = SpatialGradients::new(Matrix3::zeros(), Vector3::zeros(), None);
+        let g = Vector3::zeros();
+
+        let result = NavierStokes.calculate_acceleration(&props, &state, &gradients, g);
+        assert!(matches!(result, Err(FluidError::MissingLaplacian)));
+    }
 }

math_explorer/src/physics/fluid_dynamics/types.rs

@@ -1,6 +1,6 @@
 //! Types for Fluid Dynamics.
 
-use nalgebra::Vector3;
+use nalgebra::{Matrix3, Vector3};
 
 /// Physical properties of the fluid.
 #[derive(Debug, Clone, Copy)]
@@ -56,3 +56,29 @@ impl FlowState {
         Self { velocity, pressure }
     }
 }
+
+/// Encapsulates spatial derivatives required for momentum equations.
+#[derive(Debug, Clone, Copy)]
+pub struct SpatialGradients {
+    /// Jacobian of velocity ($\nabla \mathbf{u}$).
+    pub velocity_gradient: Matrix3<f64>,
+    /// Gradient of pressure ($\nabla p$).
+    pub pressure_gradient: Vector3<f64>,
+    /// Laplacian of velocity ($\nabla^2 \mathbf{u}$). Optional, as Euler equations don't use it.
+    pub laplacian_velocity: Option<Vector3<f64>>,
+}
+
+impl SpatialGradients {
+    /// Creates a new `SpatialGradients` structure.
+    pub fn new(
+        velocity_gradient: Matrix3<f64>,
+        pressure_gradient: Vector3<f64>,
+        laplacian_velocity: Option<Vector3<f64>>,
+    ) -> Self {
+        Self {
+            velocity_gradient,
+            pressure_gradient,
+            laplacian_velocity,
+        }
+    }
+}