qualia_core_db/specialized_libs/physics_simulation/
discretization.rs

1use super::*;
2
3/// Spatial discretizer
4pub struct SpatialDiscretizer {
5    discretization_method: SpatialDiscretizationMethod,
6    grid_generator: GridGenerator,
7    mesh_generator: MeshGenerator,
8    stencil_operators: StencilOperators,
9}
10
11/// Spatial discretization methods
12#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
13pub enum SpatialDiscretizationMethod {
14    /// Structured grid
15    Structured,
16    /// Unstructured grid
17    Unstructured,
18    /// Adaptive mesh refinement
19    AdaptiveMeshRefinement,
20    /// Moving mesh
21    MovingMesh,
22    /// Spectral element
23    SpectralElement,
24    /// Discontinuous Galerkin
25    DiscontinuousGalerkin,
26}
27
28/// Grid generator
29pub struct GridGenerator {
30    grid_type: GridType,
31    grid_parameters: GridParameters,
32    quality_metrics: GridQualityMetrics,
33}
34
35/// Grid types
36#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
37pub enum GridType {
38    /// Cartesian grid
39    Cartesian,
40    /// Curvilinear grid
41    Curvilinear,
42    /// Body-fitted grid
43    BodyFitted,
44    /// Overset grid
45    Overset,
46    /// Chimera grid
47    Chimera,
48    /// Adaptive grid
49    Adaptive,
50}
51
52/// Grid parameters
53#[derive(Debug, Clone, Serialize, Deserialize)]
54pub struct GridParameters {
55    pub domain_bounds: Vec<(f64, f64)>,
56    pub grid_spacing: Vec<f64>,
57    pub stretching_function: Option<String>,
58    pub boundary_layer: Option<BoundaryLayerConfig>,
59}
60
61/// Boundary layer configuration
62#[derive(Debug, Clone, Serialize, Deserialize)]
63pub struct BoundaryLayerConfig {
64    pub thickness: f64,
65    pub stretching_ratio: f64,
66    pub num_points: usize,
67}
68
69/// Grid quality metrics
70#[derive(Debug, Clone)]
71pub struct GridQualityMetrics {
72    pub orthogonality: f64,
73    pub skewness: f64,
74    pub aspect_ratio: f64,
75    pub smoothness: f64,
76    pub expansion_ratio: f64,
77}
78
79/// Mesh generator
80pub struct MeshGenerator {
81    mesh_type: MeshType,
82    mesh_parameters: MeshParameters,
83    quality_metrics: MeshQualityMetrics,
84}
85
86/// Mesh types
87#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
88pub enum MeshType {
89    /// Triangular mesh
90    Triangular,
91    /// Quadrilateral mesh
92    Quadrilateral,
93    /// Tetrahedral mesh
94    Tetrahedral,
95    /// Hexahedral mesh
96    Hexahedral,
97    /// Mixed mesh
98    Mixed,
99    /// Hybrid mesh
100    Hybrid,
101}
102
103/// Mesh parameters
104#[derive(Debug, Clone, Serialize, Deserialize)]
105pub struct MeshParameters {
106    pub element_size: f64,
107    pub grading_factor: f64,
108    pub refinement_regions: Vec<RefinementRegion>,
109    pub boundary_layer: Option<BoundaryLayerConfig>,
110}
111
112/// Refinement regions
113#[derive(Debug, Clone, Serialize, Deserialize)]
114pub struct RefinementRegion {
115    pub region_bounds: Vec<(f64, f64)>,
116    pub refinement_factor: f64,
117    pub element_size: f64,
118}
119
120/// Mesh quality metrics
121#[derive(Debug, Clone)]
122pub struct MeshQualityMetrics {
123    pub element_quality: f64,
124    pub node_distribution: f64,
125    pub connectivity: f64,
126    pub aspect_ratio: f64,
127}
128
129/// Stencil operators
130pub struct StencilOperators {
131    operators: HashMap<String, StencilOperator>,
132    boundary_stencils: HashMap<String, BoundaryStencil>,
133}
134
135/// Stencil operator
136#[derive(Debug, Clone)]
137pub struct StencilOperator {
138    pub operator_id: String,
139    pub operator_type: StencilType,
140    pub stencil_points: Vec<StencilPoint>,
141    pub coefficients: Vec<f64>,
142}
143
144/// Stencil types
145#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
146pub enum StencilType {
147    /// Central difference
148    Central,
149    /// Forward difference
150    Forward,
151    /// Backward difference
152    Backward,
153    /// Upwind
154    Upwind,
155    /// High-order compact
156    HighOrderCompact,
157    /// WENO scheme
158    WENO,
159    /// ENO scheme
160    ENO,
161}
162
163/// Stencil point
164#[derive(Debug, Clone)]
165pub struct StencilPoint {
166    pub relative_position: Vec<i32>,
167    pub weight: f64,
168}
169
170/// Boundary stencil
171#[derive(Debug, Clone)]
172pub struct BoundaryStencil {
173    pub stencil_id: String,
174    pub boundary_type: BoundaryType,
175    pub stencil_points: Vec<StencilPoint>,
176    pub coefficients: Vec<f64>,
177}
178
179impl SpatialDiscretizer {
180    pub fn new() -> Self {
181        Self {
182            discretization_method: SpatialDiscretizationMethod::Structured,
183            grid_generator: GridGenerator::new(),
184            mesh_generator: MeshGenerator::new(),
185            stencil_operators: StencilOperators::new(),
186        }
187    }
188
189    pub fn initialize(&mut self) -> Result<(), PhysicsError> {
190        self.grid_generator.initialize()?;
191        self.mesh_generator.initialize()?;
192        Ok(())
193    }
194
195    /// Get the discretization method.
196    pub fn get_discretization_method(&self) -> &SpatialDiscretizationMethod {
197        &self.discretization_method
198    }
199
200    /// Set the discretization method.
201    pub fn set_discretization_method(&mut self, method: SpatialDiscretizationMethod) {
202        self.discretization_method = method;
203    }
204
205    /// Get a reference to the stencil operators.
206    pub fn get_stencil_operators(&self) -> &StencilOperators {
207        &self.stencil_operators
208    }
209
210    /// Get a mutable reference to the stencil operators.
211    pub fn get_stencil_operators_mut(&mut self) -> &mut StencilOperators {
212        &mut self.stencil_operators
213    }
214}
215
216impl GridGenerator {
217    pub fn new() -> Self {
218        Self {
219            grid_type: GridType::Cartesian,
220            grid_parameters: GridParameters::new(),
221            quality_metrics: GridQualityMetrics::new(),
222        }
223    }
224
225    pub fn initialize(&mut self) -> Result<(), PhysicsError> {
226        Ok(())
227    }
228
229    /// Get the grid type.
230    pub fn get_grid_type(&self) -> &GridType {
231        &self.grid_type
232    }
233
234    /// Set the grid type.
235    pub fn set_grid_type(&mut self, grid_type: GridType) {
236        self.grid_type = grid_type;
237    }
238
239    /// Get a reference to the grid parameters.
240    pub fn get_grid_parameters(&self) -> &GridParameters {
241        &self.grid_parameters
242    }
243
244    /// Get a mutable reference to the grid parameters.
245    pub fn get_grid_parameters_mut(&mut self) -> &mut GridParameters {
246        &mut self.grid_parameters
247    }
248
249    /// Get a reference to the grid quality metrics.
250    pub fn get_quality_metrics(&self) -> &GridQualityMetrics {
251        &self.quality_metrics
252    }
253
254    /// Get a mutable reference to the grid quality metrics.
255    pub fn get_quality_metrics_mut(&mut self) -> &mut GridQualityMetrics {
256        &mut self.quality_metrics
257    }
258}
259
260impl GridParameters {
261    pub fn new() -> Self {
262        Self {
263            domain_bounds: vec![(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)],
264            grid_spacing: vec![0.01, 0.01, 0.01],
265            stretching_function: None,
266            boundary_layer: None,
267        }
268    }
269}
270
271impl GridQualityMetrics {
272    pub fn new() -> Self {
273        Self {
274            orthogonality: 1.0,
275            skewness: 0.0,
276            aspect_ratio: 1.0,
277            smoothness: 1.0,
278            expansion_ratio: 1.0,
279        }
280    }
281}
282
283impl MeshGenerator {
284    pub fn new() -> Self {
285        Self {
286            mesh_type: MeshType::Hexahedral,
287            mesh_parameters: MeshParameters::new(),
288            quality_metrics: MeshQualityMetrics::new(),
289        }
290    }
291
292    pub fn initialize(&mut self) -> Result<(), PhysicsError> {
293        Ok(())
294    }
295
296    /// Get the mesh type.
297    pub fn get_mesh_type(&self) -> &MeshType {
298        &self.mesh_type
299    }
300
301    /// Set the mesh type.
302    pub fn set_mesh_type(&mut self, mesh_type: MeshType) {
303        self.mesh_type = mesh_type;
304    }
305
306    /// Get a reference to the mesh parameters.
307    pub fn get_mesh_parameters(&self) -> &MeshParameters {
308        &self.mesh_parameters
309    }
310
311    /// Get a mutable reference to the mesh parameters.
312    pub fn get_mesh_parameters_mut(&mut self) -> &mut MeshParameters {
313        &mut self.mesh_parameters
314    }
315
316    /// Get a reference to the mesh quality metrics.
317    pub fn get_quality_metrics(&self) -> &MeshQualityMetrics {
318        &self.quality_metrics
319    }
320
321    /// Get a mutable reference to the mesh quality metrics.
322    pub fn get_quality_metrics_mut(&mut self) -> &mut MeshQualityMetrics {
323        &mut self.quality_metrics
324    }
325}
326
327impl MeshParameters {
328    pub fn new() -> Self {
329        Self {
330            element_size: 0.01,
331            grading_factor: 1.2,
332            refinement_regions: Vec::new(),
333            boundary_layer: None,
334        }
335    }
336}
337
338impl MeshQualityMetrics {
339    pub fn new() -> Self {
340        Self {
341            element_quality: 1.0,
342            node_distribution: 1.0,
343            connectivity: 1.0,
344            aspect_ratio: 1.0,
345        }
346    }
347}
348
349impl StencilOperators {
350    pub fn new() -> Self {
351        Self {
352            operators: HashMap::new(),
353            boundary_stencils: HashMap::new(),
354        }
355    }
356
357    /// Register a named stencil operator.
358    pub fn register_operator(&mut self, name: &str, stencil: StencilOperator) {
359        self.operators.insert(name.to_string(), stencil);
360    }
361
362    /// Register a named boundary stencil.
363    pub fn register_boundary_stencil(&mut self, name: &str, stencil: BoundaryStencil) {
364        self.boundary_stencils.insert(name.to_string(), stencil);
365    }
366
367    /// Apply a registered stencil to compute the spatial derivative at `index`.
368    ///
369    /// The derivative is computed as:
370    /// ```text
371    ///   sum_i( coefficients[i] * field[index + offset_i] ) / dx
372    /// ```
373    /// where `offset_i` is taken from `stencil_points[i].relative_position[0]`.
374    /// The coefficients are expected to already include the normalisation factor
375    /// (e.g. `[-0.5, 0.0, 0.5]` for a 2nd-order central difference).
376    pub fn apply_derivative(
377        &self,
378        name: &str,
379        field: &[f64],
380        dx: f64,
381        index: usize,
382    ) -> Result<f64, PhysicsError> {
383        let stencil = self.operators.get(name).ok_or_else(|| {
384            PhysicsError::SolverError(format!("Stencil operator '{}' not registered", name))
385        })?;
386
387        if stencil.stencil_points.len() != stencil.coefficients.len() {
388            return Err(PhysicsError::SolverError(format!(
389                "Stencil operator '{}' has mismatched points/coefficients",
390                name
391            )));
392        }
393
394        let n = field.len() as isize;
395        let mut sum = 0.0f64;
396        for (point, coeff) in stencil
397            .stencil_points
398            .iter()
399            .zip(stencil.coefficients.iter())
400        {
401            let offset = point.relative_position.first().copied().unwrap_or(0) as isize;
402            let idx = index as isize + offset;
403            if idx < 0 || idx >= n {
404                return Err(PhysicsError::SolverError(format!(
405                    "Stencil operator '{}' accesses out-of-bounds index {} (field len {})",
406                    name, idx, n
407                )));
408            }
409            sum += coeff * field[idx as usize];
410        }
411
412        if dx <= 0.0 {
413            return Err(PhysicsError::SolverError("dx must be positive".to_string()));
414        }
415
416        Ok(sum / dx)
417    }
418
419    /// Create a 3-point 2nd-order central difference stencil.
420    ///
421    /// Coefficients `[-0.5, 0.0, 0.5]` at offsets `[-1, 0, +1]` give
422    /// `(field[i+1] - field[i-1]) / (2*dx)`.
423    pub fn central_difference_2nd_order() -> StencilOperator {
424        StencilOperator {
425            operator_id: "central_difference_2nd_order".to_string(),
426            operator_type: StencilType::Central,
427            stencil_points: vec![
428                StencilPoint {
429                    relative_position: vec![-1],
430                    weight: 1.0,
431                },
432                StencilPoint {
433                    relative_position: vec![0],
434                    weight: 1.0,
435                },
436                StencilPoint {
437                    relative_position: vec![1],
438                    weight: 1.0,
439                },
440            ],
441            coefficients: vec![-0.5, 0.0, 0.5],
442        }
443    }
444
445    /// Create a 2-point 1st-order forward difference stencil.
446    ///
447    /// Coefficients `[-1.0, 1.0]` at offsets `[0, +1]` give
448    /// `(field[i+1] - field[i]) / dx`.
449    pub fn forward_difference_1st_order() -> StencilOperator {
450        StencilOperator {
451            operator_id: "forward_difference_1st_order".to_string(),
452            operator_type: StencilType::Forward,
453            stencil_points: vec![
454                StencilPoint {
455                    relative_position: vec![0],
456                    weight: 1.0,
457                },
458                StencilPoint {
459                    relative_position: vec![1],
460                    weight: 1.0,
461                },
462            ],
463            coefficients: vec![-1.0, 1.0],
464        }
465    }
466}
qualia_core_db/specialized_libs/physics_simulation/discretization.rs

qualia_core_db/specialized_libs/physics_simulation/
discretization.rs
