A multi-physics compiler for generating numerical solvers from differential equations

We develop a tool that enables domain experts to quickly generate numerical solvers for emerging multi-physics phenomena starting from a high-level description based on ordinary/partial differential equations and their initial and boundary conditions over a symbolic spacetime domain. This “multi-physics” compiler aims to bridge the gap between problem formulation and computation, which historically has spanned years or even decades. Specialized numerical solvers in areas such as computational fluid dynamics (CFD) often present a barrier to novice end users not well-versed in the intricacies of their underlying schemes, and requiring surgical modifications when coupling with additional physical components initially not accounted for by the solver. Through the use of an intermediate language that is neutral between classical and exterior calculus, the compiler generates correct-by-construction numerical source code that offers guarantees of immutable physical principles like conservation laws at the discrete level. We present a proof of concept for the multi-physics compiler through some examples involving compilation to OpenFOAM [1]. A specific use case that the compiler is well-suited for involves equation modification approaches, where the aim is to use simple numerical schemes such as central differencing through the additional of artificial terms to the original governing equations of the multi-physics problem [2, 3, 4].

The paper proceeds as follows. Section 2 introduces an intermediate representation called “grid templates” and explains how differential equations get expanded into these constructs. Section 3 introduces grid stencils as grid templates with added information on the geometry of the computational grid. Section 4 shows an example of grid templates in action when compiling to numerical source code readable by OpenFOAM. Finally, the appendices shows how a simple system of differential equations being specified in Python, and the C files that the compiler produces for those equations.

The grid template notation represents geometry elements using typed variables. There are currently four types of element variables: volumes (v), volume centers (c), volume faces (f), and the lines between two volume centers (l). This set can be extended to include volume edges, volume vertices, volumes centered on vertices, and the faces of volumes centered on vertices.

Time is represented in grid notation as a sequence of discrete instants. We use $t$ as the time variable, such that $t+1$ represents the next time instant, and $t-1$ represents the previous time instant.

The grid template notation supports the following geometric functions, inspired by the concept of dual cell complexes from algebraic topology [5].

• •

  $\operatorname{boundaries}(l)$ - the end points of a line (printed as $\partial(l)$)

• •

  $\operatorname{boundaries}(v)$ - the faces of a volume (printed as $\partial(v)$)

• •

  $\operatorname{centroid}(f)$ - the centroid of a face

• •

  $\operatorname{dual}(c)$ - the volume containing a center

• •

  $\operatorname{dual}(f)$ - the line between two centers that passes through a face

• •

  $\operatorname{dual}(l)$ - the face that a line passes through

• •

  $\operatorname{dual}(v)$ - the center of a volume

• •

  $\operatorname{normal}(f)$ - the normal of a face

• •

  $\operatorname{normal}(l)$ - the line as a vector

• •

  $\operatorname{orientation}(c,l)$ - the orientation of an end point relative to its line ($\pm$ 1)

• •

  $\operatorname{orientation}(f,v)$ - the orientation of a face relative to its volume ($\pm$ 1)

• •

  $\operatorname{size}(f)$ - the area of a face

• •

  $\operatorname{size}(l)$ - the length of a line

• •

  $\operatorname{size}(v)$ - the volume of a volume

In our experience so far, boundaries(), centroid(), dual(), normal(), orientation(), and size() are sufficient to give us enough flexibility for defining grid templates.

Physical variables take two arguments: the geometric element that they apply to and the time instant that they apply to. For instance, $T(c,t)$ represents temperature at the given center variable and the given time variable. Likewise, $T(\operatorname{centroid}(f),t+1)$ represents temperature at the centroid of the given face variable and the next time instant after the given time variable. By convention, if the time argument is not given (e.g. $T(c)$), then it is implicitly assumed to be $t$.

Physical variables can also have properties, such as their value type (i.e., scalar, vector, or tensor), their meaning (e.g., “temperature”), their range (e.g., “positive”), and how they are initialized. Physical constants are like physical variables except for the fact that they do not take a geometric argument or a time argument.

The grid template notation is based on SymPy [6], and uses SymPy functions where possible. For instance, +, *, -, $\curl$, $\gradient$, $\divergence$, and $\cdot$ (dot product) are implemented using SymPy functions. Other functions are implemented as SymPy unknown functions.

To compile a differential equation to a computational platform, we first expand the terms of the DEs using grid templates. In order to do the expansion, we need to know what geometric element each variable applies to. Our default assumption is that variables are defined at centers, because most CFD sofware packages (including OpenFOAM) use collocated grids where variables are stored at cell centers, and our initial testbed for the compiler uses OpenFOAM as the target language. For instance, the heat equation $\divergence\gradient T=0$ is implemented as $\divergence\gradient T(c)=0$.

