diff --git a/.jules/mason.md b/.jules/mason.md index 932a4f6..0b13074 100644 --- a/.jules/mason.md +++ b/.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. diff --git a/math_explorer/src/physics/fluid_dynamics/conservation.rs b/math_explorer/src/physics/fluid_dynamics/conservation.rs index b0a253f..51333b8 100644 --- a/math_explorer/src/physics/fluid_dynamics/conservation.rs +++ b/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, velocity: Vector3, - gradient_tensor: &nalgebra::Matrix3, + gradient_tensor: &Matrix3, ) -> Vector3 { // (\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, + ) -> Result, 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, + ) -> Result, 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, + ) -> Result, 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, gradient: &Matrix3) -> Vector3 { + // -(u . del) u + -(gradient * velocity) +} + +fn pressure_acceleration(pressure_gradient: &Vector3, density: f64) -> Vector3 { + // -(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,28 +138,24 @@ 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, + velocity_gradient: &Matrix3, pressure_gradient: Vector3, laplacian_velocity: Vector3, body_force_accel: Vector3, ) -> Vector3 { - 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). @@ -92,18 +163,25 @@ pub fn navier_stokes_time_derivative( /// $$\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, + velocity_gradient: &Matrix3, pressure_gradient: Vector3, body_force_accel: Vector3, ) -> Vector3 { - // 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") } diff --git a/math_explorer/src/physics/fluid_dynamics/error.rs b/math_explorer/src/physics/fluid_dynamics/error.rs new file mode 100644 index 0000000..681bd61 --- /dev/null +++ b/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 {} diff --git a/math_explorer/src/physics/fluid_dynamics/mod.rs b/math_explorer/src/physics/fluid_dynamics/mod.rs index a00c62a..af610c8 100644 --- a/math_explorer/src/physics/fluid_dynamics/mod.rs +++ b/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; diff --git a/math_explorer/src/physics/fluid_dynamics/tests.rs b/math_explorer/src/physics/fluid_dynamics/tests.rs index c2bf63d..510931d 100644 --- a/math_explorer/src/physics/fluid_dynamics/tests.rs +++ b/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); + } } diff --git a/math_explorer/src/physics/fluid_dynamics/types.rs b/math_explorer/src/physics/fluid_dynamics/types.rs index 100132c..b21d59e 100644 --- a/math_explorer/src/physics/fluid_dynamics/types.rs +++ b/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, + /// Gradient of pressure ($\nabla p$). + pub pressure_gradient: Vector3, + /// Laplacian of velocity ($\nabla^2 \mathbf{u}$). + /// Optional because it is not required for Euler (inviscid) flow. + pub laplacian_velocity: Option>, +} + +impl SpatialGradients { + /// Creates a new `SpatialGradients` struct. + pub fn new( + velocity_gradient: nalgebra::Matrix3, + pressure_gradient: Vector3, + laplacian_velocity: Option>, + ) -> Self { + Self { + velocity_gradient, + pressure_gradient, + laplacian_velocity, + } + } +}