Generate the model you will study
Goal: distinguish a correct discrete implementation from an accurate approximation to the PDE. You need an installed CPU SDK and its matching compiler and MPI environment. Keep that environment loaded for both compilation and execution; an SDK built with a newer C++ runtime may not work with the system compiler.
openpfc doctor --sdk /path/to/openpfc-sdk
openpfc init diffusion-course --model=diffusion
cd diffusion-course
openpfc run model.cpp input.json --sdk /path/to/openpfc-sdk --launcher=noneReplace the SDK prefix with yours. Invoke openpfc once outside an MPI launch; the CLI handles compilation and model launches. On a cluster, do this inside an appropriate allocation using the HPC guide. The model source is the canonical template generated by the installed CLI. We link the complete source used for this reference, rather than maintain another copy of the model in the website. If your SDK has a newer template, inspect its source before applying these checks.
The default input.json is small enough to inspect completely:
{
"cells": 16,
"steps": 10,
"dt": 0.1,
"diffusivity": 1.0
}Translate the units before comparing curves
The browser used \(L=2\pi\). This starter uses a periodic cube of cells points per axis, unit spacing \(h=1\), and therefore side length \(L=\text{cells}\). Its initial condition is
\[u(x,y,z,0)=1+0.1\cos(2\pi x/L).\]
It is constant in \(y\) and \(z\), so the analytic solution still contains one decaying cosine. For the default input, \(k=2\pi/16\), \(D=1\), and \(T=10\times0.1=1\). To compare with the browser’s phase coordinate, use \(\xi=2\pi x/L\); spatial derivatives acquire factors of \(2\pi/L\). Do not insert \(k=1\) into the physical decay formula for this unit grid. Changing cells alone changes the domain length. It is not the fixed-domain refinement performed in the error chapter.
Read the implementation in layers
The physics callable contains the PDE’s right-hand side:
struct Derivatives {
double xx{}, yy{}, zz{};
};
struct Diffusion {
double diffusivity;
double rhs(double /*time*/, const Derivatives &g) const {
return diffusivity * (g.xx + g.yy + g.zz);
}
};The domain describes the global grid. Decomposition assigns owned cells to ranks. A field stores those cells and a one-cell halo. The gradient operator computes second-order derivatives from neighbouring values; the time stepper advances the owned cells using the resulting rates. The full source connects them in this order:
halo.exchange();
(void)stepper.step(step * dt, u);Exchange the current field’s neighbours before evaluating the stencil. OwnedEulerStepper updates the owned values; ghost cells are refreshed by the next exchange. The starter’s seven-point 3D Laplacian gives
\[u^{n+1}_{i,j,k}=u^n_{i,j,k}+r\bigl( u^n_{i-1,j,k}+u^n_{i+1,j,k}+u^n_{i,j-1,k}+u^n_{i,j+1,k} +u^n_{i,j,k-1}+u^n_{i,j,k+1}-6u^n_{i,j,k}\bigr).\]
Here \(r=D\Delta t/h^2=0.1\). The general 3D Euler bound is \(r\leq1/6\), which the template enforces. Although our particular initial field varies only in \(x\), the validation protects the full 3D grid of possible modes.
Derive an independent discrete answer
For this mode, the \(y\) and \(z\) differences vanish and one step multiplies the cosine amplitude by
\[g=1-4D\Delta t\sin^2(\pi/L).\]
After \(M\) steps, compare the saved field against
\[u^M_{i,j,k}=1+0.1g^M\cos(2\pi i/L).\]
This formula checks the actual FD2 + Euler update. The continuous PDE instead predicts amplitude \(0.1e^{-D(2\pi/L)^2T}\). Both are useful: one checks the implementation, the other exposes its approximation error.
Inspect a checked run
We generated the unmodified starter from a clean SDK source at 02ad8c37, compiled a CPU Release SDK with GNU 15.2.0 and Open MPI 5.0.10, and ran the default input on one and two ranks on 25 October 2026. These are recorded educational reference outputs, not timing measurements.
| Quantity at T = 1 | One rank | Two ranks |
|---|---|---|
| Unique owned cells | 4,096 | 4,096 |
| Mean, computed from final fields | 1.0 | 1.0 |
| Measured cosine amplitude | 0.08577765043852273 | 0.08577765043852273 |
| Discrete formula amplitude | 0.08577765043852266 | 0.08577765043852266 |
| PDE amplitude | 0.08570898111217012 | 0.08570898111217012 |
| Maximum field error against discrete formula | \(2.23\times10^{-16}\) or less | \(2.23\times10^{-16}\) or less |
The amplitude difference from the PDE is about \(6.87\times10^{-5}\). It is much larger than the roundoff-level difference from the discrete formula. That is expected discretization error, not an MPI defect.
Download the actual one-rank diagnostics, two-rank diagnostics, final x slice, comparison report, and source/toolchain provenance with checksums. The comparison checked the full 3D fields, not only the displayed slice. The SDK build manifest identifies the source and capabilities; this run exercised the CPU finite-difference path, not spectral or GPU execution.
Check your own output
The CLI prints an Attempt: directory. Open results/summary.json and results/diagnostics.csv there. The saved l2 is a root-mean-square norm of the full field:
\[\mathrm{l2}=\sqrt{\langle u^2\rangle}=\sqrt{1+a^2/2}.\]
It tends to 1 because the mean remains 1. The fluctuation RMS is \(|a|/\sqrt2\), which tends to zero. A nearly unchanged full-field norm can hide a substantial change in the small oscillation.
Download check-results.py outside the generated model directory, then run it with Python 3. Replace the example attempt path with the one printed by your run:
python3 ../check-results.py input.json /path/to/attempt/resultsThe standard-library checker expects the unmodified starter. It checks all diagnostics against the discrete mode, every owned coordinate/value, the global mean and amplitude, and complete nonoverlapping ownership. It fails on missing or duplicate cells, non-finite values or discrepancies above an absolute tolerance of \(2\times10^{-11}\). A changed physics callable needs a changed independent reference, not simply a relaxed tolerance.
Checkpoint: a successful run with the wrong reference
The field disagrees with the PDE amplitude by 6.87 × 10⁻⁵ but agrees with the discrete formula near roundoff. What should you change?
The observed difference is consistent with FD2 + Euler discretization. First reduce \(\Delta t\) while doubling the number of steps to hold \(T\) fixed; predict the discrete amplitude before rerunning. A residual spatial error will remain. A fixed-domain grid refinement additionally needs an explicit spacing/domain change in the model; increasing cells alone here changes \(L\). Do not alter MPI communication merely to force the coarse calculation to match the continuous solution.
To keep an attempt beyond cache cleanup, export its run record. The same guide explains openpfc bundle for a compiled CPU model.