Decapodes.jl: DEC Multiphysics Solver eDSL

• Decapodes.jl is a Julia-native eDSL that uses discrete exterior calculus to symbolically compose and generate multiphysics solvers.
• It integrates combinatorial spaces, diagrammatic semantics, and operadic wiring diagrams to automatically derive solver backends for PDEs.
• The framework supports high-performance solvers like geometric multigrid and iterative methods, ensuring seamless interoperability with Julia’s SciML ecosystem.

Decapodes.jl is a Julia-native, embedded domain-specific language (eDSL) for describing and automatically generating discrete exterior calculus (DEC)-based multiphysics solvers. In the broader formulation introduced as Decapodes, it is also a diagrammatic tool for representing, composing, and solving partial differential equations by combining string diagrams, an operad of wiring diagrams, hypergraph-based solver derivation, categorical data migration, and DEC discretization, with generated Julia programs handed off to established ODE solvers (Morris et al., 2024, Morris et al., 17 Aug 2025).

1. Definition, design goals, and architectural scope

The stated design goals of Decapodes.jl are to encode physics models directly in the language of exterior calculus—forms, $d$, $\star$, $\wedge$, and related operators—to remain fully compositional and symbolic until code-generation time, to leverage the combinatorial and geometric data structures of CombinatorialSpaces.jl for meshes and geometric maps, to provide high-performance solvers, notably a built-in geometric multigrid, that respect the underlying DEC structure, and to interoperate seamlessly with the Julia/SciML ecosystem, including IterativeSolvers.jl, Krylov.jl, and OrdinaryDiffEq.jl (Morris et al., 17 Aug 2025).

At the highest level, the framework is organized into four components. CombinatorialSpaces.jl provides the core data structures for simplicial complexes, including EmbeddedDeltaSet2D, boundary operators, and discrete Hodge stars. Decapodes.jl supplies the eDSL layer in which a user writes “Decapodes,” namely diagrammatic equations in exterior-calculus notation. Solver backends determine how nodes such as the inverse Laplacian are interpreted, with options including direct LU, GMG, CG, and GMRES+ILU. Time-integrators are delegated to OrdinaryDiffEq.jl or similar packages for time-dependent systems (Morris et al., 17 Aug 2025).

The earlier Decapodes formulation emphasizes the same system from a categorical and compiler-oriented perspective. It presents Decapodes as a tool for representing the relationships between variables in a PDE system, for composing systems of PDEs using an operad of wiring diagrams, and for deriving solvers using hypergraphs and string diagrams. The resulting solvers produce numerical solutions consistent with state-of-the-art open source tools, as demonstrated by benchmark comparisons with SU2 (Morris et al., 2024).

2. Diagrammatic semantics and operadic composition

In the diagrammatic language, variables are drawn as labeled nodes of the form $x : X$, meaning that $x$ is an element of the space $X$. Operators such as gradient, divergence, Hodge star, material derivative, wedge products, and sums appear as labeled edges or small boxes with ordered input ports and output ports. Summation nodes, written as $\oplus$ or $\Sigma$, collect multiple flux components. Boundary conditions are inserted as extra unary operators res or as masked assignments on boundary variables (Morris et al., 2024).

This graphical syntax is made precise by viewing the diagram as a functor $D : J \to \mathrm{El}(C)$, where $J$ is a small category describing the shape of the diagram and $\star$0 is the category of generalized elements of an analytic category $\star$1, such as vector bundles or cochains. The categorical formulation is not merely expository: it fixes the semantics of variable typing, operator incidence, and admissible composition (Morris et al., 2024).

Composition is mediated by the operad $\star$2 of undirected wiring diagrams. The objects of $\star$3 are finite sets of typed ports, such as $\star$4, $\star$5, or named variables. A morphism from $\star$6 to $\star$7 is a wiring diagram whose inner boxes are labeled by ports $\star$8 and whose outer interface is $\star$9. Composition is the usual plugging of one wiring diagram into a box of another. At the level of Decapodes, a wiring pattern $\wedge$0 is applied by taking the structured cospan or colimit in $\wedge$1 (Morris et al., 2024).

A plausible implication is that Decapodes.jl’s compositionality has two layers: a formal compositional layer supplied by the operad of wiring diagrams, and an implementation layer in which symbolic operator graphs remain intact until code generation. The papers present both layers as central rather than incidental (Morris et al., 2024, Morris et al., 17 Aug 2025).