If the user has not provided an explicit template for a gradient like $\gradient T(c)$, then the compiler will create the default template

The first part of (2), $\operatorname{sum}(\operatorname{in}(f,\partial(\operatorname{dual}(c))),x(f,c))$, where

takes a sum over the faces of the volume that is the dual of $c$. Mathematically it is equivalent to $\sum_{f\in\partial(\operatorname{dual}(c))}{x(f,c)}$. Note that $f$ is scoped by $\operatorname{sum}$, so that if another variable named $f$ is introduced by another $\operatorname{sum}$ inside of $x(f,c)$, the embedded variable takes precedence within its scope.

The second part of (2) is a product of terms that indicates what should be summed: $T(\operatorname{centroid}(f))$ is the value of $T$ at the centroid of the face $f$ introduced by the sum, $\operatorname{orientation}(f,\operatorname{dual}(c))$ is the orientation of the face $f$ relative to its volume ($\operatorname{dual}(c)$), and $\operatorname{normal}(f))$ is the normal vector of the face $f$ scaled by the area of $f$. Each face has a globally fixed orientation. Its orientation relative to one of its volumes will be positive and its orientation relative to the adjacent volume will be negative. This means that the flux of a face will be positive for one volume and negative for the adjacent volume, as desired.

The final part of (2) is the denominator, $\operatorname{size}(\operatorname{dual}(c))$. This normalizes the sum given by the first and second part by dividing it by the volume of the cell containing $c$. The normalization is to get $\gradient T$ at $c$, the center of the volume.

To summarize, this gradient template says that $\gradient T(c)$ is implemented on a grid by taking the value of $T$ at the centroid of each face, multiplying it by the normal of the face to get a vector, summing all of the vectors to get the overall gradient, and then dividing by the volume to get the gradient at the cell center. We will discuss how to compute $T(\operatorname{centroid}(f))$ in the section on interpolation below (Section 2.10).

If the user has not provided an explicit template for a divergence such as $\divergence\mathbf{U}(c)$, then the compiler will create the default template

This template sums over face terms like the gradient template, but the face terms are the dot product of the oriented normal of the face and the divergence variable at the face instead of their multiplication. Finally, the sum is divided by the size of the volume that contains $c$ in order the get the divergence at the center of the volume.

When an operator is applied to an operator (e.g. $\divergence\gradient T(c)$), then the outer operator is expanded first. This means that the Laplacian $\divergence\gradient T(c)$ is first expanded as

Next, $\gradient T(\operatorname{centroid}(f))$ is expanded using interpolation (see Section 2.10) and substituted in. This produces

In order the determine the value of a variable at the centroid of a face (e.g. $T(\operatorname{centroid(f)})$), we interpolate the values of the variable at the centers of the neighboring volumes. On a uniform grid, this becomes

The dual of $f$ is the line between the centers of the neighboring volumes of $f$, denoted by $\partial(\operatorname{dual}(f))$). In other words, this template interpolates the variable, say $T$, at the face center by taking the average of its values at the centers of the neighboring volumes.

