scirs2_integrate/specialized/fluid_dynamics/core/
parameters.rs

1//! Parameter structures for fluid dynamics solvers.
2
3/// Navier-Stokes solver parameters
4#[derive(Debug, Clone)]
5pub struct NavierStokesParams {
6    /// Kinematic viscosity
7    pub nu: f64,
8    /// Density (for incompressible flow, usually 1.0)
9    pub rho: f64,
10    /// Time step
11    pub dt: f64,
12    /// Maximum iterations for pressure solver
13    pub max_pressure_iter: usize,
14    /// Tolerance for pressure solver
15    pub pressure_tol: f64,
16    /// Use semi-Lagrangian advection
17    pub semi_lagrangian: bool,
18}
19
20impl Default for NavierStokesParams {
21    fn default() -> Self {
22        Self {
23            nu: 0.01,
24            rho: 1.0,
25            dt: 0.01,
26            max_pressure_iter: 100,
27            pressure_tol: 1e-6,
28            semi_lagrangian: false,
29        }
30    }
31}
32
33impl NavierStokesParams {
34    /// Create new parameters with specified Reynolds number
35    /// Reynolds number = U * L / nu, where U is characteristic velocity, L is characteristic length
36    pub fn with_reynolds_number(
37        reynolds: f64,
38        characteristic_velocity: f64,
39        characteristic_length: f64,
40    ) -> Self {
41        let nu = characteristic_velocity * characteristic_length / reynolds;
42        Self {
43            nu,
44            ..Default::default()
45        }
46    }
47
48    /// Create new parameters with specified viscosity
49    pub fn with_viscosity(nu: f64) -> Self {
50        Self {
51            nu,
52            ..Default::default()
53        }
54    }
55
56    /// Create new parameters with specified time step
57    pub fn with_time_step(dt: f64) -> Self {
58        Self {
59            dt,
60            ..Default::default()
61        }
62    }
63
64    /// Calculate the Reynolds number given characteristic velocity and length
65    pub fn reynolds_number(&self, characteristic_velocity: f64, characteristiclength: f64) -> f64 {
66        characteristic_velocity * characteristiclength / self.nu
67    }
68
69    /// Calculate CFL number for stability analysis
70    /// CFL = u * dt / dx + v * dt / dy (+ w * dt / dz for 3D)
71    pub fn cfl_number_2d(&self, u_max: f64, vmax: f64, dx: f64, dy: f64) -> f64 {
72        u_max * self.dt / dx + vmax * self.dt / dy
73    }
74
75    /// Calculate CFL number for 3D case
76    pub fn cfl_number_3d(
77        &self,
78        u_max: f64,
79        v_max: f64,
80        w_max: f64,
81        dx: f64,
82        dy: f64,
83        dz: f64,
84    ) -> f64 {
85        u_max * self.dt / dx + v_max * self.dt / dy + w_max * self.dt / dz
86    }
87
88    /// Calculate viscous CFL number
89    /// Viscous CFL = nu * dt * (1/dx^2 + 1/dy^2 + 1/dz^2)
90    pub fn viscous_cfl_2d(&self, dx: f64, dy: f64) -> f64 {
91        self.nu * self.dt * (1.0 / (dx * dx) + 1.0 / (dy * dy))
92    }
93
94    /// Calculate viscous CFL number for 3D case
95    pub fn viscous_cfl_3d(&self, dx: f64, dy: f64, dz: f64) -> f64 {
96        self.nu * self.dt * (1.0 / (dx * dx) + 1.0 / (dy * dy) + 1.0 / (dz * dz))
97    }
98
99    /// Validate parameters for numerical stability
100    pub fn validate(&self) -> Result<(), String> {
101        if self.nu <= 0.0 {
102            return Err("Kinematic viscosity must be positive".to_string());
103        }
104        if self.rho <= 0.0 {
105            return Err("Density must be positive".to_string());
106        }
107        if self.dt <= 0.0 {
108            return Err("Time step must be positive".to_string());
109        }
110        if self.max_pressure_iter == 0 {
111            return Err("Maximum pressure iterations must be greater than zero".to_string());
112        }
113        if self.pressure_tol <= 0.0 {
114            return Err("Pressure tolerance must be positive".to_string());
115        }
116        Ok(())
117    }
118}