Simulation Setup
This tutorial showcases some of the other features included in the Decapodes.jl package. Currently, these features are the treatment of boundary conditions and the simulation debugger interface. To begin, we set up the same advection-diffusion problem presented in the Overview section. As before, we define the Diffusion, Advection, and Superposition components, and now include a Boundary Condition (BC) component. By convention, BCs are encoded in Decapodes by using a ∂
symbol. Below we show the graphical rendering of this boundary condition diagram, which we will use to impose a Dirichlet condition on the time derivative of concentration at the mesh boundary.
using Catlab
using DiagrammaticEquations
using Decapodes
Diffusion = @decapode begin
C::Form0
ϕ::Form1
# Fick's first law
ϕ == k(d₀(C))
end
Advection = @decapode begin
C::Form0
ϕ::Form1
V::Constant
ϕ == ∧₀₁(C,V)
end
Superposition = @decapode begin
(C, C_up)::Form0
(ϕ, ϕ₁, ϕ₂)::Form1
ϕ == ϕ₁ + ϕ₂
C_up == ⋆₀⁻¹(dual_d₁(⋆₁(ϕ)))
end
BoundaryConditions = @decapode begin
(C, C_up)::Form0
# Temporal boundary
∂ₜ(C) == Ċ
# Spatial boundary
Ċ == ∂C(C_up)
end
to_graphviz(BoundaryConditions)
As before, we compose these physics components over our wiring diagram.
compose_diff_adv = @relation (C, V) begin
diffusion(C, ϕ₁)
advection(C, ϕ₂, V)
bc(C, C_up)
superposition(ϕ₁, ϕ₂, ϕ, C_up, C)
end
draw_composition(compose_diff_adv)
DiffusionAdvection_cospan = oapply(compose_diff_adv,
[Open(Diffusion, [:C, :ϕ]),
Open(Advection, [:C, :ϕ, :V]),
Open(BoundaryConditions, [:C, :C_up]),
Open(Superposition, [:ϕ₁, :ϕ₂, :ϕ, :C_up, :C])])
DiffusionAdvection = apex(DiffusionAdvection_cospan)
to_graphviz(DiffusionAdvection)
When this is scheduled, Decapodes will apply any boundary conditions immediately after the impacted value is computed. This implementation choice ensures that this boundary condition holds true for any variables dependent on this variable, though also means that the boundary conditions on a variable have no immediate impact on the variables this variable is dependent on.
In the visualization below, we see that the final operation executed on the data is the boundary condition we are enforcing on the change in concentration.
to_graphviz(DiffusionAdvection)
Next we import the mesh we will use. In this case, we are wanting to impose boundary conditions and so we will use the plot_mesh
from the previous example instead of the mesh with periodic boundary conditions. Because the mesh is only a primal mesh, we also generate and subdivide the dual mesh.
using CombinatorialSpaces
using CairoMakie
plot_mesh = loadmesh(Rectangle_30x10())
# Generate the dual mesh
plot_mesh_dual = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3{Float64}}(plot_mesh)
# Calculate distances and subdivisions for the dual mesh
subdivide_duals!(plot_mesh_dual, Circumcenter())
fig = Figure()
ax = CairoMakie.Axis(fig[1,1], aspect = AxisAspect(3.0))
wireframe!(ax, plot_mesh)
fig
![Example block output](f8d72fc5.png)
Finally, we define our operators, generate the simulation function, and compute the simulation. Note that when we define the boundary condition operator, we hardcode the boundary indices and values into the operator itself. We also move the initial concentration to the left, so that we are able to see a constant concentration on the left boundary which will act as a source in the flow. You can find the file for boundary conditions here. The modified initial condition is shown below:
using LinearAlgebra
using ComponentArrays
using MLStyle
include("../boundary_helpers.jl")
function generate(sd, my_symbol; hodge=GeometricHodge())
op = @match my_symbol begin
:k => x -> 0.05*x
:∂C => x -> begin
boundary = boundary_inds(Val{0}, sd)
x[boundary] .= 0
x
end
x => error("Unmatched operator $my_symbol")
end
return op
end
using Distributions
c_dist = MvNormal([1, 5], [1.5, 1.5])
c = [pdf(c_dist, [p[1], p[2]]) for p in plot_mesh_dual[:point]]
fig = Figure()
ax = CairoMakie.Axis(fig[1,1], aspect = AxisAspect(3.0))
mesh!(ax, plot_mesh; color=c[1:nv(plot_mesh)])
fig
![Example block output](16edebda.png)
And the simulation result is then computed and visualized below.
using OrdinaryDiffEq
sim = eval(gensim(DiffusionAdvection))
fₘ = sim(plot_mesh_dual, generate)
velocity(p) = [-0.5, 0.0, 0.0]
v = ♭(plot_mesh_dual, DualVectorField(velocity.(plot_mesh_dual[triangle_center(plot_mesh_dual),:dual_point]))).data
u₀ = ComponentArray(C=c)
params = (V = v,)
prob = ODEProblem(fₘ, u₀, (0.0, 100.0), params)
sol = solve(prob, Tsit5());
# Plot the result
times = range(0.0, 100.0, length=150)
colors = [sol(t).C for t in times]
extrema
# Initial frame
fig = Figure()
ax = CairoMakie.Axis(fig[1,1], aspect = AxisAspect(3.0))
pmsh = mesh!(ax, plot_mesh; color=colors[1], colorrange = extrema(vcat(colors...)))
Colorbar(fig[1,2], pmsh)
framerate = 30
# Animation
record(fig, "diff_adv_right.gif", range(0.0, 100.0; length=150); framerate = 30) do t
pmsh.color = sol(t).C
end
"diff_adv_right.gif"
[ Info: Page built in 104 seconds.
[ Info: This page was last built at 2024-08-21T00:18:57.428.