Merge branch 'main' into llm/docs_fixes

Project.toml

@@ -7,7 +7,7 @@ authors = [
   "Luke Morris",
   "George Rauta",
 ]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 AlgebraicRewriting = "725a01d3-f174-5bbd-84e1-b9417bad95d9"

docs/make.jl

@@ -43,6 +43,7 @@ makedocs(
   pages     = Any[
     "Decapodes.jl" => "index.md",
     "Overview" => "overview.md",
+    "Equations" => "equations.md",
     "Misc Features" => "bc_debug.md",
     "Pipe Flow" => "poiseuille.md",
     "Glacial Flow" => "ice_dynamics.md",

docs/src/equations.md

@@ -0,0 +1,96 @@
+# Simple Equations
+
+This tutorial shows how to use Decapodes to represent simple equations. These aren't using any of the Discrete Exterior Calculus or CombinatorialSpaces features of Decapodes. They just are a reference for how to build equations with the `@decapodes` macro and see how they are stored as ACSets.
+
+```@example Harmonic
+using Catlab
+using Catlab.Graphics
+using CombinatorialSpaces
+using CombinatorialSpaces.ExteriorCalculus
+using Decapodes
+```
+
+The harmonic oscillator can be written in Decapodes in at least three different ways.
+
+
+```@example Harmonic
+oscillator = @decapode begin
+  X::Form0
+  V::Form0
+
+  ∂ₜ(X) == V
+  ∂ₜ(V) == -k(X)
+end
+```
+
+The default representation is a tabular output as an ACSet. The tables are `Var` for storing variables (`X`) and their types (`Form0`). `TVar` for identifying a subset of variables that are the tangent variables of the dynamics (`Ẋ`). The unary operators are stored in `Op1` and binary operators stored in `Op2`. If a table is empty, it doesn't get printed.
+
+Even though a diagrammatic equation is like a graph, there are no edge tables, because the arity (number of inputs) and coarity (number of outputs) is baked into the operator definitions.
+
+You can also see the output as a directed graph. The input arrows point to the state variables of the system and the output variables point from the tangent variables. You can see that I have done the differential degree reduction from  `x'' = -kx` by introducing a velocity term `v`. Decapodes has some support for derivatives in the visualization layer, so it knows that `dX/dt` should be called `Ẋ` and that `dẊ/dt` should be called `Ẋ̇`.
+
+```@example Harmonic
+to_graphviz(oscillator)
+```
+
+In the previous example, we viewed negation and transformation by `k` as operators. Notice that `k` appears as an edge in the graph and not as a vertex. You can also use a 2 argument function like multiplication (`*`). With a constant value for `k::Constant`. In this case you will see `k` enter the diagram as a vertex and multiplication with `*` as a binary operator.
+
+```@example Harmonic
+oscillator = @decapode begin
+  X::Form0
+  V::Form0
+
+  k::Constant
+
+  ∂ₜ(X) == V
+  ∂ₜ(V) == -k*(X)
+end
+```
+
+This gives you a different graphical representation as well. Now we have the cartesian product objects which represent a tupling of two values.
+
+```@example Harmonic
+to_graphviz(oscillator)
+```
+
+You can also represent negation as a multiplication by a literal `-1`.
+
+```@example Harmonic
+oscillator = @decapode begin
+  X::Form0
+  V::Form0
+
+  k::Constant
+
+  ∂ₜ(X) == V
+  ∂ₜ(V) == -1*k*(X)
+end
+```
+
+Notice that the type bubble for the literal one is `ΩL`. This means that it is a literal. The literal is also used as the variable name.
+
+```@example Harmonic
+infer_types!(oscillator)
+to_graphviz(oscillator)
+```
+
+We can allow the material properties to vary over time by changing `Constant` to `Parameter`. This is how we tell the simulator that it needs to call `k(t)` at each time step to get the updated value for `k` or if it can just reuse that constant `k` from the initial time step.
+
+```@example Harmonic
+oscillator = @decapode begin
+  X::Form0
+  V::Form0
+
+  k::Parameter
+
+  ∂ₜ(X) == V
+  ∂ₜ(V) == -1*k*(X)
+end
+```
+
+```@example Harmonic
+infer_types!(oscillator)
+to_graphviz(oscillator)
+```
+
+Often you will have a linear material where you are scaling by a constant, and a nonlinear version of that material where that scaling is replaced by a generic nonlinear function. This is why we allow Decapodes to represent both of these types of equations.

examples/chemistry/brusselator_teapot.jl

@@ -0,0 +1,79 @@
+# Imports
+using Catlab
+using Catlab.Graphics
+using CombinatorialSpaces
+using CombinatorialSpaces.ExteriorCalculus
+using Decapodes
+using MultiScaleArrays
+using MLStyle
+using OrdinaryDiffEq
+using LinearAlgebra
+using CairoMakie
+using Logging
+using JLD2
+using Printf
+using GeometryBasics: Point3
+Point3D = Point3{Float64}
+
+# Define Model
+Brusselator = @decapode begin
+  (U, V)::Form0
+  α::Constant
+  F::Parameter
+
+  U2V == (U*U) * V
+  ∂ₜ(U)== 1 + U2V - (4.4 * U) + (α * Δ(U)) + F
+  ∂ₜ(V) == (3.4 * U) - U2V + (α * Δ(U))
+end
+infer_types!(Brusselator)
+resolve_overloads!(Brusselator)
+
+# Define Domain
+download("https://graphics.stanford.edu/courses/cs148-10-summer/as3/code/as3/teapot.obj", "teapot.obj")
+s = EmbeddedDeltaSet2D("teapot.obj")
+sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s)
+subdivide_duals!(sd, Circumcenter())
+fig,ax,ob = wireframe(sd)
+save("teapot_subdivided.png", fig)
+
+# Create initial data.
+U = map(p -> abs(p[2]), point(s))
+V = map(p -> abs(p[1]), point(s))
+F₁ = map(sd[:point]) do (_,_,z)
+  z ≥ 0.8 ? 5.0 : 0.0
+end
+F₂ = zeros(nv(sd))
+
+constants_and_parameters = (
+  α = 0.001,
+  F = t -> t ≥ 1.1 ? F₂ : F₁)
+u₀ = construct(PhysicsState, [VectorForm(U), VectorForm(V)], Float64[], [:U, :V])
+
+# Run
+function generate(sd, my_symbol; hodge=GeometricHodge()) end
+sim = eval(gensim(Brusselator))
+fₘ = sim(sd, generate)
+tₑ = 1.15
+prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
+soln = solve(prob, Tsit5())
+
+@save "brusselator_teapot.jld2" soln
+
+# Create side-by-side GIF
+begin
+frames = 800
+# Initial frame
+fig = CairoMakie.Figure(resolution = (1200, 1200))
+p1 = CairoMakie.mesh(fig[1,1], s, color=findnode(soln(0), :U), colormap=:jet, colorrange=extrema(findnode(soln(0), :U)))
+p2 = CairoMakie.mesh(fig[2,1], s, color=findnode(soln(0), :V), colormap=:jet, colorrange=extrema(findnode(soln(0), :V)))
+Colorbar(fig[1,2])
+Colorbar(fig[2,2])
+
+# Animation
+record(fig, "brusselator_teapot.gif", range(0.0, tₑ; length=frames); framerate = 30) do t
+    p1.plot.color = findnode(soln(t), :U)
+    p2.plot.color = findnode(soln(t), :V)
+end
+
+end
+