🏛️ Architect: Refactor Fluid Dynamics Momentum Equations

.jules/architect.md

@@ -36,3 +36,12 @@
 **Problem:** The Reaction-Diffusion stencil logic (finite difference Laplacian + kinetics) was duplicated in both `step` (optimized Euler) and `derivative_in_place` (generic ODE support) to preserve performance, leading to a DRY violation.
 **Decision:** Extracted the stencil logic into a private `apply_reaction_diffusion_stencil` method accepting a closure.
 **Consequence:** Eliminates code duplication while preserving the performance of the hot loop (unsafe access, unrolling).
+
+## 2026-01-29 - [Fluid Dynamics Momentum Strategy]
+**Problem:** `math_explorer/src/physics/fluid_dynamics/conservation.rs` contained duplicated logic for Navier-Stokes and Euler equations, with long argument lists (Primitive Obsession) and rigid function signatures preventing easy extension (e.g., adding turbulence models).
+**Decision:** Applied "Strategy Pattern" and "Parameter Object".
+1. Defined `MomentumEquation` trait in `conservation.rs`.
+2. Extracted `SpatialGradients` struct in `types.rs` to group derivatives.
+3. Implemented `NavierStokes` and `Euler` strategies.
+4. Introduced `FluidError` for robust error handling.
+**Consequence:** New models can be added by implementing the trait. Legacy functions are now thin wrappers around the strategies, maintaining backward compatibility while allowing the internal logic to evolve. Inputs are strictly typed and validated.

math_explorer/src/physics/fluid_dynamics/conservation.rs

@@ -2,9 +2,93 @@
 //!
 //! 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;
 
+/// Defines a strategy for computing the momentum time derivative.
+pub trait MomentumEquation {
+    /// Computes $\frac{\partial \mathbf{u}}{\partial t}$ given the current state and gradients.
+    fn compute_time_derivative(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError>;
+}
+
+/// 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}$$
+pub struct NavierStokes;
+
+impl MomentumEquation for NavierStokes {
+    fn compute_time_derivative(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError> {
+        let laplacian = gradients
+            .laplacian_velocity
+            .ok_or(FluidError::MissingLaplacian)?;
+
+        let nu = properties.kinematic_viscosity();
+        let rho = properties.density;
+
+        if rho <= 0.0 {
+            return Err(FluidError::InvalidProperties(
+                "Density must be positive".to_string(),
+            ));
+        }
+
+        // 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 Momentum Equation (Inviscid Flow).
+///
+/// $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \mathbf{g}$$
+pub struct Euler;
+
+impl MomentumEquation for Euler {
+    fn compute_time_derivative(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError> {
+        let rho = properties.density;
+
+        if rho <= 0.0 {
+            return Err(FluidError::InvalidProperties(
+                "Density must be positive".to_string(),
+            ));
+        }
+
+        // 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)
+    }
+}
+
 /// Calculates the Material Derivative ($D/Dt$) of a scalar property.
 ///
 /// $$\frac{D\phi}{Dt} = \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi$$
@@ -33,10 +117,6 @@ pub fn material_derivative_vector(
     gradient_tensor: &nalgebra::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 ]
-    //     [ dA_y/dx  dA_y/dy  dA_y/dz ]
-    //     [ ...                   ]
-    // Result is J * u
     local_change + gradient_tensor * velocity
 }
 
@@ -60,9 +140,7 @@ pub fn continuity_divergence(velocity_divergence: f64) -> f64 {
 /// * `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).
+/// * `body_force_accel`: External forces (e.g., gravity $\mathbf{g}$).
 pub fn navier_stokes_time_derivative(
     properties: &FluidProperties,
     state: &FlowState,
@@ -71,20 +149,16 @@ 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;
+    let gradients = SpatialGradients {
+        velocity_gradient: *velocity_gradient,
+        pressure_gradient,
+        laplacian_velocity: Some(laplacian_velocity),
+    };
 
-    // Sum
-    convection + pressure_term + viscous_term + body_force_accel
+    let solver = NavierStokes;
+    solver
+        .compute_time_derivative(properties, state, &gradients, body_force_accel)
+        .expect("Navier-Stokes calculation failed (likely invalid properties or missing data)")
 }
 
 /// Computes the time evolution of velocity based on the Euler Equations (Inviscid).
@@ -99,11 +173,18 @@ 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);
+    let gradients = SpatialGradients {
+        velocity_gradient: *velocity_gradient,
+        pressure_gradient,
+        laplacian_velocity: None,
+    };
 
-    // Pressure term: -(1/rho) grad p
-    let pressure_term = -pressure_gradient / rho;
+    // Euler only needs density, but the strategy takes FluidProperties.
+    // We construct a dummy property struct with density and 0 viscosity.
+    let properties = FluidProperties::new(rho, 0.0);
 
-    convection + pressure_term + body_force_accel
+    let solver = Euler;
+    solver
+        .compute_time_derivative(&properties, state, &gradients, body_force_accel)
+        .expect("Euler calculation failed (likely invalid properties)")
 }

math_explorer/src/physics/fluid_dynamics/error.rs

@@ -0,0 +1,28 @@
+//! Errors for Fluid Dynamics.
+
+use std::error::Error;
+use std::fmt;
+
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub enum FluidError {
+    /// Indicates that the fluid properties are invalid (e.g., negative density).
+    InvalidProperties(String),
+    /// Indicates that the Laplacian of the velocity field is required but missing.
+    MissingLaplacian,
+}
+
+impl fmt::Display for FluidError {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            FluidError::InvalidProperties(msg) => write!(f, "Invalid Fluid Properties: {}", msg),
+            FluidError::MissingLaplacian => {
+                write!(
+                    f,
+                    "Missing Laplacian of velocity (required for Viscous flow)"
+                )
+            }
+        }
+    }
+}
+
+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/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,16 @@ impl FlowState {
         Self { velocity, pressure }
     }
 }
+
+/// Encapsulates spatial derivatives needed for momentum equations.
+///
+/// Prevents "Primitive Obsession" by grouping related gradients.
+#[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 for Inviscid flow.
+    pub laplacian_velocity: Option<Vector3<f64>>,
+}