Obtain the canonical source

Goal: run the same numerical model that the course derived. The model and independent verifier live together in OpenPFC’s nonlinear spectral example. The website explains the method; it does not maintain a second solver copy. These links pin the exact example revision used for the recorded reference.

Create a working directory and download the three files:

mkdir nonlinear-course
cd nonlinear-course
course_source=https://raw.githubusercontent.com/ahojukka5/OpenPFC/c96f370cdfd1e2e702410b8ae79033471e8c9b8f/examples/nonlinear_spectral
for file in allen_cahn.cpp input.json verify.py; do
  curl -fSLO "$course_source/$file"
done

Install an OpenPFC spectral CPU SDK with HeFFTe/FFTW and MPI. Keep the same compiler/MPI environment loaded for SDK probing, compilation and execution. Use a source-compatible SDK; the reference’s exact revision and capabilities are in its provenance.

openpfc doctor --sdk /path/to/openpfc-sdk
openpfc run allen_cahn.cpp input.json --sdk /path/to/openpfc-sdk --launcher=none

Replace the prefix with yours. Invoke the CLI once outside an MPI launch. It compiles the source and saves an input-specific run and attempt record. --launcher=none still runs a program that initializes MPI, with one rank. On a cluster, obtain a suitable allocation first; follow the HPC guide.

Read the complete input

{
  "n": 32,
  "cutoff": 7,
  "steps": 100,
  "dt": 0.01,
  "diffusivity": 0.1,
  "mean": 0.2,
  "a1": 0.35,
  "a3": 0.15
}

n is the sampling grid along x; cutoff is the retained Fourier bandwidth K. They are different controls. The global grid is n × 4 × 4 and the physical side lengths are \(2\pi\) on every axis. The field varies only along x. Unlike the unit-spacing diffusion starter, this model changes \(\Delta x=2\pi/N\) when N changes and holds the physical domain fixed.

The initial field is mean + a1*cos(x) + a3*cos(3*x). The input validation requires even \(16\leq N\leq128\), \(3\leq K\) with \(4K<N\), positive bounded time steps and diffusivity, and a bounded initial amplitude. Those bounds define this teaching program’s scope; they are not universal conditions for all Allen–Cahn solvers.

Physics: say what equation you are solving

The C++ source separates a physics callable:

struct AllenCahn {
  double diffusivity;
  double linear(double k_squared) const { return 1 - diffusivity * k_squared; }
  double nonlinear(double value) const { return -value * value * value; }
};

The linear symbol contains both diffusion and the \(+u\) reaction. Geometry maps integer Fourier modes to physical wavenumbers. Because \(L_x=2\pi\), the x wavenumber equals its signed integer mode here; that would change for a different domain length.

Representation: make the projection explicit

SpectralCPUStack owns the domain, decomposition, FFT and an inbox-sized real field. Spectral arrays are separately sized from fft.size_outbox(). The local Fourier iterator provides both physical wavevectors and global indices. The representation keeps only x modes \(0\leq i\leq K\) with \(j=k=0\) in r2c storage. The negative x partners are implicit.

The nonlinear operation is a visible sequence:

  1. Inverse transform
  2. Evaluate −u³
  3. Forward transform
  4. Project retained band

The inverse supplies physical values with its global \(1/V\) scaling. The forward produces raw coefficients again. Each cubic stage is projected, and the completed state is projected too. The initial field already lies in the retained band. No two-thirds mask or padded transform is hidden inside this example.

Time: recognize the predictor and correction

The implementation computes the two nonlinear spectra and advances each local outbox coefficient:

predictor[i] = exponential[i] * spectrum[i] + h_phi1[i] * first[i];
// Evaluate the projected cubic at predictor, producing second.
spectrum[i] = predictor[i] + h_phi2[i] * (second[i] - first[i]);

exponential is \(e^{hL}\), h_phi1 is \(h\varphi_1(hL)\) and h_phi2 is \(h\varphi_2(hL)\). OpenPFC supplies the stable diagonal exponential/ phi1 coefficient builder. The example computes phi2 with a short series near zero and expm1 elsewhere. The method’s two stages remain easy to compare with the derived update.

Output: inspect values and ownership

The CLI prints the attempt directory. Inside its results/:

File Meaning
diagnostics.csv Global mean, full-field RMS and energy density at every step
summary.json Grid, cutoff, rank count, final time and final diagnostics
field_rank_R.csv Rank R’s owned real cells, keyed by global (i,j,k)
modes_rank_R.csv Rank R’s owned r2c coefficients, normalized by global V
layout_rank_R.json Inclusive real inbox and complex outbox corners

The energy density is \(E/(2\pi)\) for the equivalent 1D field. The extra transverse copies cancel in global averages. Normalize mode output by \(V=16N\), not by the local inbox volume. Use the verifier before relying on a visually smooth field.

Checkpoint: who chooses the numerical representation?

If you change −u³ to −u⁵, can you leave the same cutoff check?

Not if you want the same alias-free retained-product claim. Degree five requires \((5+1)K<N\). The physics callable changes the PDE, but the representation must also be chosen for that new product. Revise the validation, independent convolution reference and error studies together. An FFT backend cannot infer the intended polynomial or error budget.

Next: check the program independently Previous: ETD time stepping