🧱 Mason: Refactor Momentum Equations (Fixes OCP/Primitive Obsession)

math_explorer/src/physics/fluid_dynamics/conservation.rs

@@ -2,7 +2,7 @@
 //!
 //! Implements the core Partial Differential Equations (PDEs) of Fluid Dynamics.
 
-use super::types::{FlowState, FluidProperties};
+use super::types::{FlowState, FluidError, FluidProperties, SpatialGradients};
 use nalgebra::Vector3;
 
 /// Calculates the Material Derivative ($D/Dt$) of a scalar property.
@@ -49,61 +49,108 @@ pub fn continuity_divergence(velocity_divergence: f64) -> f64 {
     velocity_divergence
 }
 
+/// Strategy for calculating the time evolution of velocity (acceleration).
+pub trait MomentumEquation {
+    fn calculate_acceleration(
+        &self,
+        properties: &FluidProperties,
+        state: &FlowState,
+        gradients: &SpatialGradients,
+        body_force_accel: Vector3<f64>,
+    ) -> Result<Vector3<f64>, FluidError>;
+}
+
+/// Navier-Stokes Equations (Viscous Flow).
+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> {
+        if properties.density <= 0.0 {
+            return Err(FluidError::InvalidDensity(
+                "Density must be positive".to_string(),
+            ));
+        }
+        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 = gradients.laplacian_velocity * nu;
+
+        Ok(convection + pressure_term + viscous_term + body_force_accel)
+    }
+}
+
+/// Euler Equations (Inviscid Flow).
+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> {
+        if properties.density <= 0.0 {
+            return Err(FluidError::InvalidDensity(
+                "Density must be positive".to_string(),
+            ));
+        }
+        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.
 ///
-/// $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{g}$$
-///
-/// Returns the local acceleration $\frac{\partial \mathbf{u}}{\partial t}$.
-///
-/// * `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).
+/// # Returns
+/// * `Result<Vector3<f64>, FluidError>` - The local acceleration.
 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, 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$.
+/// # Returns
+/// * `Result<Vector3<f64>, FluidError>` - The local acceleration.
 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> {
+    // Construct simplified properties (viscosity irrelevant for Euler)
+    let properties = FluidProperties::new(rho, 0.0);
+    // Laplacian is irrelevant for Euler
+    let gradients = SpatialGradients::new(*velocity_gradient, pressure_gradient, Vector3::zeros());
+    Euler.calculate_acceleration(&properties, state, &gradients, body_force_accel)
 }

math_explorer/src/physics/fluid_dynamics/tests.rs

@@ -99,12 +99,13 @@ 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).unwrap();
         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).unwrap();
         assert!((accel_p.x - (-2.0)).abs() < 1e-9);
     }
 

math_explorer/src/physics/fluid_dynamics/types.rs

@@ -1,6 +1,13 @@
 //! Types for Fluid Dynamics.
 
-use nalgebra::Vector3;
+use nalgebra::{Matrix3, Vector3};
+
+/// Errors that can occur during fluid dynamics calculations.
+#[derive(Debug, Clone)]
+pub enum FluidError {
+    InvalidDensity(String),
+    NumericalInstability(String),
+}
 
 /// Physical properties of the fluid.
 #[derive(Debug, Clone, Copy)]
@@ -56,3 +63,31 @@ impl FlowState {
         Self { velocity, pressure }
     }
 }
+
+/// Container for spatial derivatives required by momentum equations.
+///
+/// Encapsulates the various gradients needed to compute flow acceleration,
+/// reducing parameter bloat in solver functions.
+#[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}$).
+    pub laplacian_velocity: Vector3<f64>,
+}
+
+impl SpatialGradients {
+    pub fn new(
+        velocity_gradient: Matrix3<f64>,
+        pressure_gradient: Vector3<f64>,
+        laplacian_velocity: Vector3<f64>,
+    ) -> Self {
+        Self {
+            velocity_gradient,
+            pressure_gradient,
+            laplacian_velocity,
+        }
+    }
+}