qualia_core_db/specialized_libs/physics_simulation/
solvers.rs

1use super::*;
2
3/// Physics solver
4pub struct PhysicsSolver {
5    solver_type: SolverType,
6    linear_solver: LinearSolver,
7    nonlinear_solver: NonlinearSolver,
8    eigenvalue_solver: EigenvalueSolver,
9    optimization_solver: OptimizationSolver,
10}
11
12/// Solver types
13#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
14pub enum SolverType {
15    /// Direct solver
16    Direct,
17    /// Iterative solver
18    Iterative,
19    /// Multigrid solver
20    Multigrid,
21    /// Domain decomposition solver
22    DomainDecomposition,
23    /// Hybrid solver
24    Hybrid,
25}
26
27/// CFD (Computational Fluid Dynamics) solver
28pub struct CfdSolver {
29    solver_id: String,
30    solver_method: LinearSolverMethod,
31    preconditioner: Preconditioner,
32    convergence_criteria: ConvergenceCriteria,
33    solver_parameters: SolverParameters,
34}
35
36/// Solver result for physics computations
37pub struct SolverResult {
38    pub solver_id: String,
39    pub iterations: u64,
40    pub residual_norm: f64,
41    pub convergence_time: f64,
42    pub error_message: Option<String>,
43}
44
45/// Distribution of simulation work across mesh nodes
46pub struct NodeDistribution {
47    pub node_ids: Vec<String>,
48    pub node_loads: Vec<f64>,
49    pub communication_pattern: CommunicationPattern,
50}
51
52/// Linear solver
53pub struct LinearSolver {
54    solver_method: LinearSolverMethod,
55    preconditioner: Preconditioner,
56    convergence_criteria: ConvergenceCriteria,
57    solver_parameters: SolverParameters,
58}
59
60/// Linear solver methods
61#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
62pub enum LinearSolverMethod {
63    /// Gaussian elimination
64    GaussianElimination,
65    /// LU decomposition
66    LUDecomposition,
67    /// Cholesky decomposition
68    CholeskyDecomposition,
69    /// QR decomposition
70    QRDecomposition,
71    /// Conjugate gradient method
72    ConjugateGradient,
73    /// GMRES method
74    GMRES,
75    /// BiCGSTAB method
76    BiCGSTAB,
77    /// Multigrid method
78    Multigrid,
79}
80
81/// Preconditioner
82#[derive(Debug, Clone)]
83pub struct Preconditioner {
84    preconditioner_type: PreconditionerType,
85    preconditioner_parameters: PreconditionerParameters,
86}
87
88/// Preconditioner types
89#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
90pub enum PreconditionerType {
91    /// Jacobi preconditioner
92    Jacobi,
93    /// Gauss-Seidel preconditioner
94    GaussSeidel,
95    /// Successive over-relaxation (SOR)
96    SOR,
97    /// Incomplete LU (ILU)
98    ILU,
99    /// Algebraic multigrid (AMG)
100    AMG,
101    /// Block preconditioner
102    Block,
103}
104
105/// Preconditioner parameters
106#[derive(Debug, Clone, Serialize, Deserialize)]
107pub struct PreconditionerParameters {
108    pub relaxation_factor: f64,
109    pub fill_level: usize,
110    pub tolerance: f64,
111    pub max_iterations: usize,
112}
113
114/// Convergence criteria
115#[derive(Debug, Clone)]
116pub struct ConvergenceCriteria {
117    pub tolerance: f64,
118    pub max_iterations: usize,
119    pub relative_tolerance: f64,
120    pub absolute_tolerance: f64,
121    pub divergence_check: bool,
122}
123
124/// Solver parameters
125#[derive(Debug, Clone, Serialize, Deserialize)]
126pub struct SolverParameters {
127    pub tolerance: f64,
128    pub max_iterations: usize,
129    pub restart_frequency: usize,
130    pub orthogonalization: OrthogonalizationMethod,
131}
132
133/// Orthogonalization methods
134#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
135pub enum OrthogonalizationMethod {
136    /// Classical Gram-Schmidt
137    ClassicalGramSchmidt,
138    /// Modified Gram-Schmidt
139    ModifiedGramSchmidt,
140    /// Householder
141    Householder,
142    /// Givens rotations
143    Givens,
144}
145
146/// Nonlinear solver
147pub struct NonlinearSolver {
148    solver_method: NonlinearSolverMethod,
149    linear_solver: LinearSolver,
150    convergence_criteria: ConvergenceCriteria,
151    solver_parameters: NonlinearSolverParameters,
152}
153
154/// Nonlinear solver methods
155#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
156pub enum NonlinearSolverMethod {
157    /// Newton-Raphson method
158    NewtonRaphson,
159    /// Quasi-Newton method
160    QuasiNewton,
161    /// Fixed-point iteration
162    FixedPoint,
163    /// Picard iteration
164    Picard,
165    /// Anderson acceleration
166    Anderson,
167    /// Broyden's method
168    Broyden,
169}
170
171/// Nonlinear solver parameters
172#[derive(Debug, Clone, Serialize, Deserialize)]
173pub struct NonlinearSolverParameters {
174    pub tolerance: f64,
175    pub max_iterations: usize,
176    pub line_search: LineSearchMethod,
177    pub trust_region: TrustRegionMethod,
178}
179
180/// Line search methods
181#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
182pub enum LineSearchMethod {
183    /// Backtracking line search
184    Backtracking,
185    /// Wolfe conditions
186    Wolfe,
187    /// Goldstein conditions
188    Goldstein,
189    /// Armijo rule
190    Armijo,
191}
192
193/// Trust region methods
194#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
195pub enum TrustRegionMethod {
196    /// Dogleg method
197    Dogleg,
198    /// Double dogleg method
199    DoubleDogleg,
200    /// Powell method
201    Powell,
202    /// Levenberg-Marquardt
203    LevenbergMarquardt,
204}
205
206impl PhysicsSolver {
207    pub fn new() -> Self {
208        Self {
209            solver_type: SolverType::Iterative,
210            linear_solver: LinearSolver::new(),
211            nonlinear_solver: NonlinearSolver::new(),
212            eigenvalue_solver: EigenvalueSolver::new(),
213            optimization_solver: OptimizationSolver::new(),
214        }
215    }
216
217    pub fn initialize(&mut self) -> Result<(), PhysicsError> {
218        self.linear_solver.initialize()?;
219        self.nonlinear_solver.initialize()?;
220        self.eigenvalue_solver.initialize()?;
221        self.optimization_solver.initialize()?;
222        Ok(())
223    }
224
225    pub fn create_cfd_solver(&self, _config: &SimulationConfig) -> Result<CfdSolver, PhysicsError> {
226        let solver = CfdSolver {
227            solver_id: "cfd_solver".to_string(),
228            solver_method: LinearSolverMethod::GMRES,
229            preconditioner: Preconditioner::new(),
230            convergence_criteria: ConvergenceCriteria::new(),
231            solver_parameters: SolverParameters::new(),
232        };
233
234        Ok(solver)
235    }
236
237    pub fn solve_cfd_step(
238        &self,
239        _solver: &CfdSolver,
240        fields: &[PhysicsField],
241        _mesh: &Mesh,
242    ) -> Result<SolverResult, PhysicsError> {
243        // Real steady-state residual of the velocity field: the L2 norm of the Burgers
244        // operator ‖ν·u_xx − u·u_x‖ over the interior nodes — a genuine measure of how far
245        // the field is from a steady solution. (Previously this returned a fabricated 1e-7.)
246        let start = std::time::Instant::now();
247        let velocity = fields
248            .iter()
249            .find(|f| f.metadata.physical_quantity == "Velocity");
250        let (iterations, residual_norm) = match velocity {
251            Some(v) if v.data.len() >= 3 => {
252                let u = &v.data;
253                let n = u.len();
254                let dx = 1.0 / n as f64;
255                let nu = 1.5e-5_f64;
256                let mut sumsq = 0.0f64;
257                for i in 1..n - 1 {
258                    let u_x = (u[i + 1] - u[i - 1]) / (2.0 * dx);
259                    let u_xx = (u[i + 1] - 2.0 * u[i] + u[i - 1]) / (dx * dx);
260                    let r = nu * u_xx - u[i] * u_x;
261                    sumsq += r * r;
262                }
263                (1u64, sumsq.sqrt())
264            }
265            _ => (0u64, f64::MAX),
266        };
267
268        Ok(SolverResult {
269            solver_id: "cfd_solver".to_string(),
270            iterations,
271            residual_norm,
272            convergence_time: start.elapsed().as_secs_f64(),
273            error_message: None,
274        })
275    }
276
277    /// Get the solver type.
278    pub fn get_solver_type(&self) -> &SolverType {
279        &self.solver_type
280    }
281
282    /// Set the solver type.
283    pub fn set_solver_type(&mut self, solver_type: SolverType) {
284        self.solver_type = solver_type;
285    }
286}
287
288impl LinearSolver {
289    pub fn new() -> Self {
290        Self {
291            solver_method: LinearSolverMethod::GMRES,
292            preconditioner: Preconditioner::new(),
293            convergence_criteria: ConvergenceCriteria::new(),
294            solver_parameters: SolverParameters::new(),
295        }
296    }
297
298    pub fn initialize(&mut self) -> Result<(), PhysicsError> {
299        Ok(())
300    }
301
302    /// Get the solver method.
303    pub fn get_solver_method(&self) -> &LinearSolverMethod {
304        &self.solver_method
305    }
306
307    /// Set the solver method.
308    pub fn set_solver_method(&mut self, method: LinearSolverMethod) {
309        self.solver_method = method;
310    }
311
312    /// Get a reference to the preconditioner.
313    pub fn get_preconditioner(&self) -> &Preconditioner {
314        &self.preconditioner
315    }
316
317    /// Get a mutable reference to the preconditioner.
318    pub fn get_preconditioner_mut(&mut self) -> &mut Preconditioner {
319        &mut self.preconditioner
320    }
321
322    /// Get a reference to the convergence criteria.
323    pub fn get_convergence_criteria(&self) -> &ConvergenceCriteria {
324        &self.convergence_criteria
325    }
326
327    /// Get a mutable reference to the convergence criteria.
328    pub fn get_convergence_criteria_mut(&mut self) -> &mut ConvergenceCriteria {
329        &mut self.convergence_criteria
330    }
331
332    /// Get a reference to the solver parameters.
333    pub fn get_solver_parameters(&self) -> &SolverParameters {
334        &self.solver_parameters
335    }
336
337    /// Get a mutable reference to the solver parameters.
338    pub fn get_solver_parameters_mut(&mut self) -> &mut SolverParameters {
339        &mut self.solver_parameters
340    }
341}
342
343impl CfdSolver {
344    /// Get the solver ID.
345    pub fn get_solver_id(&self) -> &str {
346        &self.solver_id
347    }
348
349    /// Get the solver method.
350    pub fn get_solver_method(&self) -> &LinearSolverMethod {
351        &self.solver_method
352    }
353
354    /// Set the solver method.
355    pub fn set_solver_method(&mut self, method: LinearSolverMethod) {
356        self.solver_method = method;
357    }
358
359    /// Get a reference to the preconditioner.
360    pub fn get_preconditioner(&self) -> &Preconditioner {
361        &self.preconditioner
362    }
363
364    /// Get a mutable reference to the preconditioner.
365    pub fn get_preconditioner_mut(&mut self) -> &mut Preconditioner {
366        &mut self.preconditioner
367    }
368
369    /// Get a reference to the convergence criteria.
370    pub fn get_convergence_criteria(&self) -> &ConvergenceCriteria {
371        &self.convergence_criteria
372    }
373
374    /// Get a mutable reference to the convergence criteria.
375    pub fn get_convergence_criteria_mut(&mut self) -> &mut ConvergenceCriteria {
376        &mut self.convergence_criteria
377    }
378
379    /// Get a reference to the solver parameters.
380    pub fn get_solver_parameters(&self) -> &SolverParameters {
381        &self.solver_parameters
382    }
383
384    /// Get a mutable reference to the solver parameters.
385    pub fn get_solver_parameters_mut(&mut self) -> &mut SolverParameters {
386        &mut self.solver_parameters
387    }
388}
389
390impl Preconditioner {
391    pub fn new() -> Self {
392        Self {
393            preconditioner_type: PreconditionerType::ILU,
394            preconditioner_parameters: PreconditionerParameters::new(),
395        }
396    }
397
398    /// Get the preconditioner type.
399    pub fn get_preconditioner_type(&self) -> &PreconditionerType {
400        &self.preconditioner_type
401    }
402
403    /// Set the preconditioner type.
404    pub fn set_preconditioner_type(&mut self, ptype: PreconditionerType) {
405        self.preconditioner_type = ptype;
406    }
407
408    /// Get a reference to the preconditioner parameters.
409    pub fn get_preconditioner_parameters(&self) -> &PreconditionerParameters {
410        &self.preconditioner_parameters
411    }
412
413    /// Get a mutable reference to the preconditioner parameters.
414    pub fn get_preconditioner_parameters_mut(&mut self) -> &mut PreconditionerParameters {
415        &mut self.preconditioner_parameters
416    }
417}
418
419impl PreconditionerParameters {
420    pub fn new() -> Self {
421        Self {
422            relaxation_factor: 1.0,
423            fill_level: 0,
424            tolerance: 1e-6,
425            max_iterations: 100,
426        }
427    }
428}
429
430impl ConvergenceCriteria {
431    pub fn new() -> Self {
432        Self {
433            tolerance: 1e-6,
434            max_iterations: 1000,
435            relative_tolerance: 1e-6,
436            absolute_tolerance: 1e-12,
437            divergence_check: true,
438        }
439    }
440}
441
442impl SolverParameters {
443    pub fn new() -> Self {
444        Self {
445            tolerance: 1e-6,
446            max_iterations: 1000,
447            restart_frequency: 100,
448            orthogonalization: OrthogonalizationMethod::ModifiedGramSchmidt,
449        }
450    }
451}
452
453impl NonlinearSolver {
454    pub fn new() -> Self {
455        Self {
456            solver_method: NonlinearSolverMethod::NewtonRaphson,
457            linear_solver: LinearSolver::new(),
458            convergence_criteria: ConvergenceCriteria::new(),
459            solver_parameters: NonlinearSolverParameters::new(),
460        }
461    }
462
463    pub fn initialize(&mut self) -> Result<(), PhysicsError> {
464        self.linear_solver.initialize()?;
465        Ok(())
466    }
467
468    /// Get the solver method.
469    pub fn get_solver_method(&self) -> &NonlinearSolverMethod {
470        &self.solver_method
471    }
472
473    /// Set the solver method.
474    pub fn set_solver_method(&mut self, method: NonlinearSolverMethod) {
475        self.solver_method = method;
476    }
477
478    /// Get a reference to the convergence criteria.
479    pub fn get_convergence_criteria(&self) -> &ConvergenceCriteria {
480        &self.convergence_criteria
481    }
482
483    /// Get a mutable reference to the convergence criteria.
484    pub fn get_convergence_criteria_mut(&mut self) -> &mut ConvergenceCriteria {
485        &mut self.convergence_criteria
486    }
487
488    /// Get a reference to the solver parameters.
489    pub fn get_solver_parameters(&self) -> &NonlinearSolverParameters {
490        &self.solver_parameters
491    }
492
493    /// Get a mutable reference to the solver parameters.
494    pub fn get_solver_parameters_mut(&mut self) -> &mut NonlinearSolverParameters {
495        &mut self.solver_parameters
496    }
497}
498
499impl NonlinearSolverParameters {
500    pub fn new() -> Self {
501        Self {
502            tolerance: 1e-6,
503            max_iterations: 100,
504            line_search: LineSearchMethod::Backtracking,
505            trust_region: TrustRegionMethod::LevenbergMarquardt,
506        }
507    }
508}
qualia_core_db/specialized_libs/physics_simulation/solvers.rs

qualia_core_db/specialized_libs/physics_simulation/
solvers.rs