3. Meshes, forms, and DEC operators

The principal mesh abstraction highlighted in the recent implementation-oriented account is the simplicial complex type EmbeddedDeltaSet2D{T} from CombinatorialSpaces.jl. Its data include a list of vertices with their 3D coordinates, lists of edges and triangles with face maps, and from these one automatically obtains the boundary matrices $\wedge$2 and $\wedge$3 (Morris et al., 17 Aug 2025).

Discrete forms are represented in the standard DEC manner. A 0-form is a collection of real values attached to vertices, a 1-form is a collection of values attached to edges, and a 2-form is a collection of values attached to triangles or dual cells. In code, these are ordinary Vector{Float64} or SparseVector objects indexed by simplex identifiers (Morris et al., 17 Aug 2025).

Geometric maps are encoded by GMorphism{X,Y}, which stores a piecewise-affine map $\wedge$4 by its global barycentric-coordinate matrix

$\wedge$5

in compressed-sparse-column form. The construction is functorial, with

$\wedge$6

From $\wedge$7 one derives interpolation to a fine grid by right multiplication, and restriction by row-normalization followed by transpose and multiplication (Morris et al., 17 Aug 2025).

The DEC operators exposed in Decapodes.jl include $\wedge$8, $\wedge$9, diagonal Hodge stars $x : X$0 and $x : X$1, and the dual exterior derivative $x : X$2. For a 1-form $x : X$3, the codifferential is

$x : X$4

and the Hodge Laplacian on 0-forms is

$x : X$5

The wedge product is implemented for a 1-form $x : X$6 and 0-form $x : X$7 on a 2D mesh as a 1-form, with the edge value given by the average of $x : X$8 on that edge multiplied by the midpoint-average of $x : X$9. The Lie derivative for Darcy velocity $x$0 acting on temperature $x$1 is represented as

$x$2

and is implemented component-wise using existing DEC building blocks (Morris et al., 17 Aug 2025).

4. Solver derivation and compilation pipeline

Once a PDE is specified as a Decapode $x$3, the method of lines is applied so that spatial derivatives become discrete DEC operators and the remaining system is an ODE $x$4. The function $x$5 is described by a directed hypergraph in which vertices are variables, whether state or intermediate, and hyperedges are multi-input or multi-output operators. This hypergraph is stored as a relational database over the “Decapode schema,” with tables such as Var, Op₁, Op₂, and Σ, together with typing and name tables (Morris et al., 2024).

The compilation procedure then performs a topological sort of the hypergraph. The paper describes Boolean arrays visited[v], consumed₁[e], consumed₂[e], and consumedΣ[e]; it scans unary operators whose input is visited, then binary operators whose inputs are both visited, then sum nodes whose summands are all visited, emitting an ordered list of function-call expressions. This ordered list becomes Julia source for a right-hand-side function of the form rhs!(du,u,p,t) (Morris et al., 2024).

The full pipeline proceeds from a diagram in $x$6, through composition by a UWD pattern if needed, projection onto the DEC signature, instantiation of the hypergraph database, topological-sort compilation, and finally emission of a Julia function to be wrapped into an ODEProblem and solved by DifferentialEquations.jl, for example with Tsit5(). Spatial DEC operations such as $x$7, $x$8, $x$9, $X$0, $X$1, and $X$2 are implemented once for each mesh as sparse-matrix stencils or sparse data structures, so evaluation of each operator call is $X$3 on the mesh (Morris et al., 2024).

In the implementation-oriented account, this compilation model is coupled to backend selection. Decapodes.jl builds a computation graph of form spaces and DEC operators, and the user chooses how to interpret selected nodes, especially inverse elliptic operators, by binding them to direct LU, GMG, CG, GMRES+ILU, or related solver backends. For time-dependent problems, the generated right-hand side is handed off to OrdinaryDiffEq.jl or similar packages (Morris et al., 17 Aug 2025).

5. Canonical formulations: Poisson, porous convection, and geometric multigrid

A basic model problem is the Poisson equation. In the continuous setting,

$X$4

In DEC form, the 0-form Laplacian is written

$X$5

