Finished up Brusselator doc page

Also filled out CISM page with some more sections.

docs/src/brussel/brussel.md

@@ -1,5 +1,7 @@
 # Brusselator
 
+This Brusselator example is adapted from [MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/tutorials/brusselator/#brusselator). The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a autocatalytic chemical reaction that takes place between two reactants `U` and `V`.
+
 ```@setup INFO
 include(joinpath(Base.@__DIR__, ".." , "..", "docinfo.jl"))
 info = DocInfo.Info()
@@ -9,6 +11,7 @@ info = DocInfo.Info()
 using CairoMakie
 import CairoMakie: wireframe, mesh, Figure, Axis
 
+using Catlab
 using CombinatorialSpaces
 using ComponentArrays
 using DiagrammaticEquations
@@ -20,11 +23,14 @@ using OrdinaryDiffEq
 using GeometryBasics: Point2, Point3
 Point2D = Point2{Float64}
 Point3D = Point3{Float64}
+nothing # hide
 ```
 
 ## The Model
+
+We establish the model for the Brusselator, with the two reactants `U` and `V`, modelled as residing on the vertices of the mesh. The equations encode a reaction that occurs independently at each point coupled with a diffusion term as well as a source term `F` in the case of `U`. Here `α` denotes the rate of diffusion for both reactants.
 ```@example DEC
-Brusselator = @decapode begin
+BrusselatorDynamics = @decapode begin
   (U, V)::Form0
   U2V::Form0
   (U̇, V̇)::Form0
@@ -39,16 +45,54 @@ Brusselator = @decapode begin
   ∂ₜ(U) == U̇
   ∂ₜ(V) == V̇
 end
+nothing # hide
+```
+
+## Boundary Conditions
+
+We now establish the Dirichlet boundary conditions for our model. Here we intend to set some portion of the `U` variable to be a fixed value on some portion of the mesh. At this point, the conditions are only set symbolically and their actual implementation can change. Note that these values are set at the beginning of execution, as shown by the computation graph.
+```@example DEC
+BrusselatorBoundaries = @decapode begin
+  B::Constant
+end
+
+BrusselatorMorphism = @relation () begin
+  rlb(C, Cb)
+end
+
+Brusselator = collate(
+  BrusselatorDynamics,
+  BrusselatorBoundaries,
+  BrusselatorMorphism,
+  Dict(
+    :C => :U,
+    :Cb => :B))
+
+to_graphviz(Brusselator)
 ```
 
 ## The Mesh
+
+We load our triangulated mesh with horizontal and vertical resolution being `h=0.01`. `Point3D` is being used for the primal mesh `s` for ease of visualization while `Point2D` is used for the dual mesh `sd` for better memory usage. Since this conversion only drops the final z-coordinate, no information is lost.
 ```@example DEC
-s = triangulated_grid(1,1,0.008,0.008,Point3D);
+h = 0.01
+s = triangulated_grid(1,1,h,h,Point3D);
 sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s);
 subdivide_duals!(sd, Circumcenter());
+
+fig = Figure() # hide
+ax = CairoMakie.Axis(fig[1,1], aspect=1) # hide
+wf = wireframe!(ax, s; linewidth=1) # hide
+save("Brusselator_rect.png", fig) # hide
+nothing # hide
 ```
 
+!["BrusselatorRect"](Brusselator_rect.png)
+
 ## Initial data
+We assign the initial values of `U` and `V` according to a continuous function. Since they both live on the vertices of our mesh, we can simply iterate over all point coordinates, extract the coordinate values (for either x or y) and compute the desired value. `F` has some logic attached to it encoding that it will "activate" only once the simulation has reached time `t = 1.1`.
+
+Here we also decide to set our boundary conditions to be `1.0` along the left and right sides of our mesh.
 ```@example DEC
 U = map(sd[:point]) do (_,y)
   22 * (y *(1-y))^(3/2)
@@ -58,6 +102,18 @@ V = map(sd[:point]) do (x,_)
   27 * (x *(1-x))^(3/2)
 end
 
+fig = Figure() # hide
+ax = CairoMakie.Axis(fig[1,1], aspect=1, title = "Initial value of U") # hide
+msh = CairoMakie.mesh!(ax, s, color=U, colormap=:jet, colorrange=extrema(U)) # hide
+Colorbar(fig[1,2], msh) # hide
+save("initial_U.png", fig) # hide
+
+fig = Figure() # hide
+ax = CairoMakie.Axis(fig[1,1], aspect=1, title = "Initial value of V") # hide
+msh = CairoMakie.mesh!(ax, s, color=V, colormap=:jet, colorrange=extrema(V)) # hide
+Colorbar(fig[1,2], msh) # hide
+save("initial_V.png", fig) # hide
+
 F₁ = map(sd[:point]) do (x,y)
  (x-0.3)^2 + (y-0.6)^2 ≤ (0.1)^2 ? 5.0 : 0.0
 end
@@ -66,44 +122,87 @@ F₂ = zeros(nv(sd))
 
 constants_and_parameters = (
   α = 0.001,
+  B = 1.0, # Boundary value
   F = t -> t ≥ 1.1 ? F₁ : F₂)
+nothing # hide
 ```
 
+![Initial U Conditions](initial_U.png)
+![Initial V Conditions](initial_V.png)
+
+```@example DEC
+# Find left and right vertices of mesh
+min_x = minimum(x -> x[1], s[:point])
+max_x = maximum(x -> x[1], s[:point])
+left_wall_idxs = findall(x -> x[1] == min_x, s[:point])
+right_wall_idxs = findall(x -> x[1] == max_x, s[:point])
+wall_idxs = vcat(left_wall_idxs, right_wall_idxs)
+
+function generate(sd, my_symbol; hodge=GeometricHodge())
+  op = @match my_symbol begin
+    :rlb => (x,y) -> begin
+      x[wall_idxs] .= y
+      x
+    end
+    _ => error("Unmatched operator $my_symbol")
+  end
+end
+
+fig = Figure() # hide
+ax = CairoMakie.Axis(fig[1,1], aspect=1, title = "Highlighted Boundary") # hide
+value = zeros(nv(sd)) # hide
+value[wall_idxs] .= 1.0 # hide
+msh = CairoMakie.mesh!(ax, s, color=value, colormap=:jet, colorrange=(0,2)) # hide
+Colorbar(fig[1,2], msh) # hide
+save("boundary.png", fig) # hide
+nothing # hide
+```
+
+![Boundary Visualized](boundary.png)
+
+
 ## Generate the Simulation
+
+We generate our simulation code and store the function in `fₘ` and then run our simulation for `t=11.5` simulated time units.
 ```@example DEC
 sim = evalsim(Brusselator)
-fₘ = sim(sd, nothing, DiagonalHodge())
+fₘ = sim(sd, generate, DiagonalHodge())
 
 u₀ = ComponentArray(U=U, V=V)
 
 tₑ = 11.5
 
 prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
 soln = solve(prob, Tsit5())
+soln.retcode
 ```
 
 ## Visualize
-```@example DEC
-
-fig = Figure();
-ax = Axis(fig[1,1])
-mesh!(ax, s, color=soln(tₑ).U, colormap=:jet)
 
+We can then use Makie to visualize the evolution of our Brusselator model.
+```@example DEC
 function save_dynamics(save_file_name)
   time = Observable(0.0)
   u = @lift(soln($time).U)
-  f = Figure()
-  ax = CairoMakie.Axis(f[1,1], title = @lift("Brusselator U Concentration at Time $($time)"))
-  gmsh = mesh!(ax, s, color=u, colormap=:jet,
-               colorrange=extrema(soln(tₑ).U))
-  Colorbar(f[1,2], gmsh)
+  v = @lift(soln($time).V)
+  f = Figure(; size = (500, 850))
+  ax_U = CairoMakie.Axis(f[1,1], title = @lift("Concentration of U at Time $($time)"))
+  ax_V = CairoMakie.Axis(f[2,1], title = @lift("Concentration of V at Time $($time)"))
+
+  msh_U = mesh!(ax_U, s, color=u, colormap=:jet, colorrange=extrema(soln(tₑ).U))
+  Colorbar(f[1,2], msh_U)
+
+  msh_V = mesh!(ax_V, s, color=v, colormap=:jet, colorrange=extrema(soln(tₑ).V))
+  Colorbar(f[2,2], msh_V)
+
   timestamps = range(0, tₑ, step=1e-1)
   record(f, save_file_name, timestamps; framerate = 15) do t
     time[] = t
   end
 end
 
 save_dynamics("brusselator.gif")
+nothing # hide
 ```
 
 ![Brusselator_results_flat](brusselator.gif)

docs/src/cism/cism.md

@@ -194,11 +194,8 @@ constants_and_parameters = (
 
 We provide here the mapping from symbols to differential operators. As more of the differential operators of the DEC are implemented, they are upstreamed to the Decapodes and CombinatorialSpaces libraries. Of course, users can also provide their own implementations of these operators and others as they see fit.
 
+## Setting up the Simulation
 ```@example DEC
-#############################################
-# Define how symbols map to Julia functions #
-#############################################
-
 function generate(sd, my_symbol; hodge=GeometricHodge())
   # We pre-allocate matrices that encode differential operators.
   op = @match my_symbol begin
@@ -216,16 +213,14 @@ end
 The `gensim` function takes our high-level representation of the physics equations and produces compiled simulation code. It performs optimizations such as allocating memory for intermediate variables, and so on.
 
 ```@example DEC
-#######################
-# Generate simulation #
-#######################
-
 sim = eval(gensim(ice_dynamics2D))
 fₘ = sim(sd, generate)
 ```
 
 Julia is a "Just-In-Time" compiled language. That means that functions are compiled the first time they are called, and later calls to those functions skip this step. To get a feel for just how fast this simulation is, we will run the dynamics twice, once for a very short timespan to trigger pre-compilation, and then again for the actual dynamics.
 
+## Running the Simulation
+
 ```@example DEC
 # Pre-compile simulation
 
@@ -264,6 +259,8 @@ Here we save the solution information to a [file](ice_dynamics2D.jld2).
 @save "ice_dynamics2D.jld2" soln
 ```
 
+## Result Comparison and Analysis
+
 We recall that these dynamics are of the "shallow slope" and "shallow ice" approximations. So, at the edge of our parabolic dome of ice, we expect increased error as the slope increases. On the interior of the dome, we expect the dynamics to match more closely that given by the analytic model. We will see that the CISM results likewise accumulate error in the same neighborhood.
 
 !["Halfar Small Ice Approximation Quote"](halfar_quote.png)

docs/src/faq/faq.md

@@ -1,18 +1,22 @@
 # FAQ
 
-## 1. How to incorporate scalar or vector field input data where you have a function of the embedded coordinates?
+## 1. How do I incorporate scalar or vector field input data where you have a function of the embedded coordinates?
 
 We can take a look at the [Brusselator page](../brussel/brussel.md#initial-data) which sets the values of each point on its mesh to a value as determined by some function. This can be done in a similar manner in 1D and 2D. 
 
 The Brusselator also demonstrates, with the variable `F`, how one can go about changing the function by which these values are set over time.
 
-## 2. How to incorporate input data from a file with linear interpolation?
+## 2. How do I incorporate input data from a file with linear interpolation?
 
 The Grigoriev Ice Cap model has a section where after the initial data is loaded from a TIF, the data is interpolated so that it may fit over a discrete mesh of our choosing. The link for that is [here](../grigoriev/grigoriev.md#loading-a-scientific-dataset).
 
-## 3. How to set boundary conditions like fixed value, no flux, and no slip?
+## 3. How do I set boundary conditions like fixed value, no-flux, and no-slip?
 
-## 4. How to plot a derived quantity?
+Boundary conditions can be set by using "collages", which can take two variables among two different Decapodes and apply a function on the first. A general workflow would be to have the first Decapode encode the physics and have the second one encode values for boundary conditions. They can be related by a function that will mask the first variable and replace the desired values with second. An example of applying fixed boundary conditions would be in the [Brusselator page](../brussel/brussel.md#boundary-conditions).
+
+A similar workflow can be used for "no-flux" and "no-slip" conditions by fixing the value of the appropriate variable to be 0. 
+
+## 4. How do I plot derived quantities?
 
 Plotting in DECAPODES is commonly done with the [Makie](https://github.com/MakieOrg/Makie.jl) package in Julia. Makie allows for creating both still images, which are useful for visualizing the mesh itself and initial/final conditions, and videos, which can capture the full simulation from start to end.
 
@@ -23,9 +27,9 @@ Plotting in DECAPODES is commonly done with the [Makie](https://github.com/Makie
 - For [3D visualization](../ice_dynamics/ice_dynamics.md#2-manifold-in-3d)
 
 
-## 5. How to add artificial diffusion for 0- or 1-forms?
+## 5. How do I add artificial diffusion for 0- or 1-forms?
 
-Without viscosity - i.e. when ``\mu = 0`` - the incompressible (inviscid) Navier-Stokes equations can be formulated like so:
+Without viscosity - i.e. when ``μ = 0`` - the incompressible (inviscid) Navier-Stokes equations can be formulated like so:
 
 ```julia
 eq11_inviscid_poisson = @decapode begin
@@ -40,7 +44,7 @@ eq11_inviscid_poisson = @decapode begin
 end
 ```
 
-Adding a viscosity term can be accomplished by simply added the appropriate term, and declaring the ``\mu`` constant:
+Adding a viscosity term can be accomplished by simply added the appropriate term, and declaring the ``μ`` constant:
 
 ```julia
 eq11_viscid_poisson = @decapode begin
@@ -58,21 +62,25 @@ end
 
 More demonstrations on how to iterate between formulations of the same physics (the incompressible Navier-Stokes equations) is available in further detail on the [Vortices](../navier_stokes/ns.md) docs page and in the script available there.
 
-## 6. How to use a Laplacian solver / multigrid?
+## 6. How do I use multigrid methods?
 
 To use multigrid methods in the Laplacian solver, you need to create a `PrimalGeometricMapSeries` that will take a coarse mesh and apply a subdivision method to it some number of times. After that, just use this result as you would a regular mesh for simulation.
 
 ```julia
-  s = triangulated_grid(100,100,1,1,Point3D)
+s = triangulated_grid(100,100,10,10,Point3D)
+
+# Binary subdivide 4 times
+series = PrimalGeometricMapSeries(s, binary_subdivision_map, 4);
 
-  series = PrimalGeometricMapSeries(s, binary_subdivision_map, 4);
+# Retrieve highest resolution mesh
+sd = finest_mesh(series)
 
   ...
 
-  f_mg = sim_mg(series, generate);
+f_mg = sim_mg(series, generate);
 ```
 
-## 7. How to do a bunch of workflows?
+## 7. What are general workflows for DECAPODES?
 
 A common workflow is to iterate through multiple different models as is done in the [Vorticity Model page](../navier_stokes/ns.md). A formulation is first done with a direct vorticity formulation but a quick run finds that this setup is unstable. A second formulation introduces a Laplacian solve which produces nice results.