From 81583bbcc9d2614ced6dc4241f8ac0dade70e795 Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Sat, 31 Jan 2026 00:31:30 +0000 Subject: [PATCH 1/2] Refactor fluid dynamics momentum equations to use Strategy Pattern. This commit introduces the `MomentumEquation` trait to decouple the flow regime (Navier-Stokes vs Euler) from the solver logic. It also introduces a `SpatialGradients` parameter object to encapsulate the various gradient terms, reducing parameter bloat. Key changes: - Added `SpatialGradients` and `FluidError` to `types.rs`. - Defined `MomentumEquation` trait in `conservation.rs`. - Implemented `NavierStokes` and `Euler` strategies. - Updated `navier_stokes_time_derivative` and `euler_time_derivative` to return `Result` and use the strategies. - Updated tests to handle `Result` propagation. This addresses OCP and Primitive Obsession violations. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- .../physics/fluid_dynamics/conservation.rs | 128 ++++++++++++------ .../src/physics/fluid_dynamics/tests.rs | 4 +- .../src/physics/fluid_dynamics/types.rs | 37 ++++- 3 files changed, 127 insertions(+), 42 deletions(-) diff --git a/math_explorer/src/physics/fluid_dynamics/conservation.rs b/math_explorer/src/physics/fluid_dynamics/conservation.rs index b0a253f..44e9750 100644 --- a/math_explorer/src/physics/fluid_dynamics/conservation.rs +++ b/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,20 +49,81 @@ 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, + ) -> Result, 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, + ) -> Result, 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, + ) -> Result, 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, FluidError>` - The local acceleration. pub fn navier_stokes_time_derivative( properties: &FluidProperties, state: &FlowState, @@ -70,40 +131,29 @@ 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, 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, FluidError>` - The local acceleration. 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> { + // 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) } diff --git a/math_explorer/src/physics/fluid_dynamics/tests.rs b/math_explorer/src/physics/fluid_dynamics/tests.rs index c2bf63d..a05ea8d 100644 --- a/math_explorer/src/physics/fluid_dynamics/tests.rs +++ b/math_explorer/src/physics/fluid_dynamics/tests.rs @@ -99,12 +99,12 @@ 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); } diff --git a/math_explorer/src/physics/fluid_dynamics/types.rs b/math_explorer/src/physics/fluid_dynamics/types.rs index 100132c..65eab27 100644 --- a/math_explorer/src/physics/fluid_dynamics/types.rs +++ b/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, + /// Gradient of pressure ($\nabla p$). + pub pressure_gradient: Vector3, + /// Laplacian of velocity ($\nabla^2 \mathbf{u}$). + pub laplacian_velocity: Vector3, +} + +impl SpatialGradients { + pub fn new( + velocity_gradient: Matrix3, + pressure_gradient: Vector3, + laplacian_velocity: Vector3, + ) -> Self { + Self { + velocity_gradient, + pressure_gradient, + laplacian_velocity, + } + } +} From 3813b00800786f137023c53479bb1d70ef557864 Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Sat, 31 Jan 2026 00:36:08 +0000 Subject: [PATCH 2/2] Refactor fluid dynamics momentum equations to use Strategy Pattern. This commit introduces the `MomentumEquation` trait to decouple the flow regime (Navier-Stokes vs Euler) from the solver logic. It also introduces a `SpatialGradients` parameter object to encapsulate the various gradient terms, reducing parameter bloat. Key changes: - Added `SpatialGradients` and `FluidError` to `types.rs`. - Defined `MomentumEquation` trait in `conservation.rs`. - Implemented `NavierStokes` and `Euler` strategies. - Updated `navier_stokes_time_derivative` and `euler_time_derivative` to return `Result` and use the strategies. - Updated tests to handle `Result` propagation. - Formatted code with `cargo fmt`. This addresses OCP and Primitive Obsession violations. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- math_explorer/src/physics/fluid_dynamics/conservation.rs | 9 +++------ math_explorer/src/physics/fluid_dynamics/tests.rs | 3 ++- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/math_explorer/src/physics/fluid_dynamics/conservation.rs b/math_explorer/src/physics/fluid_dynamics/conservation.rs index 44e9750..f37bd08 100644 --- a/math_explorer/src/physics/fluid_dynamics/conservation.rs +++ b/math_explorer/src/physics/fluid_dynamics/conservation.rs @@ -132,7 +132,8 @@ pub fn navier_stokes_time_derivative( laplacian_velocity: Vector3, body_force_accel: Vector3, ) -> Result, FluidError> { - let gradients = SpatialGradients::new(*velocity_gradient, pressure_gradient, laplacian_velocity); + let gradients = + SpatialGradients::new(*velocity_gradient, pressure_gradient, laplacian_velocity); NavierStokes.calculate_acceleration(properties, state, &gradients, body_force_accel) } @@ -150,10 +151,6 @@ pub fn euler_time_derivative( // 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(), - ); + let gradients = SpatialGradients::new(*velocity_gradient, pressure_gradient, Vector3::zeros()); Euler.calculate_acceleration(&properties, state, &gradients, body_force_accel) } diff --git a/math_explorer/src/physics/fluid_dynamics/tests.rs b/math_explorer/src/physics/fluid_dynamics/tests.rs index a05ea8d..5097be3 100644 --- a/math_explorer/src/physics/fluid_dynamics/tests.rs +++ b/math_explorer/src/physics/fluid_dynamics/tests.rs @@ -99,7 +99,8 @@ 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).unwrap(); + 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);