Add docs fixes (#128)

* Add grid meshes script and latex halfar

* Update BC, Overview, and part of Poise docs

* Add BC and Overview fixes, partial fixes to Poiseuille

* Fix reference to MSArray, and grid mesh path

* Pass V as a state var

* Update Poise

---------

Co-authored-by: James <[email protected]>

docs/make.jl

@@ -42,10 +42,10 @@ makedocs(
   checkdocs = :none,
   pages     = Any[
     "Decapodes.jl" => "index.md",
-    #"Overview" => "overview.md",
-    #"Misc Features" => "bc_debug.md",
-    #"Pipe Flow" => "poiseuille.md",
+    "Overview" => "overview.md",
     "Equations" => "equations.md",
+    "Misc Features" => "bc_debug.md",
+    "Pipe Flow" => "poiseuille.md",
     "Glacial Flow" => "ice_dynamics.md",
     "Budyko-Sellers-Halfar" => "budyko_sellers_halfar.md",
 #    "Examples" => Any[

docs/src/bc_debug.md

@@ -13,43 +13,45 @@ impose a Dirichlet condition on the time derivative of concentration at the
 mesh boundary.
 
 ```@example Debug
-using Decapodes, Decapodes.Diagrams, Decapodes.OpenDiagrams
-using Catlab.Present, Catlab.Graphics, Catlab.Programs
+using Catlab
+using Catlab.Graphics
+using Decapodes
 
-@present DiffusionQuantities <: Decapodes2D begin
-  k::Hom(Form1(X), Form1(X))    # diffusivity (usually scalar multiplication)
-  ∂C::Hom(Form0(X), Form0(X))   # concentration boundary condition
-end;
-
-Diffusion = @decapode DiffusionQuantities begin
-  C::Form0{X}
-  ϕ::Form1{X}
+Diffusion = @decapode begin
+  C::Form0
+  ϕ::Form1
 
   # Fick's first law
-  ϕ ==  k(d₀{X}(C))
+  ϕ ==  k(d₀(C))
 end
 
-Advection = @decapode DiffusionQuantities begin
-  C::Form0{X}
-  (V, ϕ)::Form1{X}
-  ϕ == ∧₀₁{X}(C,V)
+Advection = @decapode begin
+  C::Form0
+  ϕ::Form1
+  V::Constant
+
+  ϕ == ∧₀₁(C,V)
 end
 
-Superposition = @decapode DiffusionQuantities begin
-  (C, Ċ)::Form0{X}
-  (ϕ, ϕ₁, ϕ₂)::Form1{X}
+Superposition = @decapode begin
+  (C, C_up)::Form0
+  (ϕ, ϕ₁, ϕ₂)::Form1
 
   ϕ == ϕ₁ + ϕ₂
-  Ċ == ⋆₀⁻¹{X}(dual_d₁{X}(⋆₁{X}(ϕ)))
-  ∂ₜ{Form0{X}}(C) == Ċ
+  C_up == ⋆₀⁻¹(dual_d₁(⋆₁(ϕ)))
 end
 
-BoundaryConditions = @decapode DiffusionQuantities begin
-  (Ċ, Ċ_bound)::Form0{X}
-  ∂C(Ċ) == Ċ_bound
+BoundaryConditions = @decapode begin
+  (C, C_up)::Form0
+
+  # Temporal boundary
+  ∂ₜ(C) == Ċ
+
+  # Spatial boundary
+  Ċ == ∂C(C_up)
 end
 
-draw_diagram(BoundaryConditions)
+to_graphviz(BoundaryConditions)
 ```
 
 As before, we compose these physics components over our wiring diagram.
@@ -59,24 +61,22 @@ As before, we compose these physics components over our wiring diagram.
 compose_diff_adv = @relation (C, V) begin
   diffusion(C, ϕ₁)
   advection(C, ϕ₂, V)
-  bc(Ċ)
-  superposition(ϕ₁, ϕ₂, ϕ, Ċ, C)
+  bc(C, C_up)
+  superposition(ϕ₁, ϕ₂, ϕ, C_up, C)
 end
 
-to_graphviz(compose_diff_adv, box_labels=:name, junction_labels=:variable,
-            graph_attrs=Dict(:start => "2"))
+to_graphviz(compose_diff_adv, box_labels=:name, junction_labels=:variable, prog="circo")
 ```
 
-
-
 ```@example Debug
-DiffusionAdvection = oapply(compose_diff_adv,
-                  [OpenDiagram(Diffusion, [:C, :ϕ]),
-                   OpenDiagram(Advection, [:C, :ϕ, :V]),
-                   OpenDiagram(BoundaryConditions, [:Ċ]),
-                   OpenDiagram(Superposition, [:ϕ₁, :ϕ₂, :ϕ, :Ċ, :C])])
-
-draw_diagram(DiffusionAdvection.functor)
+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
@@ -85,26 +85,23 @@ 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.
 
-Below, we see the generated schedule, which shows that the final operation
+In the visualization below, wee see that the final operation
 executed on the data is the boundary condition we are enforcing on the change
 in concentration.
 
 
 ```@example Debug
-using Decapodes.Schedules
-
-explicit_ts = diag2dwd(DiffusionAdvection.functor)
-to_graphviz(explicit_ts, orientation=LeftToRight)
+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.
 
+`Rectangle_30x10` is a default mesh that is downloaded via `Artifacts.jl` when a user installs Decapodes. Via CombinatorialSpaces.jl, we can instantiate any `.obj` file of triangulated faces as a simplicial set.
 
 ```@example Debug
-using Catlab.CategoricalAlgebra
 using CombinatorialSpaces, CombinatorialSpaces.DiscreteExteriorCalculus
 using CairoMakie
 
@@ -115,35 +112,40 @@ plot_mesh_dual = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3{Float64}}(plot
 # Calculate distances and subdivisions for the dual mesh
 subdivide_duals!(plot_mesh_dual, Circumcenter())
 
-
 fig, ax, ob = wireframe(plot_mesh)
 ax.aspect = AxisAspect(3.0)
 fig
 ```
 
 Finally, we define our operators, generate the simulation function, and compute
-the simulation. Note that when we define the boudary condition operator, we
+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. The
 modified initial condition is shown below:
 
-
 ```@example Debug
 using LinearAlgebra
-using Decapodes.Examples, Decapodes.Simulations
+using MultiScaleArrays
+using MLStyle
 using CombinatorialSpaces.DiscreteExteriorCalculus: ∧
-
-funcs = sym2func(plot_mesh_dual)
-
-funcs[:k] = Dict(:operator => 0.05 * I(ne(plot_mesh_dual)), :type => MatrixFunc())
-funcs[:⋆₁] = Dict(:operator => ⋆(Val{1}, plot_mesh_dual, hodge=DiagonalHodge()),
-                  :type => MatrixFunc());
-funcs[:∧₀₁] = Dict(:operator => (r, c,v)->r .= ∧(Tuple{0,1}, plot_mesh_dual, c, v), :type => InPlaceFunc())
-
-boundary = Examples.boundary_inds(Val{0}, plot_mesh)
-funcs[:∂C] = Dict(:operator => (∂ċ, ċ)->(∂ċ .= ċ; ∂ċ[boundary] .= 0), :type => InPlaceFunc())
-
+include("../../examples/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,y) -> begin
+      ∧(Tuple{0,1}, sd, x,y)
+    end
+    x => error("Unmatched operator $my_symbol")
+  end
+  return (args...) -> op(args...)
+end
 
 using Distributions
 c_dist = MvNormal([1, 5], [1.5, 1.5])
@@ -156,20 +158,24 @@ fig
 
 And the simulation result is then computed and visualized below.
 
-
 ```@example Debug
 using OrdinaryDiffEq
-func, code = gen_sim(explicit_ts, funcs, plot_mesh_dual; autodiff=false, params = [:V]);
+
+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
 
-prob = ODEProblem(func, c, (0.0, 100.0))
-sol = solve(prob, Tsit5(), p=v);
+u₀ = construct(PhysicsState, [VectorForm(c)], Float64[], [: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)[1:nv(plot_mesh)] for t in times]
+colors = [findnode(sol(t), :C) for t in times]
 
 # Initial frame
 fig, ax, ob = mesh(plot_mesh, color=colors[1], colorrange = extrema(vcat(colors...)))
@@ -179,104 +185,8 @@ framerate = 30
 
 # Animation
 record(fig, "diff_adv_right.gif", range(0.0, 100.0; length=150); framerate = 30) do t
-ob.color = sol(t)[1:nv(plot_mesh)]
+  ob.color = findnode(sol(t), :C)
 end
 ```
 
 ![](diff_adv_right.gif)
-
-# Debug the Simulation
-
-The benefit of having the computation graph which we generated above means that
-we have the opportunity to inspect the simulation at different stages of
-computation. Since each wire in the computation diagram has the data of a
-particular form on the mesh, this data can be visualized. The first step of
-getting at this data, though, is understanding the index associated with each
-wire. This key between indices and wires can be generated with the `sim_key`
-function.
-
-
-```@example Debug
-using Decapodes.Debug
-
-sim_key(explicit_ts)
-```
-
-Another key piece of data is what the initial state is at the time-step you are
-debugging. We choose the time step `t = 10` for this example, shown below.
-
-
-```@example Debug
-t = 10.0
-fig, ax, ob = mesh(plot_mesh; color=sol(t)[1:nv(plot_mesh_dual)])
-xlims!(ax, (0,10.0))
-fig
-```
-
-Now we can see what the value is of the result of the product between velocity
-and the concentration.
-
-
-```@example Debug
-fig, ax, ob = draw_wire(plot_mesh, plot_mesh_dual, explicit_ts, func, sol(10.0), 2;p = v, xlim=(0, 10.0), ylim=(0.0,10.0))
-fig
-```
-
-Now, this information doesn't seem all that useful at first glance. The result
-is basically a heatmap of the magnitude of the vector field across each edge of
-the mesh. Here, we see that this product has a greater magnitude on edges that
-are both aligned with the flow and have a high concentraiton at that point.
-
-Something which is a little more useful is to includ arrows! By using the
-`use_arrows` keyword argument, we are able to get arrows, directed along the
-mesh edge elements, which show the direction of the vectorfield at that point.
-The color of these arrows shows the magnitude of the vectorfield in that
-direction. Note that this method of plotting arrows with CairoMakie is fairly
-computationally expensive, so as the `n_arrows` argument increases, so will the
-time it takes to render.
-
-
-```@example Debug
-fig, ax, ob = draw_wire(plot_mesh, plot_mesh_dual, explicit_ts, func, sol(10.0), 2;p = v, xlim=(0, 10.0), ylim=(0.0,10.0), use_arrows=true, n_arrows=1200)
-fig
-```
-
-Note that all of the arrows from the result of the product are pointing in the
-direction of flow. This intuitively makes sense for how the concentration flow
-vectorfield should look.
-
-Below, we show the concentration flow resulting just from the gradient.
-
-
-```@example Debug
-fig, ax, ob = draw_wire(plot_mesh, plot_mesh_dual, explicit_ts, func, sol(10.0), 7;p = v, xlim=(0, 10.0), ylim=(0.0,10.0), use_arrows=true, n_arrows=1200)
-fig
-```
-
-Here, we see that the arrows are pointing away from areas of high concentration
-to areas of low concentration, and that the magnitude of the arrows grows
-greater as the slope of concentration is greater.
-
-Next, we can show what happens after the effects of diffusion and advection are
-added together.
-
-
-```@example Debug
-fig, ax, ob = draw_wire(plot_mesh, plot_mesh_dual, explicit_ts, func, sol(10.0), 4;p = v, xlim=(0, 10.0), ylim=(0.0,10.0), use_arrows=true, n_arrows=1200)
-fig
-```
-
-Finally, in order to compute the change of concentration, we visualize what the
-resulting change in concentration is.
-
-
-```@example Debug
-
-fig, ax, ob = draw_wire(plot_mesh, plot_mesh_dual, explicit_ts, func, sol(10.0), 3;p = v, xlim=(0, 10.0), ylim=(0.0,10.0), use_arrows=true, n_arrows=1200)
-fig
-```
-
-Here, we can see the positive change which is weighted heavily to the right of
-the existing distribution, and the slight negative change which follows from
-where the distribution is flowing.
-

docs/src/ice_dynamics.md

@@ -33,6 +33,14 @@ We use the `@decapode` macro to interpret the equations. Here, we have equation
 ```
 We'll change the term out front to Γ so we can demonstrate composition in a moment.
 
+In the exterior calculus, we could write the above equations like so:
+```math
+\partial_t(h) = \circ(\star, d, \star)(\Gamma\quad d(h)\quad \text{avg}_{01}|d(h)^\sharp|^{n-1} \quad \text{avg}_{01}(h^{n+2})).
+```
+
+`avg` here is an operator that performs the midpoint rule, setting the value at an edge to be the average of the values at its two vertices.
+
+
 ``` @example DEC
 halfar_eq2 = @decapode begin
   h::Form0