diff --git a/.jules/systems_core.md b/.jules/systems_core.md index 183b3ca..c0a1c39 100644 --- a/.jules/systems_core.md +++ b/.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). diff --git a/math_explorer/src/physics/fluid_dynamics/conservation.rs b/math_explorer/src/physics/fluid_dynamics/conservation.rs index b0a253f..e57e2e4 100644 --- a/math_explorer/src/physics/fluid_dynamics/conservation.rs +++ b/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,20 +50,84 @@ 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, + ) -> Result, 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, + ) -> Result, 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, + ) -> Result, 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, @@ -70,40 +135,27 @@ pub fn navier_stokes_time_derivative( 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; - - // Sum - convection + pressure_term + viscous_term + body_force_accel +) -> Result, 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, pressure_gradient: Vector3, body_force_accel: Vector3, -) -> Vector3 { - // 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, 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) } 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..9360b2c --- /dev/null +++ b/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 {} 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..d785793 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, 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))); + } } diff --git a/math_explorer/src/physics/fluid_dynamics/types.rs b/math_explorer/src/physics/fluid_dynamics/types.rs index 100132c..c7be0df 100644 --- a/math_explorer/src/physics/fluid_dynamics/types.rs +++ b/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, + /// Gradient of pressure ($\nabla p$). + pub pressure_gradient: Vector3, + /// Laplacian of velocity ($\nabla^2 \mathbf{u}$). Optional, as Euler equations don't use it. + pub laplacian_velocity: Option>, +} + +impl SpatialGradients { + /// Creates a new `SpatialGradients` structure. + pub fn new( + velocity_gradient: Matrix3, + pressure_gradient: Vector3, + laplacian_velocity: Option>, + ) -> Self { + Self { + velocity_gradient, + pressure_gradient, + laplacian_velocity, + } + } +}