Measure time accuracy with a nonlinear field

Flow / Let the nonlinear stages do some work

Measure time accuracy with a nonlinear field

Fix the grid. Refine the step. Keep the comparison above its floor.

Download runnable lesson (.zip)

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

The exact vortex has no advective tendency, so its exact viscous decay cannot establish nonlinear fourth-order time accuracy. This retained driver offers a two-mode field whose interactions do activate the nonlinear stages:

\[ \begin{aligned} \omega(x,y,0)&=\sin x\sin y\\ &\quad+\tfrac12\sin(2x)\sin(2y). \end{aligned} \]

It lives on \([0,2\pi)^2\). Its two Fourier shells have different Laplacian eigenvalues, so vorticity is not a single constant multiple of its streamfunction and advection need not vanish.

The separate shear option uses the unit square. The driver keeps those geometries explicit:

Held make_solver(CaseKind kind, int n, double dt, double nu, int rank, int nproc) {
  const auto domain =
      (kind == CaseKind::Shear) ? ns2d::make_unit_slab(n) : ns2d::make_twopi_slab(n);
  Held h;
  h.stack = std::make_unique<pfc::sim::stacks::SpectralCPUStack>(domain, rank, nproc,
                                                                 MPI_COMM_WORLD);
  h.solver =
      std::make_unique<ns2d::VorticityStreamCPU>(*h.stack, ns2d::Params{nu, dt});
  return h;
}

Keep the viscous part in an integrating factor

For one spectral mode, write \(\dot w=Lw+N(w)\) with \(L=-\nu|\mathbf k|^2\). The variable \(v=e^{-Lt}w\) satisfies

\[ \dot v=e^{-Lt}N(e^{Lt}v). \]

Classical RK4 applied to this transformed equation produces four nonlinear stages. In the implementation, exp_dt is \(e^{hL}\) and exp_half is \(e^{hL/2}\); the final combination is

const double dt6 = dt / 6.0;
for (std::size_t i = 0; i < u_hat.size(); ++i) {
  u_hat[i] =
      exp_dt[i] * u_hat[i] +
      dt6 * (exp_dt[i] * n1[i] + 2.0 * exp_half[i] * (n2[i] + n3[i]) + n4[i]);
}
sp.project_hat(u_hat);

The two-thirds state projection also applies before each stage’s products. When \(N=0\), only the exact viscous exponential remains. When \(L=0\), this reduces to classical RK4. With both parts active, the nonlinear stages matter; a stable run alone does not establish their order.

Compare with one finest-step trajectory

The grid, viscosity, initial field and physical sample times remain fixed. Only \(h\) changes. Let

\[ E(h)=\|\omega_h(T)-\omega_{h_{\mathrm{ref}}}(T)\|_{\mathrm{RMS}}. \]

The comparison reuses the common-grid Fourier helper from the spatial lesson, with equal grid sizes. There is no spatial rescaling between these trajectories. The reference is still a finite-step numerical solution, not an exact answer.

For adjacent steps related by a factor two, the reported estimate is

\[ p=\log_2\!\bigl(E(h)/E(h/2)\bigr). \]

[[nodiscard]] inline std::optional<double>
observed_order(double e_coarse, double e_fine, double floor) {
  if (!(e_coarse > floor && e_fine > floor)) return std::nullopt;
  return std::log2(e_coarse / e_fine);
}

Choose a factor-of-two ladder when interpreting this helper’s log2 formula.

The driver withholds a value when either error is at or below the configured floor, when a row failed, or when the finer comparison is the reference itself. Its roundoff label means that the floor prevented the estimate; it does not prove that floating point error is the only source of the small difference. Even above that floor, reference-step error can bias a ratio. Refine the reference before interpreting a stable estimate as a method’s asymptotic order.

Run a small nonlinear ladder

From the repository root, with a compatible installed spectral CPU SDK:

mkdir -p results
openpfc compile examples/flow/temporal-refinement/model.cpp --source-root examples/flow -o results/temporal-refinement
results/temporal-refinement --case two-mode --n 16 --dts 0.2,0.1,0.05,0.025 --times 0.8 --nu 0.05 --outdir results/flow/temporal-refinement

Use one MPI rank, as required by this comparison driver. Each step must be an integer multiple of the finest step. Sample times are rounded to finest-step indices, and those indices must be divisible by every coarser subdivision. Choose observation times that are exact multiples of every step and check the recorded time column.

The model uses the shared solver, spectral operations, initial fields, comparison helper and order helper. Keep their relative layout when using downloads; the declared source root captures their quoted includes. The CMake target remains ns2d_temporal_refine at its existing location. Its original short smoke case and numerical tests are retained.

The supplied small ladder gives the following representative differences against its finest step:

\(h\)	Vorticity RMS difference	Reported \(p\)
\(0.2\)	\(1.54\times10^{-7}\)	\(3.99\)
\(0.1\)	\(9.69\times10^{-9}\)	\(4.08\)
\(0.05\)	\(5.73\times10^{-10}\)	—
\(0.025\)	Reference	—

The last comparison is withheld because its finer member is the reference. A value slightly above four can reflect reference-step bias; it does not establish a higher-order method.

Read temporal_errors.csv and run.json under the result directory. Inspect status before interpreting error or order columns. Rows with nonfinite fields or excessive CFL are recorded as failures; their data must not be fitted into a convergence line. The table also records mean drift, divergence, energy differences and the outer-band diagnostic.

Add a finer reference step, keeping the physical times fixed. Then repeat at a finer grid. The first experiment tests sensitivity to reference time error; the second asks whether spatial truncation affects the time study. A two-level smoke test checks execution and output structure, while a meaningful order estimate needs multiple valid levels above the floor.