Write the expression once, inspect the generated C++

Model development / Algebra before execution

Write the expression once, inspect the generated C++

A small symbolic expression can become ordinary C++ with an ordinary numerical test.

Download runnable lesson (.zip)

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

The coupled-gradient lesson keeps two pointwise equations in a short rhs function. A symbolic description is another useful way to author that same local law, especially when terms recur. It should leave the resulting computation easy to inspect.

Keep the earlier coefficients and add an explicit source \(s(t,x,y,z)\) to both rates:

\[ \begin{aligned} f_u &= v+a u_{xy}+s,\\ f_v &= b\Delta u+c v_{xz}\\ &\quad +d u_{yz}+s,\\ s&=r\,t(x+y+z). \end{aligned} \tag{1}\]

Here \(\Delta u=u_{xx}+u_{yy}+u_{zz}\) and \(r\) is source_rate. Setting \(r=0\) recovers the previous law. A nonzero value defines a different, forced problem. In this chapter we evaluate the law on prescribed exact fields; we do not claim those fields solve the forced PDE. The linear source is sampled at the displayed grid points and is not a periodic function across the cube’s seam. It is never Fourier-differentiated here.

Make the shared work explicit

SymPy’s common-subexpression elimination identifies algebra shared by several output expressions (SymPy Development Team, n.d.). Our fixed generator asks for canonical ordering and names its temporary cse0. The committed output is small enough to read completely:

// SPDX-FileCopyrightText: 2026 VTT Technical Research Centre of Finland Ltd
// SPDX-License-Identifier: AGPL-3.0-or-later

// Generated by generate.py with SymPy 1.14.0. Do not edit.
#pragma once
#include <openpfc/kernel/data/host_device.hpp>

namespace coupled::generated {
template <class Increments, class Position, class Local, class Parameters>
OPENPFC_INLINE_HD Increments rhs(double t, const Position &position, const Local &g,
                                 const Parameters &p,
                                 double source_rate = 0.0) noexcept {
  const double a = p.a;
  const double b = p.b;
  const double c = p.c;
  const double d = p.d;
  const double u_xx = g.u.xx;
  const double u_xy = g.u.xy;
  const double u_yy = g.u.yy;
  const double u_yz = g.u.yz;
  const double u_zz = g.u.zz;
  const double v_value = g.v.value;
  const double v_xz = g.v.xz;
  const double x = position[0];
  const double y = position[1];
  const double z = position[2];
  const double cse0 = source_rate * t * (x + y + z);
  return Increments{a * u_xy + cse0 + v_value,
                    b * (u_xx + u_yy + u_zz) + c * v_xz + cse0 + d * u_yz};
}
} // namespace coupled::generated

The one temporary contains \(r\,t(x+y+z)\), reused by both outputs. Local, Parameters, Position and Increments remain model-owned types. The generated function receives derivative values; it does not compute derivatives, allocate fields, exchange halos or choose a timestep. No symbolic library is needed when compiling or running this header.

The dependency manifest records the requested field members, parameters, position coordinates, arithmetic operation counts and header SHA-256. A smaller symbolic operation count describes this expression. It does not measure runtime, register use or GPU speed.

Wrap the expression in the existing model contract

struct GeneratedModel {
  coupled::Parameters parameters;
  double source_rate;

  coupled::Increments rhs(double time, const Located &local) const noexcept {
    return coupled::generated::rhs<coupled::Increments>(
        time, local.position, local.derivatives, parameters, source_rate);
  }
};

Located combines the derivative aggregate with physical coordinates. The wrapper supplies those arguments to the generated function and returns the same host tuple increments as the previous lesson. It keeps the framework’s rhs(time, local) contract intact.

For this experiment, a small analytic evaluator supplies exact derivatives of \(u=\cos t\cos\phi\) and \(v=(-\sin t+0.6\cos t)\cos\phi\), where \(\phi=x+2y+3z\). It has no preparation work, because it reads no neighbouring samples. This deliberately isolates generated algebra from discretization error; the previous chapter supplies actual FD and spectral evaluators. The native owned-cell sweep still performs the model call and tuple scatter:

ExactEvaluator evaluator{grid, time};
std::vector<double> du(grid.size()), dv(grid.size());
pfc::sim::for_each_interior(GeneratedModel{{}, rate}, evaluator,
                            std::make_tuple(du.data(), dv.data()), time);

An independent check is immediate:

\[ \begin{aligned} f_u&=-\sin t\cos\phi+s,\\ f_v&=(-\cos t-0.6\sin t)\\ &\quad\cdot\cos\phi+s. \end{aligned} \tag{2}\]

Figure 1: Computed generated-law rates on a grid line, followed by the identical source increments in both outputs.

Compile the C++, then try the generator

Keep equations.hpp one directory above the chapter, and generated_rhs.hpp beside model.cpp. From the chapter directory run

openpfc run model.cpp input.json --source-root .. --launcher=none
openpfc run model.cpp homogeneous.json --source-root .. --launcher=none

The starting input uses \(r=0.1\); the homogeneous input uses zero. Every owned point appears in the attempt’s results/rates.csv. Subtract the homogeneous output from the forced output: both rate differences must equal \(s\) at every point. At time=0 the source vanishes even when source_rate is nonzero.

Generation is a separate, optional Python step. Download generate.py and requirements.txt to the chapter directory. From the checkout root, for example:

python3 -m venv results/symbolic-env
results/symbolic-env/bin/pip install -r examples/coupling/symbolic-rhs/requirements.txt
results/symbolic-env/bin/python examples/coupling/symbolic-rhs/generate.py --check
results/symbolic-env/bin/python examples/coupling/symbolic-rhs/generate.py --output-dir results/generated-rhs

The requirements pin SymPy 1.14.0 and mpmath 1.3.0. The generator accepts only the two bounded expressions in its own source; it is not a parser for arbitrary C++ or a general PDE compiler. Repeating generation with different Python hash seeds gives byte-identical artifacts. --check fails if the committed header or manifest differs. Edit the expression and regenerate deliberately rather than editing the generated header.

The preserved C++ regressions exercise the host tuple protocol and compare 273 arbitrary input cases against an independent 80-digit decimal oracle. Pinned-environment generation checks reconstruct the CSE expressions and perturb each derivative input to check its sensitivity. Optional HIP tests use the existing fixed-arity device increment pack; they do not call host std::tie methods inside a device kernel.

A generated Hessian term such as \(b u_{xx}\) is not a conservative discretization of \(\partial_x(b(x)u_x)\). Spatially varying fluxes need the appropriate flux and divergence operators. Symbolic simplification does not decide a discretization or establish a physical model.

References

SymPy Development Team. n.d. Term Rewriting: Common Subexpression Detection and Collection. https://docs.sympy.org/latest/modules/rewriting.html.