Boundary Conditions as Collage-Decapodes

docs/make.jl

@@ -44,6 +44,7 @@ makedocs(
     "Decapodes.jl" => "index.md",
     "Overview" => "overview.md",
     "Equations" => "equations.md",
+    "BCs & ICs" => "collages.md",
     "ASCII Operators" => "ascii.md",
     "Misc Features" => "bc_debug.md",
     "Pipe Flow" => "poiseuille.md",

docs/src/collages.md

@@ -0,0 +1,29 @@
+# Collages: Formalizing Boundary and Initial Conditions
+
+Representing physical equations via formal diagrams (i.e. a Decapode) solves many challenges in automatic simulation generation. In practice, however, setting initial and boundary conditions is a major part of the workflow, and this process itself must be formalized. We do so by implementing the notion of a *collage* of Decapodes. This is directly inspired by the [Diagrammatic view of differential equations in physics](https://arxiv.org/abs/2204.01843). See §4 for relevant theoretical discussion.
+
+## Layman's Definition: Semantics
+
+We can quickly introduce a *collage* as simply taking two Decapodes, specifying a *morphism* between the two (which point from quantities in one to quantities in the other), and interpreting this entire collection of arrows as itself a single large Decapode diagram. In contrast to the composition patterns used to specify composition between different physical components, we will not merge any quantities here. Rather, we interpret these new arrows as providing information on boundary and initial conditions. This information is passed directly to our simulation generator to help handle loading initial conditions and boundary conditions.
+
+## The Workflow
+
+### Goals
+We have a few goals for this collage-based workflow:
+- Make it quick and easy to try out different initial conditions
+- Make it quick and easy to try out different boundary conditions
+- Make it quick and easy to swap out different meshes
+- Separate the specification of the "location" of boundaries from the "values" those boundaries should take
+- Separate the specification of distributions from the properties of the mesh
+- Formalize these representations for the sake of future-proofness and consistency
+
+And of course, if this workflow is not ergonomic for your particular circumstances, you can hand-craft your initial and boundary conditions, while still using the other simulation automation that the Decapodes framework offers.
+
+### Higher-Order Functions
+To enforce a consistent workflow, we employ the programming technique of using higher-order functions. In particular, we will create a collage, and call `simulation_helper` on it. The simulation helper returns 3 things: a boundary condition **loader**, an initial condition **loader**, and a simulation **generator**.
+
+The boundary condition **loader** accepts a dictionary which maps the names of physical quantities to functions which evaluate them on an arbitrary mesh. This loader then returns a **generator** which takes in the mesh that you want to run your simulation on, and any other regular constants and parameters that you want passed to your simulation. When you execute the generator on a particular mesh, you get a named tuple of parameters that is suitable to be used with the `ODEProblem` interface.
+
+Similarly, the initial condition **loader** accepts a dictionary which maps the names of physical quantities to functions which evaluate them on an arbitrary mesh. This loader then returns a **generator** which takes in the mesh that you want to run your simulation on. When you execute the generator on a particular mesh, you get a state vector that is suitable to be used with the `ODEProblem` interface.
+
+The simulation **generator** is the output of `gensim`, with operators that apply your boundary conditions included for you. Executing the output of `gensim` on a mesh returns a simulation function that is suitable to be used with the `ODEProblem` interface.

examples/chemistry/brusselator_bounded.jl

@@ -0,0 +1,163 @@
+# Demonstrate how to encode boundary conditions in Decapodes using the "collage"
+# technique.
+using Catlab
+using Catlab.Graphics
+using CombinatorialSpaces
+using CombinatorialSpaces.ExteriorCalculus
+using Decapodes
+using MultiScaleArrays
+using MLStyle
+using OrdinaryDiffEq
+using LinearAlgebra
+using GLMakie
+using Logging
+using JLD2
+using Printf
+
+using GeometryBasics: Point2, Point3
+Point2D = Point2{Float64}
+Point3D = Point3{Float64}
+
+BrusselatorDynamics = @decapode begin
+  (U, V)::Form0
+  U2V::Form0
+  α::Constant
+  F::Parameter
+  U2V == (U .* U) .* V
+  ∂ₜ(U) == 1 + U2V - (4.4 * U) + (α * Δ(U)) + F
+  ∂ₜ(V) == (3.4 * U) - U2V + (α * Δ(U))
+end
+# Visualize. You must have graphviz installed.
+to_graphviz(BrusselatorDynamics)
+
+# This is a "discrete" Decapode, with no morphisms.
+BrusselatorBoundaries = @decapode begin
+  L::Constant
+end
+
+# Specify the BC morphism between Decapodes.
+BrusselatorBCMorphism = BCMorphism(ACSetTransformation(
+  BrusselatorBoundaries, BrusselatorDynamics,
+  Var = [1]))
+
+# This is a "discrete" Decapode, with no morphisms.
+BrusselatorInitialConditions = @decapode begin
+  (U₀, V₀)::Form0
+end
+
+# Specify the IC morphism between Decapodes.
+BrusselatorICMorphism = ICMorphism(ACSetTransformation(
+  BrusselatorInitialConditions, BrusselatorDynamics,
+  Var = [1,2]))
+
+# Wrap these morphisms into a single collage.
+BrusselatorCollage = Collage(
+  BrusselatorBCMorphism,
+  BrusselatorICMorphism)
+
+# Create the BC loader, IC loader and simulation generator.
+bc_loader, ic_loader, sim =
+  simulation_helper(BrusselatorCollage)
+
+# Specify a function that finds boundaries and values to assign, generic to any
+# mesh.
+function left_wall_constant_1(sd)
+  min_x = minimum(x -> x[1], sd[:point])
+  return (
+    findall(x -> x[1] == min_x, sd[:point]),
+    ones(count(x -> x[1] == min_x, sd[:point])))
+end
+
+# Specify a function that assigns initial values, generic to any mesh.
+# This could be an explicit function, as below, or a function that
+# interpolates data from e.g. a netcdf file.
+function parabola_in_y(sd)
+  map(sd[:point]) do (_,y)
+    22 * (y *(1-y))^(3/2)
+  end
+end
+function parabola_in_x(sd)
+  map(sd[:point]) do (x,_)
+    27 * (x *(1-x))^(3/2)
+  end
+end
+
+# Specify some mesh.
+s = loadmesh(Rectangle_30x10())
+scaling_mat = Diagonal([1/maximum(x->x[1], s[:point]),
+                        1/maximum(x->x[2], s[:point]),
+                        1.0])
+s[:point] = map(x -> scaling_mat*x, s[:point])
+s[:edge_orientation] = false
+orient!(s)
+# Visualize the mesh.
+GLMakie.wireframe(s)
+sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s)
+subdivide_duals!(sd, Circumcenter())
+
+# Demonstrate the boundary condition loader.
+# This parameter data will be passed normally for demonstration.
+F₁ = map(sd[:point]) do (x,y)
+ (x-0.3)^2 + (y-0.6)^2 ≤ (0.1)^2 ? 5.0 : 0.0
+end
+F₂ = zeros(nv(sd))
+
+bc_generator = bc_loader(Dict(
+    :L => left_wall_constant_1))
+
+p = bc_generator(sd,
+  (α = 0.001,
+   F = t -> t ≥ 1.1 ? F₂ : F₁))
+
+# Demonstrate the initial condition loader.
+ic_generator = ic_loader(Dict(
+  :U₀ => parabola_in_x,
+  :V₀ => parabola_in_y))
+u₀ = ic_generator(sd)
+
+# Generate the simulation.
+fₘ = eval(sim)(sd, default_dec_generate)
+
+# Create problem and run sim for t ∈ [0,tₑ).
+tₑ = 11.5
+
+@info("Precompiling Solver")
+prob = ODEProblem(fₘ, u₀, (0, 1e-4), p)
+soln = solve(prob, Tsit5())
+soln.retcode != :Unstable || error("Solver was not stable")
+@info("Solving")
+prob = ODEProblem(fₘ, u₀, (0, tₑ), p)
+soln = solve(prob, Tsit5())
+@info("Done")
+
+@save "brusselator_bounded.jld2" soln
+
+# Visualize the final conditions.
+GLMakie.mesh(s, color=findnode(soln(tₑ), :U), colormap=:jet)
+
+begin # BEGIN Gif creation
+frames = 100
+# Initial frame
+fig = GLMakie.Figure(resolution = (1200, 800))
+p1 = GLMakie.mesh(fig[1,2], s, color=findnode(soln(0), :U), colormap=:jet, colorrange=extrema(findnode(soln(0), :U)))
+p2 = GLMakie.mesh(fig[1,4], s, color=findnode(soln(0), :V), colormap=:jet, colorrange=extrema(findnode(soln(0), :V)))
+ax1 = Axis(fig[1,2], width = 400, height = 400)
+ax2 = Axis(fig[1,4], width = 400, height = 400)
+hidedecorations!(ax1)
+hidedecorations!(ax2)
+hidespines!(ax1)
+hidespines!(ax2)
+Colorbar(fig[1,1])
+Colorbar(fig[1,5])
+Label(fig[1,2,Top()], "U")
+Label(fig[1,4,Top()], "V")
+lab1 = Label(fig[1,3], "")
+
+# Animation
+record(fig, "brusselator_bounded.gif", range(0.0, tₑ; length=frames); framerate = 15) do t
+    p1.plot.color = findnode(soln(t), :U)
+    p2.plot.color = findnode(soln(t), :V)
+    lab1.text = @sprintf("%.2f",t)
+end
+
+end # END Gif creation

src/Decapodes.jl

@@ -22,7 +22,10 @@ export normalize_unicode, DerivOp, append_dot, unicode!, vec_to_dec!,
   AbstractMeshKey, loadmesh, Icosphere, Rectangle_30x10, Torus_30x10, Point_Map,
   Open, OpenSummationDecapodeOb, OpenSummationDecapode, unique_by, unique_by!, oapply,
   CartesianPoint, SpherePoint, r, theta, phi, TangentBasis, θhat, ϕhat,
-  #average_rewrite,
+  average_rewrite, recursive_delete_parents, contract_operators,
+  default_dec_matrix_generate, default_dec_generate, 
+  AbstractDecapodeMorphism, BCMorphism, ICMorphism,
+  AbstractCollage, Collage, simulation_helper,
   recursive_delete_parents, contract_operators,
   default_dec_matrix_generate, default_dec_generate, Term,
   @decapode
@@ -35,6 +38,7 @@ append_dot(s::Symbol) = Symbol(string(s)*'\U0307')
 include("decapodeacset.jl")
 include("language.jl")
 include("composition.jl")
+include("collages.jl")
 include("coordinates.jl")
 include("visualization.jl")
 include("simulation.jl")