bem/core/solver/
fmm_interface.rs

1//! FMM-solver interface
2//!
3//! This module provides a unified interface for solving BEM systems using
4//! either direct methods (TBEM) or iterative methods with FMM acceleration.
5//!
6//! The key abstraction is the `LinearOperator` trait, which allows iterative
7//! solvers to work with any matrix representation.
8//!
9//! ## Preconditioning
10//!
11//! BEM systems are typically ill-conditioned. Simple diagonal (Jacobi) preconditioning
12//! is **not sufficient** for BEM. Use ILU preconditioning instead:
13//!
14//! ```ignore
15//! use bem::core::solver::{IluMethod, IluScanningDegree, solve_with_ilu};
16//!
17//! // For TBEM systems
18//! let solution = solve_with_ilu(
19//!     &matrix,
20//!     &rhs,
21//!     IluMethod::Tbem,
22//!     IluScanningDegree::Fine,
23//!     &cgs_config,
24//! );
25//! ```
26
27use ndarray::{Array1, Array2};
28use num_complex::Complex64;
29
30use crate::core::assembly::mlfmm::MlfmmSystem;
31use crate::core::assembly::slfmm::SlfmmSystem;
32use crate::core::assembly::sparse::CsrMatrix;
33
34pub use super::ilu_preconditioner::{IluMethod, IluPreconditioner, IluScanningDegree, IluSetup};
35pub use super::preconditioner::Preconditioner;
36
37/// Trait for linear operators that can perform matrix-vector products
38///
39/// This abstraction allows iterative solvers to work with:
40/// - Dense matrices (TBEM)
41/// - FMM systems (SLFMM, MLFMM)
42/// - Sparse matrices (CSR)
43/// - Preconditioners
44pub trait LinearOperator: Send + Sync {
45    /// Number of rows
46    fn num_rows(&self) -> usize;
47
48    /// Number of columns
49    fn num_cols(&self) -> usize;
50
51    /// Apply the operator: y = A * x
52    fn apply(&self, x: &Array1<Complex64>) -> Array1<Complex64>;
53
54    /// Apply the transpose operator: y = A^T * x
55    fn apply_transpose(&self, x: &Array1<Complex64>) -> Array1<Complex64>;
56
57    /// Check if the operator is square
58    fn is_square(&self) -> bool {
59        self.num_rows() == self.num_cols()
60    }
61}
62
63/// Dense matrix linear operator
64pub struct DenseOperator {
65    matrix: ndarray::Array2<Complex64>,
66}
67
68impl DenseOperator {
69    /// Create a new dense operator from a matrix
70    pub fn new(matrix: ndarray::Array2<Complex64>) -> Self {
71        Self { matrix }
72    }
73}
74
75impl LinearOperator for DenseOperator {
76    fn num_rows(&self) -> usize {
77        self.matrix.nrows()
78    }
79
80    fn num_cols(&self) -> usize {
81        self.matrix.ncols()
82    }
83
84    fn apply(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
85        self.matrix.dot(x)
86    }
87
88    fn apply_transpose(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
89        self.matrix.t().dot(x)
90    }
91}
92
93/// SLFMM linear operator
94pub struct SlfmmOperator {
95    system: SlfmmSystem,
96}
97
98impl SlfmmOperator {
99    /// Create a new SLFMM operator
100    pub fn new(system: SlfmmSystem) -> Self {
101        Self { system }
102    }
103
104    /// Get a reference to the underlying system
105    pub fn system(&self) -> &SlfmmSystem {
106        &self.system
107    }
108
109    /// Get the RHS vector
110    pub fn rhs(&self) -> &Array1<Complex64> {
111        &self.system.rhs
112    }
113}
114
115impl LinearOperator for SlfmmOperator {
116    fn num_rows(&self) -> usize {
117        self.system.num_dofs
118    }
119
120    fn num_cols(&self) -> usize {
121        self.system.num_dofs
122    }
123
124    fn apply(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
125        self.system.matvec(x)
126    }
127
128    fn apply_transpose(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
129        self.system.matvec_transpose(x)
130    }
131}
132
133/// MLFMM linear operator
134pub struct MlfmmOperator {
135    system: MlfmmSystem,
136}
137
138impl MlfmmOperator {
139    /// Create a new MLFMM operator
140    pub fn new(system: MlfmmSystem) -> Self {
141        Self { system }
142    }
143
144    /// Get a reference to the underlying system
145    pub fn system(&self) -> &MlfmmSystem {
146        &self.system
147    }
148
149    /// Get the RHS vector
150    pub fn rhs(&self) -> &Array1<Complex64> {
151        &self.system.rhs
152    }
153}
154
155impl LinearOperator for MlfmmOperator {
156    fn num_rows(&self) -> usize {
157        self.system.num_dofs
158    }
159
160    fn num_cols(&self) -> usize {
161        self.system.num_dofs
162    }
163
164    fn apply(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
165        self.system.matvec(x)
166    }
167
168    fn apply_transpose(&self, _x: &Array1<Complex64>) -> Array1<Complex64> {
169        // MLFMM transpose not yet implemented
170        unimplemented!("MLFMM transpose not yet implemented")
171    }
172}
173
174/// CSR sparse matrix linear operator
175pub struct CsrOperator {
176    matrix: CsrMatrix,
177}
178
179impl CsrOperator {
180    /// Create a new CSR operator
181    pub fn new(matrix: CsrMatrix) -> Self {
182        Self { matrix }
183    }
184
185    /// Get a reference to the underlying matrix
186    pub fn matrix(&self) -> &CsrMatrix {
187        &self.matrix
188    }
189}
190
191impl LinearOperator for CsrOperator {
192    fn num_rows(&self) -> usize {
193        self.matrix.num_rows
194    }
195
196    fn num_cols(&self) -> usize {
197        self.matrix.num_cols
198    }
199
200    fn apply(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
201        self.matrix.matvec(x)
202    }
203
204    fn apply_transpose(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
205        self.matrix.matvec_transpose(x)
206    }
207}
208
209/// Solve a linear system using CGS with the given linear operator
210///
211/// # Arguments
212/// * `operator` - Linear operator implementing the matrix-vector product
213/// * `b` - Right-hand side vector
214/// * `config` - CGS solver configuration
215///
216/// # Returns
217/// Solution vector and convergence information
218pub fn solve_cgs<O: LinearOperator>(
219    operator: &O,
220    b: &Array1<Complex64>,
221    config: &super::cgs::CgsConfig,
222) -> super::cgs::CgsSolution {
223    let matvec = |x: &Array1<Complex64>| operator.apply(x);
224    super::cgs::cgs_solve(matvec, b, None, config)
225}
226
227/// Solve a linear system using BiCGSTAB with the given linear operator
228///
229/// # Arguments
230/// * `operator` - Linear operator implementing the matrix-vector product
231/// * `b` - Right-hand side vector
232/// * `config` - BiCGSTAB solver configuration
233///
234/// # Returns
235/// Solution vector and convergence information
236pub fn solve_bicgstab<O: LinearOperator>(
237    operator: &O,
238    b: &Array1<Complex64>,
239    config: &super::bicgstab::BiCgstabConfig,
240) -> super::bicgstab::BiCgstabSolution {
241    let matvec = |x: &Array1<Complex64>| operator.apply(x);
242    super::bicgstab::bicgstab_solve(matvec, b, None, config)
243}
244
245/// Solve a linear system using the appropriate method based on system size
246///
247/// For small systems (< 1000 DOFs), uses direct LU factorization.
248/// For larger systems, uses iterative CGS with the provided operator.
249///
250/// # Arguments
251/// * `operator` - Linear operator (FMM or dense)
252/// * `b` - Right-hand side vector
253/// * `use_iterative` - Force iterative solver even for small systems
254///
255/// # Returns
256/// Solution vector
257pub fn solve_adaptive<O: LinearOperator>(
258    operator: &O,
259    b: &Array1<Complex64>,
260    use_iterative: bool,
261) -> Array1<Complex64> {
262    let n = operator.num_rows();
263
264    if !use_iterative && n < 1000 {
265        // For small systems, collect to dense and use direct solver
266        // This is a fallback - ideally the caller provides a dense matrix directly
267        eprintln!(
268            "Warning: solve_adaptive called with operator on small system (n={}). \
269             Consider using direct solver.",
270            n
271        );
272    }
273
274    // Use iterative solver
275    let config = super::cgs::CgsConfig {
276        max_iterations: n.max(100),
277        tolerance: 1e-8,
278        print_interval: 0,
279    };
280
281    let solution = solve_cgs(operator, b, &config);
282
283    if !solution.converged {
284        eprintln!(
285            "Warning: CGS did not converge after {} iterations (residual = {:.2e})",
286            solution.iterations, solution.residual
287        );
288    }
289
290    solution.x
291}
292
293/// Diagonal preconditioner from a linear operator
294///
295/// Extracts the diagonal and uses its inverse as a preconditioner.
296pub struct DiagonalPreconditioner {
297    inv_diag: Array1<Complex64>,
298}
299
300impl DiagonalPreconditioner {
301    /// Create from a CSR matrix
302    pub fn from_csr(matrix: &CsrMatrix) -> Self {
303        let diag = matrix.diagonal();
304        let inv_diag = diag.mapv(|d| {
305            if d.norm() > 1e-15 {
306                Complex64::new(1.0, 0.0) / d
307            } else {
308                Complex64::new(1.0, 0.0)
309            }
310        });
311        Self { inv_diag }
312    }
313
314    /// Create from a diagonal vector
315    pub fn from_diagonal(diag: Array1<Complex64>) -> Self {
316        let inv_diag = diag.mapv(|d| {
317            if d.norm() > 1e-15 {
318                Complex64::new(1.0, 0.0) / d
319            } else {
320                Complex64::new(1.0, 0.0)
321            }
322        });
323        Self { inv_diag }
324    }
325
326    /// Apply the preconditioner: y = M^{-1} * x
327    pub fn apply(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
328        x * &self.inv_diag
329    }
330}
331
332impl LinearOperator for DiagonalPreconditioner {
333    fn num_rows(&self) -> usize {
334        self.inv_diag.len()
335    }
336
337    fn num_cols(&self) -> usize {
338        self.inv_diag.len()
339    }
340
341    fn apply(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
342        x * &self.inv_diag
343    }
344
345    fn apply_transpose(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
346        // Diagonal matrix is symmetric
347        x * &self.inv_diag
348    }
349}
350
351/// ILU preconditioner wrapper implementing LinearOperator
352///
353/// This wraps the ILU preconditioner to work with the LinearOperator interface.
354pub struct IluOperator {
355    preconditioner: IluPreconditioner,
356    n: usize,
357}
358
359impl IluOperator {
360    /// Create from an ILU preconditioner
361    pub fn new(preconditioner: IluPreconditioner, n: usize) -> Self {
362        Self { preconditioner, n }
363    }
364
365    /// Create from a dense matrix with default settings for TBEM
366    pub fn from_tbem_matrix(matrix: &Array2<Complex64>) -> Self {
367        let n = matrix.nrows();
368        let precond =
369            IluPreconditioner::from_matrix(matrix, IluMethod::Tbem, IluScanningDegree::Fine);
370        Self {
371            preconditioner: precond,
372            n,
373        }
374    }
375
376    /// Create from a dense matrix with specified method and degree
377    pub fn from_matrix(
378        matrix: &Array2<Complex64>,
379        method: IluMethod,
380        degree: IluScanningDegree,
381    ) -> Self {
382        let n = matrix.nrows();
383        let precond = IluPreconditioner::from_matrix(matrix, method, degree);
384        Self {
385            preconditioner: precond,
386            n,
387        }
388    }
389
390    /// Get fill ratio (nnz(L+U) / n^2)
391    pub fn fill_ratio(&self) -> f64 {
392        self.preconditioner.fill_ratio()
393    }
394}
395
396impl LinearOperator for IluOperator {
397    fn num_rows(&self) -> usize {
398        self.n
399    }
400
401    fn num_cols(&self) -> usize {
402        self.n
403    }
404
405    fn apply(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
406        self.preconditioner.apply(x)
407    }
408
409    fn apply_transpose(&self, x: &Array1<Complex64>) -> Array1<Complex64> {
410        // ILU transpose application is more complex - for now just use forward
411        // This is acceptable for left preconditioning
412        self.preconditioner.apply(x)
413    }
414}
415
416/// Solve a linear system using CGS with ILU preconditioning
417///
418/// This is the **recommended** method for solving BEM systems.
419/// Diagonal preconditioning is NOT sufficient for BEM.
420///
421/// # Arguments
422/// * `matrix` - The coefficient matrix
423/// * `b` - Right-hand side vector
424/// * `method` - BEM method type (TBEM, SLFMM, MLFMM)
425/// * `degree` - ILU scanning degree (affects accuracy vs speed)
426/// * `config` - CGS solver configuration
427///
428/// # Returns
429/// Solution vector and convergence information
430///
431/// # Example
432/// ```ignore
433/// let solution = solve_with_ilu(
434///     &matrix,
435///     &rhs,
436///     IluMethod::Tbem,
437///     IluScanningDegree::Fine,
438///     &CgsConfig::default(),
439/// );
440/// assert!(solution.converged);
441/// ```
442pub fn solve_with_ilu(
443    matrix: &Array2<Complex64>,
444    b: &Array1<Complex64>,
445    method: IluMethod,
446    degree: IluScanningDegree,
447    config: &super::cgs::CgsConfig,
448) -> super::cgs::CgsSolution {
449    // Set up ILU with row scaling
450    let setup = IluPreconditioner::setup_system(matrix, method, degree);
451
452    // Scale the RHS
453    let scaled_b: Array1<Complex64> = b
454        .iter()
455        .zip(setup.row_scale.iter())
456        .map(|(&bi, &si)| bi * si)
457        .collect();
458
459    // Create matvec using scaled matrix
460    let matvec = |x: &Array1<Complex64>| setup.scaled_matrix.dot(x);
461
462    // Create preconditioner application
463    let precond_solve = |r: &Array1<Complex64>| setup.preconditioner.apply(r);
464
465    // Solve with preconditioning
466    super::cgs::cgs_solve_preconditioned(matvec, precond_solve, &scaled_b, None, config)
467}
468
469/// Solve a linear system using CGS with ILU preconditioning (using LinearOperator)
470///
471/// This version accepts any LinearOperator and builds the ILU preconditioner
472/// from the near-field matrix (for FMM systems) or the full matrix (for dense).
473///
474/// # Arguments
475/// * `operator` - Linear operator for matrix-vector products
476/// * `nearfield_matrix` - Matrix to build ILU from (near-field for FMM, full for TBEM)
477/// * `b` - Right-hand side vector
478/// * `method` - BEM method type
479/// * `degree` - ILU scanning degree
480/// * `config` - CGS configuration
481pub fn solve_with_ilu_operator<O: LinearOperator>(
482    operator: &O,
483    nearfield_matrix: &Array2<Complex64>,
484    b: &Array1<Complex64>,
485    method: IluMethod,
486    degree: IluScanningDegree,
487    config: &super::cgs::CgsConfig,
488) -> super::cgs::CgsSolution {
489    // Build ILU from near-field matrix
490    let setup = IluPreconditioner::setup_system(nearfield_matrix, method, degree);
491
492    // Scale the RHS
493    let scaled_b: Array1<Complex64> = b
494        .iter()
495        .zip(setup.row_scale.iter())
496        .map(|(&bi, &si)| bi * si)
497        .collect();
498
499    // Create scaled matvec: uses operator but applies row scaling
500    let matvec = |x: &Array1<Complex64>| {
501        let y = operator.apply(x);
502        // Apply row scaling to output
503        y.iter()
504            .zip(setup.row_scale.iter())
505            .map(|(&yi, &si)| yi * si)
506            .collect()
507    };
508
509    // Create preconditioner application
510    let precond_solve = |r: &Array1<Complex64>| setup.preconditioner.apply(r);
511
512    super::cgs::cgs_solve_preconditioned(matvec, precond_solve, &scaled_b, None, config)
513}
514
515/// Solve a TBEM system with ILU preconditioning
516///
517/// Convenience function that uses recommended settings for TBEM.
518pub fn solve_tbem_with_ilu(
519    matrix: &Array2<Complex64>,
520    b: &Array1<Complex64>,
521    config: &super::cgs::CgsConfig,
522) -> super::cgs::CgsSolution {
523    solve_with_ilu(matrix, b, IluMethod::Tbem, IluScanningDegree::Fine, config)
524}
525
526/// Information about ILU setup for diagnostics
527#[derive(Debug, Clone)]
528pub struct IluDiagnostics {
529    /// Number of nonzeros in L factor
530    pub nnz_l: usize,
531    /// Number of nonzeros in U factor
532    pub nnz_u: usize,
533    /// Fill ratio (nnz(L+U) / n^2)
534    pub fill_ratio: f64,
535    /// Threshold used for dropping
536    pub threshold_used: f64,
537}
538
539/// Create ILU diagnostics from a preconditioner
540pub fn ilu_diagnostics(
541    matrix: &Array2<Complex64>,
542    method: IluMethod,
543    degree: IluScanningDegree,
544) -> IluDiagnostics {
545    let precond = IluPreconditioner::from_matrix(matrix, method, degree);
546
547    // Get threshold for reporting
548    let threshold = match method {
549        IluMethod::Tbem => match degree {
550            IluScanningDegree::Coarse => 1.2,
551            IluScanningDegree::Medium => 1.0,
552            IluScanningDegree::Fine => 0.8,
553            IluScanningDegree::Finest => 0.6,
554        },
555        IluMethod::Slfmm => match degree {
556            IluScanningDegree::Coarse => 0.9,
557            IluScanningDegree::Medium => 0.35,
558            IluScanningDegree::Fine => 0.07,
559            IluScanningDegree::Finest => 0.01,
560        },
561        IluMethod::Mlfmm => match degree {
562            IluScanningDegree::Coarse => 0.65,
563            IluScanningDegree::Medium => 0.15,
564            IluScanningDegree::Fine => 0.05,
565            IluScanningDegree::Finest => 0.005,
566        },
567    };
568
569    IluDiagnostics {
570        nnz_l: precond.nnz_l(),
571        nnz_u: precond.nnz_u(),
572        fill_ratio: precond.fill_ratio(),
573        threshold_used: threshold,
574    }
575}
576
577// ============================================================================
578// GMRES Integration
579// ============================================================================
580
581/// Solve a linear system using GMRES with the given linear operator
582///
583/// # Arguments
584/// * `operator` - Linear operator implementing the matrix-vector product
585/// * `b` - Right-hand side vector
586/// * `config` - GMRES solver configuration
587///
588/// # Returns
589/// Solution vector and convergence information
590pub fn solve_gmres<O: LinearOperator>(
591    operator: &O,
592    b: &Array1<Complex64>,
593    config: &super::gmres::GmresConfig,
594) -> super::gmres::GmresSolution {
595    let matvec = |x: &Array1<Complex64>| operator.apply(x);
596    super::gmres::gmres_solve(matvec, b, None, config)
597}
598
599/// Solve a linear system using GMRES with ILU preconditioning
600///
601/// This is the **recommended** method for large BEM systems.
602/// GMRES provides monotonic convergence and handles non-symmetric matrices well.
603///
604/// # Arguments
605/// * `matrix` - The coefficient matrix
606/// * `b` - Right-hand side vector
607/// * `method` - BEM method type (TBEM, SLFMM, MLFMM)
608/// * `degree` - ILU scanning degree (affects accuracy vs speed)
609/// * `config` - GMRES solver configuration
610///
611/// # Returns
612/// Solution vector and convergence information
613///
614/// # Example
615/// ```ignore
616/// let config = GmresConfig::for_large_bem();
617/// let solution = gmres_solve_with_ilu(
618///     &matrix,
619///     &rhs,
620///     IluMethod::Tbem,
621///     IluScanningDegree::Fine,
622///     &config,
623/// );
624/// assert!(solution.converged);
625/// ```
626pub fn gmres_solve_with_ilu(
627    matrix: &Array2<Complex64>,
628    b: &Array1<Complex64>,
629    method: IluMethod,
630    degree: IluScanningDegree,
631    config: &super::gmres::GmresConfig,
632) -> super::gmres::GmresSolution {
633    // Set up ILU with row scaling
634    let setup = IluPreconditioner::setup_system(matrix, method, degree);
635
636    // Scale the RHS
637    let scaled_b: Array1<Complex64> = b
638        .iter()
639        .zip(setup.row_scale.iter())
640        .map(|(&bi, &si)| bi * si)
641        .collect();
642
643    // Create matvec using scaled matrix
644    let matvec = |x: &Array1<Complex64>| setup.scaled_matrix.dot(x);
645
646    // Create preconditioner application
647    let precond_solve = |r: &Array1<Complex64>| setup.preconditioner.apply(r);
648
649    // Solve with GMRES + preconditioning
650    super::gmres::gmres_solve_preconditioned(matvec, precond_solve, &scaled_b, None, config)
651}
652
653/// Solve a TBEM system with GMRES + ILU preconditioning
654///
655/// Convenience function that uses recommended settings for large TBEM problems.
656pub fn gmres_solve_tbem_with_ilu(
657    matrix: &Array2<Complex64>,
658    b: &Array1<Complex64>,
659    config: &super::gmres::GmresConfig,
660) -> super::gmres::GmresSolution {
661    gmres_solve_with_ilu(matrix, b, IluMethod::Tbem, IluScanningDegree::Fine, config)
662}
663
664/// Solve using GMRES with ILU preconditioning (using LinearOperator)
665///
666/// This version accepts any LinearOperator and builds the ILU preconditioner
667/// from the near-field matrix (for FMM systems) or the full matrix (for dense).
668///
669/// # Arguments
670/// * `operator` - Linear operator for matrix-vector products
671/// * `nearfield_matrix` - Matrix to build ILU from (near-field for FMM, full for TBEM)
672/// * `b` - Right-hand side vector
673/// * `method` - BEM method type
674/// * `degree` - ILU scanning degree
675/// * `config` - GMRES configuration
676pub fn gmres_solve_with_ilu_operator<O: LinearOperator>(
677    operator: &O,
678    nearfield_matrix: &Array2<Complex64>,
679    b: &Array1<Complex64>,
680    method: IluMethod,
681    degree: IluScanningDegree,
682    config: &super::gmres::GmresConfig,
683) -> super::gmres::GmresSolution {
684    // Build ILU from near-field matrix
685    let setup = IluPreconditioner::setup_system(nearfield_matrix, method, degree);
686
687    // Scale the RHS
688    let scaled_b: Array1<Complex64> = b
689        .iter()
690        .zip(setup.row_scale.iter())
691        .map(|(&bi, &si)| bi * si)
692        .collect();
693
694    // Create scaled matvec: uses operator but applies row scaling
695    let matvec = |x: &Array1<Complex64>| {
696        let y = operator.apply(x);
697        // Apply row scaling to output
698        y.iter()
699            .zip(setup.row_scale.iter())
700            .map(|(&yi, &si)| yi * si)
701            .collect()
702    };
703
704    // Create preconditioner application
705    let precond_solve = |r: &Array1<Complex64>| setup.preconditioner.apply(r);
706
707    super::gmres::gmres_solve_preconditioned(matvec, precond_solve, &scaled_b, None, config)
708}
709
710// ============================================================================
711// Hierarchical FMM Preconditioner Integration
712// ============================================================================
713
714/// Solve using GMRES with hierarchical FMM preconditioner
715///
716/// This method avoids the O(N²) dense matrix assembly for preconditioning.
717/// Instead, it uses block-wise LU factorization of the FMM near-field blocks.
718///
719/// # Advantages over ILU
720/// - O(N) setup cost (only diagonal blocks)
721/// - Parallel LU factorization of each block
722/// - No dense matrix extraction needed
723///
724/// # Arguments
725/// * `fmm_system` - The SLFMM system with near-field blocks
726/// * `b` - Right-hand side vector
727/// * `config` - GMRES configuration
728pub fn gmres_solve_with_hierarchical_precond(
729    fmm_system: &crate::core::assembly::slfmm::SlfmmSystem,
730    b: &Array1<Complex64>,
731    config: &super::gmres::GmresConfig,
732) -> super::gmres::GmresSolution {
733    use super::preconditioner::HierarchicalFmmPreconditioner;
734
735    // Build hierarchical preconditioner from near-field blocks
736    let precond = HierarchicalFmmPreconditioner::from_slfmm(fmm_system);
737
738    // Create matvec closure
739    let matvec = |x: &Array1<Complex64>| fmm_system.matvec(x);
740
741    // Create preconditioner application closure
742    let precond_solve = |r: &Array1<Complex64>| precond.apply(r);
743
744    // Solve with preconditioned GMRES
745    super::gmres::gmres_solve_preconditioned(matvec, precond_solve, b, None, config)
746}
747
748/// Solve using GMRES with hierarchical preconditioner (operator interface)
749///
750/// This version takes ownership of the FMM system and provides a cleaner interface.
751pub fn gmres_solve_fmm_hierarchical(
752    fmm_operator: &SlfmmOperator,
753    config: &super::gmres::GmresConfig,
754) -> super::gmres::GmresSolution {
755    use super::preconditioner::HierarchicalFmmPreconditioner;
756
757    // Build preconditioner
758    let precond = HierarchicalFmmPreconditioner::from_slfmm(fmm_operator.system());
759
760    // Get RHS
761    let b = fmm_operator.rhs();
762
763    // Create closures
764    let matvec = |x: &Array1<Complex64>| fmm_operator.apply(x);
765    let precond_solve = |r: &Array1<Complex64>| precond.apply(r);
766
767    super::gmres::gmres_solve_preconditioned(matvec, precond_solve, b, None, config)
768}
769
770// ============================================================================
771// Batched BLAS GMRES Solver
772// ============================================================================
773
774/// Solve using GMRES with batched BLAS operations
775///
776/// This version uses pre-allocated workspace and batched matrix operations
777/// for improved performance on large FMM systems.
778///
779/// # Advantages over standard matvec
780/// - Pre-allocated workspace avoids allocations in hot path
781/// - Batched operations for better cache locality
782/// - Reuses workspace across GMRES iterations
783///
784/// # Arguments
785/// * `fmm_system` - The SLFMM system
786/// * `b` - Right-hand side vector
787/// * `config` - GMRES configuration
788pub fn gmres_solve_fmm_batched(
789    fmm_system: &crate::core::assembly::slfmm::SlfmmSystem,
790    b: &Array1<Complex64>,
791    config: &super::gmres::GmresConfig,
792) -> super::gmres::GmresSolution {
793    use super::batched_blas::SlfmmMatvecWorkspace;
794    use std::cell::RefCell;
795
796    // Pre-allocate workspace in RefCell for interior mutability
797    let workspace = RefCell::new(SlfmmMatvecWorkspace::new(
798        fmm_system.num_clusters,
799        fmm_system.num_sphere_points,
800        fmm_system.num_dofs,
801    ));
802
803    // Create matvec closure using batched operations
804    // The RefCell allows mutation from within an Fn closure
805    let matvec = |x: &Array1<Complex64>| {
806        super::batched_blas::slfmm_matvec_batched(fmm_system, x, &mut workspace.borrow_mut())
807    };
808
809    super::gmres::gmres_solve(matvec, b, None, config)
810}
811
812/// Solve using GMRES with batched BLAS and ILU preconditioning
813///
814/// Combines batched matvec with ILU preconditioning for optimal performance.
815pub fn gmres_solve_fmm_batched_with_ilu(
816    fmm_system: &crate::core::assembly::slfmm::SlfmmSystem,
817    b: &Array1<Complex64>,
818    method: IluMethod,
819    degree: IluScanningDegree,
820    config: &super::gmres::GmresConfig,
821) -> super::gmres::GmresSolution {
822    use super::batched_blas::SlfmmMatvecWorkspace;
823    use std::cell::RefCell;
824
825    // Extract near-field matrix for ILU
826    let nearfield_matrix = fmm_system.extract_near_field_matrix();
827
828    // Build ILU from near-field matrix
829    let setup = IluPreconditioner::setup_system(&nearfield_matrix, method, degree);
830
831    // Scale the RHS
832    let scaled_b: Array1<Complex64> = b
833        .iter()
834        .zip(setup.row_scale.iter())
835        .map(|(&bi, &si)| bi * si)
836        .collect();
837
838    // Pre-allocate workspace for batched matvec with RefCell for interior mutability
839    let workspace = RefCell::new(SlfmmMatvecWorkspace::new(
840        fmm_system.num_clusters,
841        fmm_system.num_sphere_points,
842        fmm_system.num_dofs,
843    ));
844
845    // Create scaled matvec closure using batched operations
846    let matvec = |x: &Array1<Complex64>| {
847        let y = super::batched_blas::slfmm_matvec_batched(fmm_system, x, &mut workspace.borrow_mut());
848        // Apply row scaling to output
849        y.iter()
850            .zip(setup.row_scale.iter())
851            .map(|(&yi, &si)| yi * si)
852            .collect()
853    };
854
855    // Create preconditioner application closure
856    let precond_solve = |r: &Array1<Complex64>| setup.preconditioner.apply(r);
857
858    super::gmres::gmres_solve_preconditioned(matvec, precond_solve, &scaled_b, None, config)
859}
860
861// ============================================================================
862// Frequency-Adaptive Mesh Utilities
863// ============================================================================
864
865/// Calculate recommended mesh resolution for a given frequency
866///
867/// Based on the Nyquist criterion for BEM: mesh element size should be
868/// at most λ/6 to λ/10 for accurate results.
869///
870/// # Arguments
871/// * `frequency` - Simulation frequency in Hz
872/// * `speed_of_sound` - Speed of sound in m/s (default 343)
873/// * `elements_per_wavelength` - Number of elements per wavelength (default 6)
874///
875/// # Returns
876/// Recommended elements per meter
877pub fn recommended_mesh_resolution(
878    frequency: f64,
879    speed_of_sound: f64,
880    elements_per_wavelength: usize,
881) -> f64 {
882    let wavelength = speed_of_sound / frequency;
883    elements_per_wavelength as f64 / wavelength
884}
885
886/// Calculate mesh resolution for a frequency range
887///
888/// Uses the highest frequency to determine minimum element size.
889pub fn mesh_resolution_for_frequency_range(
890    min_freq: f64,
891    max_freq: f64,
892    speed_of_sound: f64,
893    elements_per_wavelength: usize,
894) -> f64 {
895    // Use max frequency (shortest wavelength) to determine resolution
896    recommended_mesh_resolution(max_freq, speed_of_sound, elements_per_wavelength)
897}
898
899/// Estimate element count for a room at given mesh resolution
900///
901/// Rough estimate based on room surface area.
902pub fn estimate_element_count(
903    room_dimensions: (f64, f64, f64), // (width, depth, height)
904    mesh_resolution: f64,
905) -> usize {
906    let (w, d, h) = room_dimensions;
907
908    // Total surface area of rectangular room
909    let surface_area = 2.0 * (w * d + w * h + d * h);
910
911    // Element size based on resolution (element per meter)
912    let element_size = 1.0 / mesh_resolution;
913    let element_area = element_size * element_size;
914
915    (surface_area / element_area).ceil() as usize
916}
917
918/// Adaptive mesh configuration based on frequency and room size
919pub struct AdaptiveMeshConfig {
920    /// Base mesh resolution (elements per meter)
921    pub base_resolution: f64,
922    /// Refinement factor near sources (1.0 = no refinement)
923    pub source_refinement: f64,
924    /// Refinement radius around sources (meters)
925    pub source_refinement_radius: f64,
926}
927
928impl AdaptiveMeshConfig {
929    /// Create configuration for a frequency range
930    pub fn for_frequency_range(min_freq: f64, max_freq: f64) -> Self {
931        let speed_of_sound = 343.0;
932        let base = mesh_resolution_for_frequency_range(min_freq, max_freq, speed_of_sound, 6);
933
934        Self {
935            base_resolution: base,
936            source_refinement: 1.5,
937            source_refinement_radius: 0.5,
938        }
939    }
940
941    /// Create from explicit resolution
942    pub fn from_resolution(resolution: f64) -> Self {
943        Self {
944            base_resolution: resolution,
945            source_refinement: 1.0,
946            source_refinement_radius: 0.0,
947        }
948    }
949}
950
951#[cfg(test)]
952mod tests {
953    use super::*;
954    use ndarray::array;
955
956    #[test]
957    fn test_dense_operator() {
958        let matrix = array![
959            [Complex64::new(2.0, 0.0), Complex64::new(1.0, 0.0)],
960            [Complex64::new(1.0, 0.0), Complex64::new(3.0, 0.0)],
961        ];
962
963        let op = DenseOperator::new(matrix);
964        let x = array![Complex64::new(1.0, 0.0), Complex64::new(2.0, 0.0)];
965
966        let y = op.apply(&x);
967
968        // [2 1] * [1] = [4]
969        // [1 3]   [2]   [7]
970        assert!((y[0] - Complex64::new(4.0, 0.0)).norm() < 1e-10);
971        assert!((y[1] - Complex64::new(7.0, 0.0)).norm() < 1e-10);
972    }
973
974    #[test]
975    fn test_csr_operator() {
976        let matrix = CsrMatrix::from_triplets(
977            2,
978            2,
979            vec![
980                (0, 0, Complex64::new(2.0, 0.0)),
981                (0, 1, Complex64::new(1.0, 0.0)),
982                (1, 0, Complex64::new(1.0, 0.0)),
983                (1, 1, Complex64::new(3.0, 0.0)),
984            ],
985        );
986
987        let op = CsrOperator::new(matrix);
988        let x = array![Complex64::new(1.0, 0.0), Complex64::new(2.0, 0.0)];
989
990        let y = op.apply(&x);
991
992        assert!((y[0] - Complex64::new(4.0, 0.0)).norm() < 1e-10);
993        assert!((y[1] - Complex64::new(7.0, 0.0)).norm() < 1e-10);
994    }
995
996    #[test]
997    fn test_solve_cgs_with_operator() {
998        let matrix = array![
999            [Complex64::new(4.0, 0.0), Complex64::new(1.0, 0.0)],
1000            [Complex64::new(1.0, 0.0), Complex64::new(3.0, 0.0)],
1001        ];
1002
1003        let op = DenseOperator::new(matrix.clone());
1004        let b = array![Complex64::new(1.0, 0.0), Complex64::new(2.0, 0.0)];
1005
1006        let config = super::super::cgs::CgsConfig {
1007            max_iterations: 100,
1008            tolerance: 1e-10,
1009            print_interval: 0,
1010        };
1011
1012        let solution = solve_cgs(&op, &b, &config);
1013
1014        assert!(solution.converged);
1015
1016        // Verify Ax = b
1017        let ax = matrix.dot(&solution.x);
1018        let error: f64 = (&ax - &b).iter().map(|e| e.norm_sqr()).sum::<f64>().sqrt();
1019        assert!(error < 1e-8);
1020    }
1021
1022    #[test]
1023    fn test_diagonal_preconditioner() {
1024        let diag = array![
1025            Complex64::new(2.0, 0.0),
1026            Complex64::new(4.0, 0.0),
1027        ];
1028
1029        let precond = DiagonalPreconditioner::from_diagonal(diag);
1030        let x = array![Complex64::new(1.0, 0.0), Complex64::new(2.0, 0.0)];
1031
1032        let y = precond.apply(&x);
1033
1034        // y = D^{-1} * x = [1/2, 2/4] = [0.5, 0.5]
1035        assert!((y[0] - Complex64::new(0.5, 0.0)).norm() < 1e-10);
1036        assert!((y[1] - Complex64::new(0.5, 0.0)).norm() < 1e-10);
1037    }
1038}
bem/core/solver/fmm_interface.rs

bem/core/solver/
fmm_interface.rs