If we take the divergence of a product (e.g., $\divergence(\mathbf{U}(c)p(c))$), then we end up with a product of two variables at the face center (e.g., $\mathbf{U}(\operatorname{centroid}(f))p(\operatorname{centroid}(f)$). We could expand these variables separately and multiply the result, but we could also expand them together as

This is one of the ambiguities that arise when trying to discretize continuous DEs. It is a nontrivial choice, since in some cases one of the expansions leads to a conservative discretization and the other one does not.

If we take the divergence of a gradient, then we end up with a gradient at the face center (e.g., $\gradient T(\operatorname{centroid}(f))$). We could interpolate the gradient term like we did in (7). But we could also calculate the gradient directly by looking at the difference between $T(c)$ in the neighboring volumes,

In this template, the summation calculates the difference between $T$ at one end point of the line and $T$ at the other end point, where the line is between the centers of the neighboring volumes. It is a difference because the orientations of the two end points have opposite signs. The orientations of the end points are determined by the orientation of the line, and the orientation of the line is the same as the orientation of the face that is its dual. So the difference is consistent with the orientation of the face. The difference is divided by $\operatorname{size}(\operatorname{dual}(f))$ in order to determine the rate of change of $T$. Finally, the rate of change is multiplied by the unit vector of the line ($\operatorname{normal}(\operatorname{dual}(f))/\operatorname{size}(\operatorname{dual}(f))$) to turn it into a gradient vector (the normal of a line is just the line itself).

Taking the gradient directly on the face is more efficient than calculating the gradient at the cell centers and interpolating the results to the faces, and often is also the better choice for generating conservative discretizations. Hence, by default the system takes the gradient directly on the face, but the user can override this default heuristic by defining a template according to (7),

The Hodge star operator is a construct from exterior calculus. It moves a variable from one geometric element to its dual by multiplying by the size of the original geometric element and dividing by the size of its dual (note that size(c) is equal to 1). For instance, $\operatorname{Hodge}(T(v))$ approximates $T(v)$ in terms of $T(c)$ using

Similarly, $\operatorname{Hodge}(T(l))$ approximates $T(l)$ in terms of $T(f)$ using

The default expansion of the exterior derivative of a variable applied to a line ($d(T(l))$) is given by the template

The exterior derivative of a variable on a line is the difference between the variable’s values at the end points. The template calculates this by summing the end point values oriented relative to the line. Hence, one of orientations will be positive and the other one will be negative. The orientations of the end points correspond to the orientation of the line, which is defined to be consistent with the orientation of the face that is the dual of the line. The exterior derivative of a variable at a cell center ($d(T(c))$) is given by the template

The compiler translates grid templates using specialized code, which involve default patterns that can be overridden by the user. For instance, (17) shows a pattern for translating a simple interpolation,

Patterns are applied to the grid templates bottom up. If the left part of a pattern matches, then it is replaced by the right part. Any symbol that begins with a question mark is a variable that can match any expression.

Template variables (e.g., $T(c)$) get translated into C variables. Where possible, the C variables are given the same name as the template variable that they come from. The compiler also introduces variables to represent internal computations when that makes things more clear. It tries to name these variables such that they can be easily understood by the user.

After the grid templates have been translated to OpenFOAM expressions using patterns or specialized code, the compiler generates C code from the OpenFOAM expressions. For instance, “interpolation(T, linear)” is generated as “fvc::interpolate(T)” with an entry in the fvSchemes case file saying that the interpolation should be linear,

interpolationSchemes
{
    interpolate(T)    linear;
}

OpenFOAM variables are initialized in a createFields.H file that is used by the solver. The user specifies that a variable should be initialized by reading in values from the OpenFOAM case by setting the variable’s “init” property to T. The compiler then generates

Info<< "Initializing temperature T\n" << endl;

volScalarField T
        (
            IOobject
            (
                "T",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh
        );

The user can also set the “init” property to an expression, such as $p/(R\rho)$. The compiler will then generate code to initialize the variable with this value,

Info<< "Initializing temperature T\n" << endl;

volScalarField T
        (
            IOobject
            (
                "T",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionedScalar("T", p.dimensions() /
              (R.dimensions() * rho.dimensions()), scalar(0))
        );

T = p / (R * rho + dimensionedScalar("eps",
                   R.dimensions() * rho.dimensions(), 1e-10));
T.correctBoundaryConditions();

If an expression includes a divisor that is known to be positive because of user declarations (such as R * rho above), the compiler will automatically add a small value “epsilon” to the divisor in order to avoid division-by-zero errors.

The compiler automatically generates code to correct boundary conditions. If a variable like $T$ is set to an expression, then the compiler will generate “T.correctBoundaryConditions()” immediately after it. If a variable is solved for, then the generator will generate boundary correction code for the “boundary_correction” property if it exists, and otherwise for the “init” property if it is an expression. For instance, if the “init” property for rhoU is $U*rho$, then the boundary correction code looks like

rhoU.boundaryFieldRef() == U.boundaryField() * rho.boundaryField();

The compiler currently supports all temporal integration schemes built into OpenFOAM (including forward/backward Euler and Crank-Nicolson) as well as the standard fourth-order Runge-Kutta method (RK4). Rather than requiring that the user change the grid templates to reflect the desired time scheme, the compiler lets the user specify the desired time scheme directly and then the compiler generates the appropriate code for that time scheme.

For instance, if the user specifies that they want to use a built-in explicit time scheme, then the compiler will always use “fvc::” calls (explicit) rather than “fvm::” calls (implicit) for spatial operators. If the user wants to use a built-in implicit time scheme, then the compiler will generate “fvm::” calls for spatial operators on variables that are being solved for and “fvc::” otherwise. If the user specifies that they want to use the non-built-in scheme RK4, then the compiler will explicitly generate code to iterate over the various stages of the Runge-Kutta algorithm.