lobatto_fft/lib.rs
1//! # lobatto-fft
2//!
3//! High-order FFT (**HoFFT**) on Gauss–Lobatto grids and corresponding high-order solver for Poisson problems.
4//!
5//! The solver solves $(\alpha I - \beta \Delta)~u = f$ for complex-valued right-hand side on a tensor-product domain
6//! $\prod_{d=0}^{N-1} [x_{l,d}, ~ x_{l,d} + L_d]$
7//! using a grid based on Gauss–Lobatto points with $n_d$ elements of polynomial
8//! order $r_d$ per direction $d$.
9//!
10//! The method is based on the **High-order FFT (HoFFT)** of Caforio & Imperiale (2019):
11//! A single FFT of length $n_d$
12//! (or $2n_d$ for non-periodic BCs) per direction decouples the Fourier modes.
13//! Each decoupled mode tuple $(k_0, \ldots, k_{N-1})$ then requires inverting an
14//! independent $r \times r$ Hermitian system whose eigendecomposition is precomputed
15//! at setup time, giving an $\mathcal{O}(N_\mathrm{dof} \log N_\mathrm{dof})$
16//! solver with small constants.
17//!
18//! ## Modules
19//!
20//! | Module | Description |
21//! |--------|-------------|
22//! | [`hofft`] | HoFFT transform engine: [`hofft::Engine`] (1D) and [`hofft::EngineND`] (ND).|
23//! | [`solver`] | Screened Poisson solvers: [`solver::Poisson`] (1D) and [`solver::PoissonND`] (ND). Supports `Periodic`, `Dirichlet`, and `Neumann` boundary conditions per direction. |
24//!
25//! ## Quick start
26//!
27//! ```rust
28//! use lobatto_fft::solver::*;
29//! use rustfft::num_complex::Complex;
30//! use std::f64::consts::PI;
31//!
32//! // Solve (α I - β Δ) u = f on [0,1]^2 with n=8 elements of order r=4 per direction
33//! // with α = 1, β = 1 and periodic boundary conditions.
34//! let solver = PoissonND::<2>::new(
35//! [1.0; 2], // domain lengths L_d
36//! [0.0; 2], // left boundaries x_{l,d}
37//! [8; 2], // elements per direction n_d
38//! [4; 2], // polynomial orders r_d
39//! 1.0, 1.0, // α, β
40//! [BoundaryCondition::Periodic; 2], // boundary conditions
41//! );
42//!
43//! let x = solver.get_x(); // Gauss-Lobatto grid, length (8·4)² = 1024
44//!
45//! let mut f: Vec<Complex<f64>> = x
46//! .iter()
47//! .map(|&[xi, yi]| {
48//! // evaluate the right-hand side f(x,y) here (= sin(pi x) sin(pi y) for the example)
49//! Complex::new((PI * xi).sin() * (PI * yi).sin(), 0.0)
50//! })
51//! .collect();
52//!
53//! solver.solve_in_place(&mut f); // f now contains the solution u at the Gauss-Lobatto nodes
54//! ```
55
56pub mod hofft;
57pub mod solver;