The recent account presents illustrative pseudo-code in which the DEC operators are assembled once from a mesh and the resulting matrix is solved either by GMG or by direct LU (Morris et al., 17 Aug 2025).

For porous convection, the vector-calculus form is given as

$X$6

The corresponding DEC form is

$X$7

In Decapodes.jl one defines a “Decapode” graph, declares forms and constants, writes equations in exterior-calculus style, binds the inverse Laplacian or pressure solve to a selected backend, and time-integrates. The description explicitly states that the eDSL automatically pulls out the $X$8 node and uses the chosen solver backend to form $X$9 calls when solving the incompressibility condition (Morris et al., 17 Aug 2025).

The built-in geometric multigrid solver is exposed as a type GMG with parameters nlevels, smoother, and ν1, ν2. Its internal structure tracks a mesh hierarchy

$\oplus$0

with barycentric-coordinate matrices $\oplus$1. Interpolation is multiplication by $\oplus$2; restriction is obtained by row-normalizing $\oplus$3 and transposing; the coarse solve is a direct LU or a small iterative solve; and recursive composition yields V-, W-, or F-cycles. These building blocks are described as mirroring the diagrammatic formalism and being automatically assembled once solver=GMG(...) is chosen (Morris et al., 17 Aug 2025).

6. Numerical results, ecosystem integration, and present limitations

The numerical results reported for the geometric multigrid implementation are problem-specific and comparatively favorable in elliptic and porous-media settings. For the Poisson problem, the reported behavior includes linear scaling in V-cycle count and number of vertices, residual decay exponential in cycle count, and a time per V-cycle of approximately $\oplus$4 overhead on an M2 MacBook Air. As a stand-alone solver, using 5 W-cycles with Gauss-Seidel smoothing, the total time is approximately $\oplus$5 with relative residual approximately $\oplus$6 on 7 levels. As a preconditioner for CG, the total solve time is approximately $\oplus$7 versus ILU0-CG at approximately $\oplus$8, with residual $\oplus$9 versus $\Sigma$0, described as a 78% reduction (Morris et al., 17 Aug 2025).

For the porous convection problem, the reported domain is a rectangular disk with finest mesh 16 641 vertices. The pressure solver backends are compared against direct LU as reference: direct LU has runtime $\Sigma$1, GMRES+ILU0 has runtime $\Sigma$2 with RMSE approximately $\Sigma$3 versus LU, and DEC GMG has runtime $\Sigma$4 with RMSE approximately $\Sigma$5, described as 97.8% lower error. Over time, the temperature RMSE of DEC GMG remains an order of magnitude below GMRES+ILU0 (Morris et al., 17 Aug 2025).

The 2024 paper places these implementation results in a broader benchmark context against SU2. Its runtime table reports, in seconds of compute time per second of simulated time, 50400 for conjugate heat transfer with Decapodes versus 94.3 for SU2, and 28600 for buoyancy with Decapodes versus 9.08 for SU2. The accompanying interpretation is explicit: both SU2 and Decapodes reproduce the key physics, including wakes and convection rolls, with low average error, less than 5% in most fields, but SU2’s implicit timestepping gives it much larger stable $\Sigma$6, hence SU2 is $\Sigma$7–$\Sigma$8 faster per simulated second than Decapodes, which currently only has explicit methods (Morris et al., 2024).

Integration with Julia iterative methods is presented as a deliberate architectural feature. Each linear-solve node, such as an inverse Hodge Laplacian, can be bound to any linear solver implementing the standard Julia LinearAlgebra or IterativeSolvers.jl interface. The examples include preconditioned CG and preconditioned GMRES via Krylov.jl, and when used as a preconditioner the DEC GMG object itself implements the LinearMap interface so that it can be passed directly to CG, GMRES, or related Krylov routines (Morris et al., 17 Aug 2025).

The published literature also distinguishes between the formal compilation framework and the stability of a user-facing API. One account states that the paper does not give the final public API and sketches plausible code snippets, whereas the later paper provides illustrative pseudo-code centered on concrete DEC and GMG workflows. This suggests that the most mature, fully specified contributions lie in the categorical semantics, the compiler architecture, and the DEC solver backends, while the public programming surface continues to be documented through worked examples rather than a finalized interface specification (Morris et al., 2024, Morris et al., 17 Aug 2025).