From b09415e7639fa0ee5a79e62078c0c6fc041935df Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Thu, 29 Jan 2026 22:26:58 +0000 Subject: [PATCH 1/2] Refactor Fluid Dynamics Momentum Equations to use Strategy Pattern. - Introduced `SpatialGradients` parameter object in `types.rs`. - Defined `MomentumEquation` trait in `conservation.rs`. - Implemented `NavierStokes` and `Euler` strategies. - Updated `navier_stokes_time_derivative` and `euler_time_derivative` as backward-compatible wrappers. - Added `FluidError` in `error.rs` for explicit error handling. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- .jules/architect.md | 9 ++ .../physics/fluid_dynamics/conservation.rs | 133 ++++++++++++++---- .../src/physics/fluid_dynamics/error.rs | 25 ++++ .../src/physics/fluid_dynamics/mod.rs | 1 + .../src/physics/fluid_dynamics/types.rs | 15 +- 5 files changed, 156 insertions(+), 27 deletions(-) create mode 100644 math_explorer/src/physics/fluid_dynamics/error.rs diff --git a/.jules/architect.md b/.jules/architect.md index 816d790..ff91537 100644 --- a/.jules/architect.md +++ b/.jules/architect.md @@ -36,3 +36,12 @@ **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-01-29 - [Fluid Dynamics Momentum Strategy] +**Problem:** `math_explorer/src/physics/fluid_dynamics/conservation.rs` contained duplicated logic for Navier-Stokes and Euler equations, with long argument lists (Primitive Obsession) and rigid function signatures preventing easy extension (e.g., adding turbulence models). +**Decision:** Applied "Strategy Pattern" and "Parameter Object". +1. Defined `MomentumEquation` trait in `conservation.rs`. +2. Extracted `SpatialGradients` struct in `types.rs` to group derivatives. +3. Implemented `NavierStokes` and `Euler` strategies. +4. Introduced `FluidError` for robust error handling. +**Consequence:** New models can be added by implementing the trait. Legacy functions are now thin wrappers around the strategies, maintaining backward compatibility while allowing the internal logic to evolve. Inputs are strictly typed and validated. diff --git a/math_explorer/src/physics/fluid_dynamics/conservation.rs b/math_explorer/src/physics/fluid_dynamics/conservation.rs index b0a253f..0329ae7 100644 --- a/math_explorer/src/physics/fluid_dynamics/conservation.rs +++ b/math_explorer/src/physics/fluid_dynamics/conservation.rs @@ -2,9 +2,93 @@ //! //! 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; +/// Defines a strategy for computing the momentum time derivative. +pub trait MomentumEquation { + /// Computes $\frac{\partial \mathbf{u}}{\partial t}$ given the current state and gradients. + fn compute_time_derivative( + &self, + properties: &FluidProperties, + state: &FlowState, + gradients: &SpatialGradients, + body_force_accel: Vector3, + ) -> Result, FluidError>; +} + +/// Navier-Stokes Momentum Equation (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}$$ +pub struct NavierStokes; + +impl MomentumEquation for NavierStokes { + fn compute_time_derivative( + &self, + properties: &FluidProperties, + state: &FlowState, + gradients: &SpatialGradients, + body_force_accel: Vector3, + ) -> Result, FluidError> { + let laplacian = gradients + .laplacian_velocity + .ok_or(FluidError::MissingLaplacian)?; + + let nu = properties.kinematic_viscosity(); + let rho = properties.density; + + if rho <= 0.0 { + return Err(FluidError::InvalidProperties( + "Density must be positive".to_string(), + )); + } + + // 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 Momentum Equation (Inviscid Flow). +/// +/// $$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \mathbf{g}$$ +pub struct Euler; + +impl MomentumEquation for Euler { + fn compute_time_derivative( + &self, + properties: &FluidProperties, + state: &FlowState, + gradients: &SpatialGradients, + body_force_accel: Vector3, + ) -> Result, FluidError> { + let rho = properties.density; + + if rho <= 0.0 { + return Err(FluidError::InvalidProperties( + "Density must be positive".to_string(), + )); + } + + // 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) + } +} + /// Calculates the Material Derivative ($D/Dt$) of a scalar property. /// /// $$\frac{D\phi}{Dt} = \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi$$ @@ -33,10 +117,6 @@ pub fn material_derivative_vector( gradient_tensor: &nalgebra::Matrix3, ) -> Vector3 { // (\mathbf{u} \cdot \nabla) \mathbf{A} corresponds to Jacobian * velocity vector - // J = [ dA_x/dx dA_x/dy dA_x/dz ] - // [ dA_y/dx dA_y/dy dA_y/dz ] - // [ ... ] - // Result is J * u local_change + gradient_tensor * velocity } @@ -60,9 +140,7 @@ pub fn continuity_divergence(velocity_divergence: f64) -> f64 { /// * `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). +/// * `body_force_accel`: External forces (e.g., gravity $\mathbf{g}$). pub fn navier_stokes_time_derivative( properties: &FluidProperties, state: &FlowState, @@ -71,20 +149,16 @@ pub fn navier_stokes_time_derivative( 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 + let solver = NavierStokes; + solver + .compute_time_derivative(properties, state, &gradients, body_force_accel) + .expect("Navier-Stokes calculation failed (likely invalid properties or missing data)") } /// Computes the time evolution of velocity based on the Euler Equations (Inviscid). @@ -99,11 +173,18 @@ pub fn euler_time_derivative( pressure_gradient: Vector3, body_force_accel: Vector3, ) -> Vector3 { - // Convective term: -(u . del) u - let convection = -(velocity_gradient * state.velocity); + let gradients = SpatialGradients { + velocity_gradient: *velocity_gradient, + pressure_gradient, + laplacian_velocity: None, + }; - // Pressure term: -(1/rho) grad p - let pressure_term = -pressure_gradient / rho; + // Euler only needs density, but the strategy takes FluidProperties. + // We construct a dummy property struct with density and 0 viscosity. + let properties = FluidProperties::new(rho, 0.0); - convection + pressure_term + body_force_accel + let solver = Euler; + solver + .compute_time_derivative(&properties, state, &gradients, body_force_accel) + .expect("Euler calculation failed (likely invalid properties)") } 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..9478b1f --- /dev/null +++ b/math_explorer/src/physics/fluid_dynamics/error.rs @@ -0,0 +1,25 @@ +//! Errors for Fluid Dynamics. + +use std::error::Error; +use std::fmt; + +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum FluidError { + /// Indicates that the fluid properties are invalid (e.g., negative density). + InvalidProperties(String), + /// Indicates that the Laplacian of the velocity field is required but missing. + MissingLaplacian, +} + +impl fmt::Display for FluidError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + FluidError::InvalidProperties(msg) => write!(f, "Invalid Fluid Properties: {}", msg), + FluidError::MissingLaplacian => { + write!(f, "Missing Laplacian of velocity (required for Viscous flow)") + } + } + } +} + +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/types.rs b/math_explorer/src/physics/fluid_dynamics/types.rs index 100132c..3245a51 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,16 @@ impl FlowState { Self { velocity, pressure } } } + +/// Encapsulates spatial derivatives needed for momentum equations. +/// +/// Prevents "Primitive Obsession" by grouping related gradients. +#[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 for Inviscid flow. + pub laplacian_velocity: Option>, +} From ee530bc0f311a384c4d063e07ce462317cb4d979 Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Thu, 29 Jan 2026 22:33:26 +0000 Subject: [PATCH 2/2] Refactor Fluid Dynamics Momentum Equations to use Strategy Pattern. - Introduced `SpatialGradients` parameter object in `types.rs`. - Defined `MomentumEquation` trait in `conservation.rs`. - Implemented `NavierStokes` and `Euler` strategies. - Updated `navier_stokes_time_derivative` and `euler_time_derivative` as backward-compatible wrappers. - Added `FluidError` in `error.rs` for explicit error handling. - Ensured formatting compliance with `cargo fmt`. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- math_explorer/src/physics/fluid_dynamics/error.rs | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/math_explorer/src/physics/fluid_dynamics/error.rs b/math_explorer/src/physics/fluid_dynamics/error.rs index 9478b1f..b86f6b4 100644 --- a/math_explorer/src/physics/fluid_dynamics/error.rs +++ b/math_explorer/src/physics/fluid_dynamics/error.rs @@ -16,7 +16,10 @@ impl fmt::Display for FluidError { match self { FluidError::InvalidProperties(msg) => write!(f, "Invalid Fluid Properties: {}", msg), FluidError::MissingLaplacian => { - write!(f, "Missing Laplacian of velocity (required for Viscous flow)") + write!( + f, + "Missing Laplacian of velocity (required for Viscous flow)" + ) } } }