⚙️ Core: Strategy Pattern for Fluid Conservation

.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-04-10 - Strategy Pattern for Fluid Momentum
+**Violation:** `conservation.rs` contained duplicate logic for Navier-Stokes and Euler time derivatives, with hardcoded physics regimes (viscous vs inviscid) and parameter soup (`velocity_gradient`, `pressure_gradient`, `laplacian`), violating OCP and DRY.
+**Refactor:** Applied the **Strategy Pattern**. Extracted `MomentumEquation` trait. Implemented `NavierStokes` and `Euler` strategies. Introduced `SpatialGradients` struct to encapsulate spatial derivatives.
+**Trade-off:** Minimal boilerplate for trait definitions, but achieved full separation of physics models and eliminated argument soup.

math_explorer/src/physics/fluid_dynamics/conservation.rs

@@ -2,9 +2,81 @@
 //!
 //! Implements the core Partial Differential Equations (PDEs) of Fluid Dynamics.
 
-use super::types::{FlowState, FluidProperties};
+use super::types::{FlowState, FluidProperties, SpatialGradients};
 use nalgebra::Vector3;
 
+/// Defines the strategy for computing the momentum conservation term.
+pub trait MomentumEquation {
+    /// Computes the local acceleration $\frac{\partial \mathbf{u}}{\partial t}$.
+    fn compute_time_derivative(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Vector3<f64>;
+}
+
+/// Navier-Stokes Momentum Equation (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 compute_time_derivative(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Vector3<f64> {
+        let nu = properties.kinematic_viscosity();
+        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;
+
+        // Viscous term: nu * del^2 u
+        let viscous_term = match gradients.laplacian_velocity {
+            Some(laplacian) => laplacian * nu,
+            None => Vector3::zeros(),
+        };
+
+        // Sum
+        convection + pressure_term + viscous_term + body_force_accel
+    }
+}
+
+/// Euler Momentum Equation (Inviscid Flow).
+///
+/// $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \mathbf{g}$$
+#[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_accel: Vector3<f64>,
+    ) -> Vector3<f64> {
+        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;
+
+        convection + pressure_term + body_force_accel
+    }
+}
+
 /// Calculates the Material Derivative ($D/Dt$) of a scalar property.
 ///
 /// $$\frac{D\phi}{Dt} = \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi$$
@@ -71,20 +143,12 @@ pub fn navier_stokes_time_derivative(
     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
+    let gradients = SpatialGradients::new(
+        *velocity_gradient,
+        pressure_gradient,
+        Some(laplacian_velocity),
+    );
+    NavierStokes.compute_time_derivative(properties, state, &gradients, body_force_accel)
 }
 
 /// Computes the time evolution of velocity based on the Euler Equations (Inviscid).
@@ -99,11 +163,8 @@ pub fn euler_time_derivative(
     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
+    // Construct dummy properties with the given rho
+    let properties = FluidProperties::new(rho, 0.0); // Viscosity ignored by Euler
+    let gradients = SpatialGradients::new(*velocity_gradient, pressure_gradient, None);
+    Euler.compute_time_derivative(&properties, state, &gradients, body_force_accel)
 }

math_explorer/src/physics/fluid_dynamics/tests.rs

@@ -3,10 +3,11 @@ 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,
         },
         regimes::{FlatPlateClassifier, FlowClassifier, FlowRegime, PipeFlowClassifier},
-        types::{FlowState, FluidProperties},
+        types::{FlowState, FluidProperties, SpatialGradients},
     };
     use nalgebra::{Matrix3, Vector3};
 
@@ -113,4 +114,27 @@ mod tests {
         assert_eq!(continuity_divergence(0.0), 0.0);
         assert_eq!(continuity_divergence(0.5), 0.5);
     }
+
+    #[test]
+    fn test_momentum_strategies() {
+        let props = FluidProperties::new(1.0, 1.0); // rho=1, mu=1
+        let state = FlowState::new(Vector3::new(1.0, 0.0, 0.0), 0.0);
+
+        let vel_grad = Matrix3::zeros();
+        let p_grad = Vector3::zeros();
+        let lap_vel = Vector3::new(1.0, 0.0, 0.0);
+        let g = Vector3::zeros();
+
+        // Check Navier Stokes: should include viscosity
+        // Term = nu * lap_vel = 1.0 * 1.0 = 1.0
+        let gradients_ns = SpatialGradients::new(vel_grad, p_grad, Some(lap_vel));
+        let ns_accel = NavierStokes.compute_time_derivative(&props, &state, &gradients_ns, g);
+        assert!((ns_accel.x - 1.0).abs() < 1e-9);
+
+        // Check Euler: should ignore viscosity even if laplacian is provided
+        let gradients_euler_with_lap = SpatialGradients::new(vel_grad, p_grad, Some(lap_vel));
+        let euler_accel_2 =
+            Euler.compute_time_derivative(&props, &state, &gradients_euler_with_lap, g);
+        assert!((euler_accel_2.x - 0.0).abs() < 1e-9);
+    }
 }

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,28 @@ impl FlowState {
         Self { velocity, pressure }
     }
 }
+
+/// Encapsulates spatial gradients required for momentum equations.
+#[derive(Debug, Clone)]
+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 for inviscid flow.
+    pub laplacian_velocity: Option<Vector3<f64>>,
+}
+
+impl SpatialGradients {
+    pub fn new(
+        velocity_gradient: Matrix3<f64>,
+        pressure_gradient: Vector3<f64>,
+        laplacian_velocity: Option<Vector3<f64>>,
+    ) -> Self {
+        Self {
+            velocity_gradient,
+            pressure_gradient,
+            laplacian_velocity,
+        }
+    }
+}