From cc626209dfd41aab780bf6d50efb03560ec61759 Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Thu, 29 Jan 2026 19:23:06 +0000 Subject: [PATCH 1/2] Refactor fluid dynamics conservation to use Strategy Pattern - Introduced `MomentumEquation` trait to define momentum conservation strategies. - Implemented `NavierStokes` and `Euler` strategies, adhering to OCP. - Introduced `SpatialGradients` struct to encapsulate spatial derivatives, solving "Tuple Soup". - Refactored `navier_stokes_time_derivative` and `euler_time_derivative` to use these strategies internally, preserving backward compatibility. - Added regression tests verifying strategy implementations. - Updated Systems Core Engineering Log. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- .jules/systems_core.md | 5 + .../physics/fluid_dynamics/conservation.rs | 109 ++++++++++++++---- .../src/physics/fluid_dynamics/tests.rs | 25 +++- .../src/physics/fluid_dynamics/types.rs | 27 ++++- 4 files changed, 142 insertions(+), 24 deletions(-) diff --git a/.jules/systems_core.md b/.jules/systems_core.md index 183b3ca..9a2c8b3 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-04-10 - Strategy Pattern for Fluid Momentum +**Violation:** `conservation.rs` contained duplicate logic for Navier-Stokes and Euler time derivatives, with hardcoded physics regimes (viscous vs inviscid) and parameter soup (`velocity_gradient`, `pressure_gradient`, `laplacian`), violating OCP and DRY. +**Refactor:** Applied the **Strategy Pattern**. Extracted `MomentumEquation` trait. Implemented `NavierStokes` and `Euler` strategies. Introduced `SpatialGradients` struct to encapsulate spatial derivatives. +**Trade-off:** Minimal boilerplate for trait definitions, but achieved full separation of physics models and eliminated argument soup. diff --git a/math_explorer/src/physics/fluid_dynamics/conservation.rs b/math_explorer/src/physics/fluid_dynamics/conservation.rs index b0a253f..0d6e265 100644 --- a/math_explorer/src/physics/fluid_dynamics/conservation.rs +++ b/math_explorer/src/physics/fluid_dynamics/conservation.rs @@ -2,9 +2,81 @@ //! //! Implements the core Partial Differential Equations (PDEs) of Fluid Dynamics. -use super::types::{FlowState, FluidProperties}; +use super::types::{FlowState, FluidProperties, SpatialGradients}; use nalgebra::Vector3; +/// Defines the strategy for computing the momentum conservation term. +pub trait MomentumEquation { + /// Computes the local acceleration $\frac{\partial \mathbf{u}}{\partial t}$. + fn compute_time_derivative( + &self, + properties: &FluidProperties, + state: &FlowState, + gradients: &SpatialGradients, + body_force_accel: Vector3, + ) -> Vector3; +} + +/// 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}$$ +#[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_accel: Vector3, + ) -> Vector3 { + 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 = match gradients.laplacian_velocity { + Some(laplacian) => laplacian * nu, + None => Vector3::zeros(), + }; + + // Sum + 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}$$ +#[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_accel: Vector3, + ) -> Vector3 { + 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; + + 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$$ @@ -71,20 +143,12 @@ 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; - - // Sum - convection + pressure_term + viscous_term + body_force_accel + let gradients = SpatialGradients::new( + *velocity_gradient, + pressure_gradient, + Some(laplacian_velocity), + ); + NavierStokes.compute_time_derivative(properties, state, &gradients, body_force_accel) } /// Computes the time evolution of velocity based on the Euler Equations (Inviscid). @@ -99,11 +163,12 @@ pub fn euler_time_derivative( 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 + // Construct dummy properties with the given rho + let properties = FluidProperties::new(rho, 0.0); // Viscosity ignored by Euler + let gradients = SpatialGradients::new( + *velocity_gradient, + pressure_gradient, + None, + ); + Euler.compute_time_derivative(&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..89ed4c6 100644 --- a/math_explorer/src/physics/fluid_dynamics/tests.rs +++ b/math_explorer/src/physics/fluid_dynamics/tests.rs @@ -4,9 +4,10 @@ mod tests { analysis::{bernoulli_constant, reynolds_number, shear_stress}, conservation::{ continuity_divergence, material_derivative_scalar, navier_stokes_time_derivative, + MomentumEquation, NavierStokes, Euler, }, regimes::{FlatPlateClassifier, FlowClassifier, FlowRegime, PipeFlowClassifier}, - types::{FlowState, FluidProperties}, + types::{FlowState, FluidProperties, SpatialGradients}, }; use nalgebra::{Matrix3, Vector3}; @@ -113,4 +114,26 @@ 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); // rho=1, mu=1 + let state = FlowState::new(Vector3::new(1.0, 0.0, 0.0), 0.0); + + let vel_grad = Matrix3::zeros(); + let p_grad = Vector3::zeros(); + let lap_vel = Vector3::new(1.0, 0.0, 0.0); + let g = Vector3::zeros(); + + // Check Navier Stokes: should include viscosity + // Term = nu * lap_vel = 1.0 * 1.0 = 1.0 + let gradients_ns = SpatialGradients::new(vel_grad, p_grad, Some(lap_vel)); + let ns_accel = NavierStokes.compute_time_derivative(&props, &state, &gradients_ns, g); + assert!((ns_accel.x - 1.0).abs() < 1e-9); + + // Check Euler: should ignore viscosity even if laplacian is provided + let gradients_euler_with_lap = SpatialGradients::new(vel_grad, p_grad, Some(lap_vel)); + let euler_accel_2 = Euler.compute_time_derivative(&props, &state, &gradients_euler_with_lap, g); + assert!((euler_accel_2.x - 0.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..e1400b6 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,28 @@ impl FlowState { Self { velocity, pressure } } } + +/// Encapsulates spatial gradients required for momentum equations. +#[derive(Debug, Clone)] +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>, +} + +impl SpatialGradients { + pub fn new( + velocity_gradient: Matrix3, + pressure_gradient: Vector3, + laplacian_velocity: Option>, + ) -> Self { + Self { + velocity_gradient, + pressure_gradient, + laplacian_velocity, + } + } +} From 1a536e87db99ae158c079c60c2a1556f1bc99548 Mon Sep 17 00:00:00 2001 From: fderuiter <127706008+fderuiter@users.noreply.github.com> Date: Thu, 29 Jan 2026 19:32:25 +0000 Subject: [PATCH 2/2] Fix formatting in fluid dynamics module - Ran `cargo fmt` to address CI failure in `conservation.rs` and `tests.rs`. - Ensured compliance with project style guidelines. Co-authored-by: google-labs-jules[bot] <161369871+google-labs-jules[bot]@users.noreply.github.com> --- math_explorer/src/physics/fluid_dynamics/conservation.rs | 6 +----- math_explorer/src/physics/fluid_dynamics/tests.rs | 7 ++++--- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/math_explorer/src/physics/fluid_dynamics/conservation.rs b/math_explorer/src/physics/fluid_dynamics/conservation.rs index 0d6e265..cb77217 100644 --- a/math_explorer/src/physics/fluid_dynamics/conservation.rs +++ b/math_explorer/src/physics/fluid_dynamics/conservation.rs @@ -165,10 +165,6 @@ pub fn euler_time_derivative( ) -> Vector3 { // Construct dummy properties with the given rho let properties = FluidProperties::new(rho, 0.0); // Viscosity ignored by Euler - let gradients = SpatialGradients::new( - *velocity_gradient, - pressure_gradient, - None, - ); + let gradients = SpatialGradients::new(*velocity_gradient, pressure_gradient, None); Euler.compute_time_derivative(&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 89ed4c6..df781af 100644 --- a/math_explorer/src/physics/fluid_dynamics/tests.rs +++ b/math_explorer/src/physics/fluid_dynamics/tests.rs @@ -3,8 +3,8 @@ mod tests { use crate::physics::fluid_dynamics::{ analysis::{bernoulli_constant, reynolds_number, shear_stress}, conservation::{ - continuity_divergence, material_derivative_scalar, navier_stokes_time_derivative, - MomentumEquation, NavierStokes, Euler, + Euler, MomentumEquation, NavierStokes, continuity_divergence, + material_derivative_scalar, navier_stokes_time_derivative, }, regimes::{FlatPlateClassifier, FlowClassifier, FlowRegime, PipeFlowClassifier}, types::{FlowState, FluidProperties, SpatialGradients}, @@ -133,7 +133,8 @@ mod tests { // Check Euler: should ignore viscosity even if laplacian is provided let gradients_euler_with_lap = SpatialGradients::new(vel_grad, p_grad, Some(lap_vel)); - let euler_accel_2 = Euler.compute_time_derivative(&props, &state, &gradients_euler_with_lap, g); + let euler_accel_2 = + Euler.compute_time_derivative(&props, &state, &gradients_euler_with_lap, g); assert!((euler_accel_2.x - 0.0).abs() < 1e-9); } }