3D Diffusion and Advection

[GeorgeR227]

This PR is meant to serve as an initial lift to supporting 3D simulations in Decapodes. The initial commit includes an example simulation showcasing heat diffusion in a 3D cube, result shown below. This is based off a test created by @lukem12345 in CombinatorialSpaces.jl.

3d_heat.mp4

Changes to the Decapodes compiler and various changes/additions to DiagX will have to be made to support simulation generation. This PR can be considered complete once the above simulation can be run in a workflow that is essentially the same as the 2D heat diffusion simulation shown below.

using Decapodes
using DiagrammaticEquations
using CombinatorialSpaces
using GeometryBasics
using MLStyle
using ComponentArrays
using OrdinaryDiffEq
using Distributions
using CairoMakie

lx = ly = 10
s = triangulated_grid(lx, ly, 0.1, 0.1, Point3D)
sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s)
subdivide_duals!(sd, Circumcenter())

Heat = @decapode begin
    T::Form0
    D::Constant
    ∂ₜ(T) == D * Δ(T)
end

simulate = evalsim(Heat)
  
fₘ = simulate(sd, nothing)

T_dist = MvNormal([lx/2, ly/2], [1, 1])
T = [pdf(T_dist, [p[1], p[2]]) for p in sd[:point]]

D = 0.1

u₀ = ComponentArray(T=T)

constants_and_parameters = (D = D,)

tₑ = 100

prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
soln = solve(prob, Tsit5())

fig = Figure();
ax = CairoMakie.Axis(fig[1,1])
msh = CairoMakie.mesh!(ax, s, color=soln.u[end].T, colormap=:jet)
Colorbar(fig[1,2], msh)
save("2d_heat.png", fig)

[GeorgeR227]

Videos of 3D constant advection working along with a video of a projected down slab following the center of the disturbance. Note the relative lack of outward diffusion in the mass. There is some loss of mass of the advected ball as it moves along the domain, as per the bottom graph, which can be seen by the trailing wisps in the first video. However, this loss of mass seems not too great and could likely be improved with higher resolution.

3d_adv.mp4

3d_adv_slab.mp4

[DavidRobertsSBL]

What is the matter with the tests with this? I'm interested in expanding the physics library of standard equations, and something like Maxwell's equation(s) seems a fairly standard 'open' problem to add, in order to gain some experience with the system (perhaps by dropping myself in the deep end!). EM in 2+1-dimensions is a bit weird, the magnetic field is only a pseudoscalar, and there is no Gauss' law.

[GeorgeR227]

It's not so much the tests but the full range of features we would want for this PR have not yet been completed. There was an initial push a while ago to get 3D support but priorities changed and this got pushed onto the back-burner. However we're still interested in eventually getting 3D simulations fully supported so I'm glad we line up there!

I should be able to work on finishing support for this fairly soon, but if you'd like you can follow the example script in this PR and work on the equations purely using CombinatorialSpaces. There are a few things that need to be merged in for that as well but those should be completed for the most part.

[DavidRobertsSBL]

@GeorgeR227 ok, so just taking the following

using Decapodes
using DiagrammaticEquations
using CombinatorialSpaces
using GeometryBasics
using MLStyle
using ComponentArrays
using OrdinaryDiffEq
using GLMakie
using Distributions

lx = 100
ly = lz = 30

s = parallelepiped(lx = lx, ly = ly, lz = lz; tetcmd = "vpVq2a5")
sd = EmbeddedDeltaDualComplex3D{Bool, Float64, Point3D}(s)
subdivide_duals!(sd, Circumcenter())

# GLMakie.wireframe(s)

wdg_01 = dec_wedge_product(Tuple{0,1}, sd)

with the standard everything being imported (so no branches), the last line there fails:

MethodError: no method matching cache_wedge(::Type{Tuple{0, 1}}, ::EmbeddedDeltaDualComplex3D{Bool, Float64, Point{3, Float64}}, ::Type{Val{:CPU}}, ::typeof(identity), ::Nothing)
The function `cache_wedge` exists, but no method is defined for this combination of argument types.

And it suggests the 1-d and 2- versions from CombinatorialSpaces. However, I see DiscreteExteriorCalculus.jl does export this function! So not sure what's going on there.

Perhaps we should continue the discussion on my specific needs in a different place, or is here ok?

[GeorgeR227]

To get this to work, you'd have to switch over to use the gr/3d_dec_ops branch on CombinatorialSpaces.jl. The issue is that the current version doesn't allow 3D meshes to be used but the updates on that branch fix that. This is one of those things I'm planning to merge soon.

If you're in the Julia Zulip we can continue there or if not, you can email me at grauta@ufl.edu.

