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