From 7340e65fc7a8e0a05f7b346ac8bd9c93e08b4798 Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Fri, 30 Jan 2026 22:13:10 +0000 Subject: [PATCH] Refactor fluid dynamics conservation to use Strategy Pattern - Introduced `FluidError` for robust error handling. - Introduced `SpatialGradients` parameter object. - Defined `MomentumEquation` trait with `NavierStokes` and `Euler` strategies. - Refactored `navier_stokes_time_derivative` and `euler_time_derivative` to use the new strategies and return `Result`. - Updated tests and added ADR. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- .jules/architect.md | 5 + .../physics/fluid_dynamics/conservation.rs | 135 ++++++++++++++---- .../src/physics/fluid_dynamics/error.rs | 29 ++++ .../src/physics/fluid_dynamics/mod.rs | 1 + .../src/physics/fluid_dynamics/tests.rs | 42 +++++- .../src/physics/fluid_dynamics/types.rs | 32 ++++- 6 files changed, 211 insertions(+), 33 deletions(-) create mode 100644 math_explorer/src/physics/fluid_dynamics/error.rs diff --git a/.jules/architect.md b/.jules/architect.md index 816d790..a4a5965 100644 --- a/.jules/architect.md +++ b/.jules/architect.md @@ -36,3 +36,8 @@ **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-06-25 - [Fluid Dynamics Conservation Strategy] +**Problem:** `math_explorer/src/physics/fluid_dynamics/conservation.rs` contained loosely coupled functions `navier_stokes_time_derivative` and `euler_time_derivative` with "Data Clump" signatures (passing multiple gradient vectors individually) and duplicated logic. Error handling (e.g., zero density, missing Laplacian) was absent. +**Decision:** Applied **Strategy Pattern** and **Parameter Object**. Defined `SpatialGradients` struct to encapsulate derivatives. Defined `MomentumEquation` trait with `NavierStokes` and `Euler` strategies. Introduced `FluidError` for robust error handling. +**Consequence:** The API is now type-safe and extensible. New flow regimes can be added by implementing the trait. Error states are explicitly handled via `Result`. diff --git a/math_explorer/src/physics/fluid_dynamics/conservation.rs b/math_explorer/src/physics/fluid_dynamics/conservation.rs index b0a253f..8962646 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,95 @@ pub fn continuity_divergence(velocity_divergence: f64) -> f64 { velocity_divergence } +/// Strategy for calculating momentum conservation (flow acceleration). +/// +/// This trait allows for swapping between different flow regimes (e.g., Navier-Stokes, Euler) +/// without changing the integration logic. +pub trait MomentumEquation { + /// Computes the local acceleration $\frac{\partial \mathbf{u}}{\partial t}$. + /// + /// # Arguments + /// * `properties` - Fluid density and viscosity. + /// * `state` - Current velocity and pressure. + /// * `grads` - Spatial gradients of velocity and pressure. + /// * `body_force` - External body force (acceleration, e.g., gravity). + fn compute_acceleration( + &self, + properties: &FluidProperties, + state: &FlowState, + grads: &SpatialGradients, + body_force: Vector3, + ) -> Result, FluidError>; +} + +/// Navier-Stokes momentum equation for viscous flow. +/// +/// Requires the Laplacian of velocity. +#[derive(Debug, Clone, Copy, Default)] +pub struct NavierStokes; + +impl MomentumEquation for NavierStokes { + fn compute_acceleration( + &self, + properties: &FluidProperties, + state: &FlowState, + grads: &SpatialGradients, + body_force: Vector3, + ) -> Result, FluidError> { + if properties.density.abs() < f64::EPSILON { + return Err(FluidError::ZeroDensity); + } + + let laplacian = grads + .laplacian_velocity + .ok_or(FluidError::MissingLaplacian)?; + + let nu = properties.kinematic_viscosity(); + let rho = properties.density; + + // Convective term: -(u . del) u + let convection = -(grads.velocity_gradient * state.velocity); + + // Pressure term: -(1/rho) grad p + let pressure_term = -grads.pressure_gradient / rho; + + // Viscous term: nu * del^2 u + let viscous_term = laplacian * nu; + + Ok(convection + pressure_term + viscous_term + body_force) + } +} + +/// Euler momentum equation for inviscid flow. +/// +/// Ignores viscosity and velocity Laplacian. +#[derive(Debug, Clone, Copy, Default)] +pub struct Euler; + +impl MomentumEquation for Euler { + fn compute_acceleration( + &self, + properties: &FluidProperties, + state: &FlowState, + grads: &SpatialGradients, + body_force: Vector3, + ) -> Result, FluidError> { + if properties.density.abs() < f64::EPSILON { + return Err(FluidError::ZeroDensity); + } + + let rho = properties.density; + + // Convective term: -(u . del) u + let convection = -(grads.velocity_gradient * state.velocity); + + // Pressure term: -(1/rho) grad p + let pressure_term = -grads.pressure_gradient / rho; + + Ok(convection + pressure_term + body_force) + } +} + /// 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}$$ @@ -66,25 +156,17 @@ pub fn continuity_divergence(velocity_divergence: f64) -> f64 { 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; - - // Sum - convection + pressure_term + viscous_term + body_force_accel +) -> Result, FluidError> { + let grads = SpatialGradients::new( + *velocity_gradient, + pressure_gradient, + Some(laplacian_velocity), + ); + NavierStokes.compute_acceleration(properties, state, &grads, body_force_accel) } /// Computes the time evolution of velocity based on the Euler Equations (Inviscid). @@ -95,15 +177,12 @@ pub fn navier_stokes_time_derivative( 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); - - // Pressure term: -(1/rho) grad p - let pressure_term = -pressure_gradient / rho; - - convection + pressure_term + body_force_accel +) -> Result, FluidError> { + // Create dummy properties for Euler (viscosity not used) + let properties = FluidProperties::new(rho, 0.0); + let grads = SpatialGradients::new(*velocity_gradient, pressure_gradient, None); + Euler.compute_acceleration(&properties, state, &grads, 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..8648691 --- /dev/null +++ b/math_explorer/src/physics/fluid_dynamics/error.rs @@ -0,0 +1,29 @@ +use std::fmt; + +/// Errors that can occur during fluid dynamics calculations. +#[derive(Debug, Clone, PartialEq)] +pub enum FluidError { + /// The Laplacian of the velocity field is required but was not provided. + MissingLaplacian, + /// Fluid density is zero, causing division by zero. + ZeroDensity, + /// Invalid configuration of fluid properties or state. + InvalidConfiguration(String), +} + +impl fmt::Display for FluidError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + FluidError::MissingLaplacian => { + write!( + f, + "Navier-Stokes equation requires the Laplacian of the velocity field." + ) + } + FluidError::ZeroDensity => write!(f, "Fluid density cannot be zero."), + FluidError::InvalidConfiguration(msg) => write!(f, "Invalid 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..c683632 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}; @@ -99,12 +101,15 @@ 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); + // Update: Expecting Result + 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 +118,33 @@ 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); + let state = FlowState::new(Vector3::zeros(), 0.0); + let g = Vector3::zeros(); + + // 1. Test Navier-Stokes Missing Laplacian + let grads_no_lap = SpatialGradients::new(Matrix3::zeros(), Vector3::zeros(), None); + let ns_result = NavierStokes.compute_acceleration(&props, &state, &grads_no_lap, g); + assert_eq!(ns_result, Err(FluidError::MissingLaplacian)); + + // 2. Test Euler ignores Laplacian (works even if None) + let euler_result = Euler.compute_acceleration(&props, &state, &grads_no_lap, g); + assert!(euler_result.is_ok()); + + // 3. Test Zero Density + let bad_props = FluidProperties::new(0.0, 1.0); + let grads = + SpatialGradients::new(Matrix3::zeros(), Vector3::zeros(), Some(Vector3::zeros())); + assert_eq!( + NavierStokes.compute_acceleration(&bad_props, &state, &grads, g), + Err(FluidError::ZeroDensity) + ); + assert_eq!( + Euler.compute_acceleration(&bad_props, &state, &grads, g), + Err(FluidError::ZeroDensity) + ); + } } diff --git a/math_explorer/src/physics/fluid_dynamics/types.rs b/math_explorer/src/physics/fluid_dynamics/types.rs index 100132c..96d9dcb 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,33 @@ impl FlowState { Self { velocity, pressure } } } + +/// Encapsulates spatial derivatives of the flow field. +/// +/// This struct groups the gradient terms required for momentum conservation equations, +/// reducing parameter clutter and improving type safety. +#[derive(Debug, Clone, Copy)] +pub struct SpatialGradients { + /// Jacobian of the velocity field ($\nabla \mathbf{u}$). + pub velocity_gradient: Matrix3, + /// Gradient of the pressure field ($\nabla p$). + pub pressure_gradient: Vector3, + /// Laplacian of the velocity field ($\nabla^2 \mathbf{u}$). + /// Only required for viscous flows (Navier-Stokes). + pub laplacian_velocity: Option>, +} + +impl SpatialGradients { + /// Creates a new `SpatialGradients`. + pub fn new( + velocity_gradient: Matrix3, + pressure_gradient: Vector3, + laplacian_velocity: Option>, + ) -> Self { + Self { + velocity_gradient, + pressure_gradient, + laplacian_velocity, + } + } +}