[Copilot]

Pull request overview

Initial lift toward 3D simulation support in Decapodes by extending operator expansion/codegen paths to dimension=3 and adding example scripts + tests around 3D DEC operators.

Changes:

• Enable gensim(...; dimension=3) and extend form allocation/type handling for 3D forms.
• Add 3D operator opening support (Lie derivatives, codifferential, de Rham Laplacians) and extend default DEC operator generation mappings.
• Add new 2D/3D diffusion/advection example scripts plus 3D operator-opening tests.

Reviewed changes

Copilot reviewed 9 out of 9 changed files in this pull request and generated 11 comments.

Show a summary per file

File	Description
src/simulation.jl	Allows dimension=3 in gensim, adds 3D form allocation mapping, and switches form detection to FORM_TYPES.
src/operators.jl	Extends default operator generation and operator opening to 3D (⋆₃, d₂, dual_d₂, new wedges, 3D Laplacians/Lie derivatives).
test/simulation_core.jl	Updates dimension fuzzing test to reflect that dimension=3 is now supported.
test/operators.jl	Adds a new testset asserting 3D operator-opening expansions match expected ACSet structure.
examples/diff_adv/heat.jl	Updates 2D heat example to use parameterized diffusion constant + Makie visualization.
examples/diff_adv/adv.jl	Adds a 2D advection example with visualization/recording.
examples/diff_adv/3d_heat.jl	Adds a 3D heat diffusion example on a tetrahedralized cube with GLMakie visualization/recording.
examples/diff_adv/3d_adv.jl	Adds a 3D advection example (DualForm0 transport) with visualization/recording.
examples/diff_adv/test_adv.jl	Adds an exploratory script comparing a codifferential-style construction in 2D vs 3D.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

is_form is defined using = return ..., which is invalid Julia syntax in a short-form function definition and will prevent this file from parsing. Drop the return and just return the boolean expression.

[Copilot]

The gensim docstring immediately above still states that dimension “Must be 1 or 2”, but this guard now allows 3. Please update the documentation so the public API contract matches the implementation.

[Copilot]

default_dec_matrix_generate can emit :⋆₁⁻¹ via dec_pair_inv_hodge(Val{1}, sd, hodge), but the only GeometricHodge method is specialized to AbstractDeltaDualComplex2D. With the new 3D operator expansions (which introduce ⋆₁⁻¹), running a 3D simulation with the default hodge=GeometricHodge() will throw a MethodError. Add a 3D GeometricHodge implementation (or change the default/emit a clearer error) so dimension=3 workflows don’t fail by default.

[Copilot]

dec_pair_wedge_product(Tuple{1,1}, ...) is only implemented for HasDeltaSet2D, but the new 3D operator expansions produce ∧₁₁ (e.g. from L₁/L₂). This will cause a MethodError during 3D codegen/execution when default_dec_matrix_generate tries to generate :∧₁₁ on a 3D complex. Add a HasDeltaSet3D (or more general HasDeltaSet) method for the (1,1) wedge product, assuming cache_wedge supports it.

[Copilot]

Point3D is used as the point type for triangulated_grid/EmbeddedDeltaDualComplex2D, but it isn’t defined anywhere in this file (most other examples define Point3D = Point3{Float64} or use Point3d). As written, this script will fail to run.

[Copilot]

The 3D Gaussian initial condition uses MvNormal([lx/2, ly/2, lx/2], ...); the third component should be centered using lz/2 (or another z-center), not lx/2. This currently mis-centers the distribution along z.

[Copilot]

This example uses Point3D (and later Point3d) but does not import GeometryBasics or define the Point3D alias. Most other examples define Point3D = Point3{Float64} (or consistently use Point3d). As written it will error when constructing EmbeddedDeltaDualComplex3D / evaluating the constant form.

[Copilot]

This example uses Point3D as the point type in triangulated_grid/EmbeddedDeltaDualComplex2D, but Point3D is never defined in this file. Other examples define it as Point3{Float64} (or use Point3d). As written, the example will error on load.

[Copilot]

This script uses Point3d and Point3D but only imports CombinatorialSpaces. Point3d/Point3{Float64} come from GeometryBasics, and Point3D isn’t defined here. As written, the example will error on load; import GeometryBasics and define the point type consistently.

[Copilot]

Dimension=3 is now accepted by gensim, but the test suite only exercises 3D at the operator-expansion level (and updates the “unsupported dimension” fuzzing test). Consider adding at least one end-to-end 3D gensim/evalsim smoke test (e.g., compile a simple 3D heat decapode and ensure the generated RHS runs on a tiny 3D complex) to catch missing runtime operator generators like wedges/inverse hodges.