Decomposition and FFT

OpenPFC spectral workflows partition a Domain across MPI ranks and transform local data through an implementation of the FFT contract. Read the spectral-stack concept before using the exact API below.

For uniform bounded axes, see boundary transforms: explicit DCT/DST grids and periodic x/y with bounded z.

Chebyshev Helmholtz solves

The rank-local functions in openpfc/kernel/fft/chebyshev.hpp solve u'' - lambda*u = f on the Lobatto grid of [-1,1]. These are not uniform-grid DCT/DST Laplacian solves. Boundary slopes are derivatives with respect to increasing z; an outward Robin slope has positive weight at +1 and negative weight at -1.

For nonzero lambda, let J integrate a Chebyshev series with integration constant zero. The solver enforces a[k] - lambda*(J*J*a)[k] = (J*J*f)[k], k=2,…,N, and two exact coefficient boundary rows. Both extra modes of the first integral are retained through the second integration before truncating the equations. This integral formulation replaces the downward recurrence and its cancellation between exponentially growing boundary-matching columns. The finite truncation is an integrated coefficient formulation, not bitwise the previous tau recurrence or pointwise collocation equations.

Row equilibration, partial pivoting and extended-precision elimination are internal; the public values remain doubles. No small forcing modes are discarded. Singular or numerically unresolved boundary systems and nonfinite data are rejected. Negative lambda remains allowed, but pure Dirichlet/Neumann resonances at lambda=-(m*pi/2)², m>=1, are rejected explicitly. General Robin poles are only guarded by the pivot test, not a condition-number estimate. Near-resonance sensitivity and insufficient spatial resolution are not cured by this formulation. Numerical differentiation of a returned high-degree field can still amplify double-precision sampling roundoff.

Zero lambda retains the existing Poisson integrations, Neumann flux compatibility check and zero-integral gauge. Fourier-Chebyshev callers reuse this solve with lambda plus the periodic squared wavenumbers. The current direct solve uses O(N²) storage and O(N³) worst-case work per bounded line; factorization/plan caching is not provided.

The Chebyshev line transform calls pfc::fft::ComplexLineFFT, an ordinary unnormalized complex DFT. For a Lobatto series of degree n the even extension has length 2n. FFTW scaling applies: the backward transform of the forward transform is length times the input. CudaComplexLineFFT is the same contract on one rank through cufftPlan1d Z2Z. IHostFFT and the cuFFT/rocFFT pencil plans remain real-to-complex volume transforms. rocFFT and a distributed or HeFFTe line DFT are not part of this step (#141).

pfc::decomposition::Decomposition

struct Decomposition

Global geometry and rank-owned boxes for a rectangular domain.

Geometry is stored by value, independently of the source Domain’s lifetime. The process grid uses x-fastest rank order. No communication is performed by this data type; MPI halo and FFT operations use its layout separately.

The grid is const, but m_domain and m_local_boxes are public mutable data. Treat a published layout as stable while fields, halo plans or FFT plans refer to it; changing these members does not rebuild their dependent state.

Public Functions

Decomposition(const Domain &domain, const Int3 &grid)

Locally partition domain with a positive, x-fastest process grid.

Throws:

  • std::invalid_argument – If the grid is invalid or would make empty boxes.

  • std::runtime_error – If the splitter violates x-fastest rank ordering. No MPI calls are made. Geometry and boxes are owned by value.

Public Members

const std::array<int, 3> m_grid

The number of parts in each dimension.

std::vector<pfc::Box3i> m_local_boxes

Local subdomain boxes.

pfc::Domain m_domain

Global Domain this partitions.

pfc::fft::IHostFFT

struct IHostFFT : public pfc::fft::IFFTQueries

Host-container FFT (std::vector / host Field). Device backends are not this.

Subclassed by pfc::fft::LocalFFT

Public Functions

virtual void forward(const RealVector &in, ComplexVector &out) = 0

Collective, unscaled real-to-complex transform on the plan communicator.

Every rank must call the same transform in the same order with valid local buffers. The output uses get_outbox_bounds(), with x-fastest ordering; its partition can differ from the input box.

Parameters:

  • in – Real buffer; must satisfy in.size() == size_inbox()

  • out – Complex buffer; must satisfy out.size() == size_outbox()

Throws:

std::invalid_argument – when either buffer size mismatches the inbox/outbox contract (implementations report expected and actual).

virtual void backward(const ComplexVector &in, RealVector &out) = 0

Collective complex-to-real transform, scaled by 1/N globally.

N is the product of the full real-grid dimensions (not the local size). This normalization makes backward(forward(x)) recover x to roundoff. All plan ranks must participate in the same transform order.

Parameters:

  • in – Complex buffer; must satisfy in.size() == size_outbox()

  • out – Real buffer; must satisfy out.size() == size_inbox()

Throws:

std::invalid_argument – when either buffer size mismatches the outbox/inbox contract (implementations report expected and actual).

pfc::fft::CPUFFT

class CPUFFT : public pfc::fft::FFT_Impl<heffte::backend::fftw>