Fourier / A trustworthy transform

Trust an FFT by checking its coefficients

A correct total norm can conceal incorrect phases. Compare the complex answers themselves.

Includes the model, inputs, required headers and a short run guide.

An FFT is an efficient algorithm for a specified discrete Fourier transform. On a small grid we can afford the defining sum and check the algorithm against it. This refactors example 03: the independent DFT now uses the complete global signal and global mode indices, rather than treating each MPI rank’s local values as an independent periodic signal.

Choose a signal whose phases matter

Use \(N_x=8\), \(N_y=4\), \(N_z=2\) and

\[ \begin{gathered} u(x,y,z)=1+\tfrac12\cos(2x)-\tfrac14\sin x, \\ x_j=2\pi j/8. \end{gathered} \tag{1}\]

The transverse directions repeat the same signal. For the unnormalized forward convention, the only nonzero coefficients in the stored x-half spectrum are \(\widehat u_{000}=64\), \(\widehat u_{100}=8i\), and \(\widehat u_{200}=16\). The imaginary coefficient checks the sine’s sign. The Fourier viewpoint is the same one used for differentiation in Trefethen (2000), program 5.

Transform the distributed field

pfc::sim::stacks::SpectralCPUStack stack(domain, comm.rank(), comm.size());
auto &u = stack.u();
auto &fft = stack.fft();
u.apply([](double x, double, double) { return signal(x); });
const auto original = u.vec();
std::vector<std::complex<double>> spectrum(fft.size_outbox());
fft.forward(u.vec(), spectrum);

The stack chooses compatible real-space and Fourier-space ownership. spectrum is sized to this rank’s outbox. Its length need not match the length of u or the outbox length on another rank.

Now evaluate the defining global sum for every locally owned coefficient:

pfc::fft::kspace::for_each_kpoint(
    fft.get_outbox_bounds(), domain,
    [&](std::size_t index, double, double, double, int i, int j, int k) {
      std::complex<double> reference{};
      for (int z = 0; z < nz; ++z) {
        for (int y = 0; y < ny; ++y) {
          for (int x = 0; x < nx; ++x) {
            const double angle =
                -two_pi *
                (double(i * x) / nx + double(j * y) / ny + double(k * z) / nz);
            reference += signal(two_pi * x / nx) * std::polar(1.0, angle);
          }
        }
      }
      coefficient_error =
          std::max(coefficient_error, std::abs(spectrum[index] - reference));
      output << i << ',' << j << ',' << k << ',' << spectrum[index].real() << ','
             << spectrum[index].imag() << ',' << reference.real() << ','
             << reference.imag() << '\n';
    });

This reference is deliberately expensive and small. It constructs the known signal independently at every global node. An arbitrary distributed input would require collecting or otherwise communicating those global values before the same reference could be used.

Figure 1: Real and imaginary Fourier coefficients on the transverse zero-mode line. The native FFT agrees with the independent global DFT, including the positive imaginary sine coefficient.

Complete the round trip

fft.backward(spectrum, u.vec());
double round_trip_error = 0.0;
for (std::size_t i = 0; i < original.size(); ++i) {
  round_trip_error =
      std::max(round_trip_error, std::abs(u.vec()[i] - original[i]));
}
const double local = std::max(coefficient_error, round_trip_error);
double error = 0.0;
MPI_Allreduce(&local, &error, 1, MPI_DOUBLE, MPI_MAX, comm);

The inverse applies the normalization required to recover the original field. Dividing again by 64 would be an error. Conversely, a round trip alone could miss a pair of consistently wrong transforms. The coefficient comparison and the round trip answer different questions.

Try to break the check

openpfc run model.cpp input.json --ranks=2 --launcher=mpirun

The input is empty; edit the explicit signal in model.cpp. The attempt’s results/spectrum_rank_<rank>.csv records mode indices, computed complex values, and the independent reference. Repeat on one rank: the ownership changes but the global coefficients should agree within roundoff.

Change the sine sign. Predict which stored coefficient changes, then update its independent expected value. Compare complex differences, not only squared magnitudes: \(8i\) and \(-8i\) have the same magnitude.

The executable checks a \(10^{-10}\) absolute maximum discrepancy; the test additionally verifies all 40 stored mode indices and the three known coefficients. This qualifies this small signal, not every transform size or hardware backend.

Continue with choosing an FFT implementation and its layout.

The complete C++ model and starting input are available with this lesson.

References

Trefethen, Lloyd N. 2000. Spectral Methods in MATLAB. Society for Industrial; Applied Mathematics. https://doi.org/10.1137/1.9780898719598.