Add FAQ page to Decapodes

docs/make.jl

@@ -52,6 +52,7 @@ makedocs(
     "Meshes" => "meshes/meshes.md",
     "Vortices" => "navier_stokes/ns.md",
     "Harmonics" => "harmonics/harmonics.md",
+    "Brusselator" => "brussel/brussel.md",
     "Cahn-Hilliard" => "ch/cahn-hilliard.md",
     "Klausmeier" => "klausmeier/klausmeier.md",
     "CISM v2.1" => "cism/cism.md",
@@ -63,6 +64,7 @@ makedocs(
     "NHS" => "nhs/nhs_lite.md",
     "Pipe Flow" => "poiseuille/poiseuille.md",
     "Misc Features" => "bc/bc_debug.md", # Requires overview
+    "FAQ" => "faq/faq.md",
     "ASCII Operators" => "ascii.md",
     "Canonical Models" => "canon.md",
     "Library Reference" => "api.md"

docs/src/brussel/brussel.md

@@ -0,0 +1,213 @@
+# Brusselator
+
+This Brusselator example is adapted from [MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/tutorials/brusselator/#brusselator)'s page on the same topic. 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()
+```
+## Dependencies
+```@example DEC
+using CairoMakie
+import CairoMakie: wireframe, mesh, Figure, Axis
+
+using Catlab
+using CombinatorialSpaces
+using ComponentArrays
+using DiagrammaticEquations
+using Decapodes
+using LinearAlgebra
+using MLStyle
+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`, modeled 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
+BrusselatorDynamics = @decapode begin
+  (U, V)::Form0
+  U2V::Form0
+  (U̇, V̇)::Form0
+
+  (α)::Constant
+  F::Parameter
+
+  U2V == (U .* U) .* V
+
+  U̇ == 1 + U2V - (4.4 * U) + (α * Δ(U)) + F
+  V̇ == (3.4 * U) - U2V + (α * Δ(V))
+  ∂ₜ(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 boundary 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
+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)
+end
+
+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
+
+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, 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
+
+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)
+  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)
+
+```@example INFO
+DocInfo.get_report(info) # hide
+```
+

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

@@ -0,0 +1,89 @@
+# FAQ
+
+## 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 both 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 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 and the data is interpolated so that it may fit over a discrete mesh of our choosing. You may view that [here](../grigoriev/grigoriev.md#loading-a-scientific-dataset).
+
+## 3. How do I set boundary conditions like fixed value, no-flux, and no-slip?
+
+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 a second Decapode encode values for the boundary conditions. They can be related by a function that will mask the first variable and replace the desired values with those of the second's. 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.
+
+- For [1D visualization](../ice_dynamics/ice_dynamics.md#visualize)
+
+- For [2D visualization](../ice_dynamics/ice_dynamics.md#visualize-2d)
+
+- For [3D visualization](../ice_dynamics/ice_dynamics.md#2-manifold-in-3d)
+
+
+## 5. How do I add artificial diffusion for 0- or 1-forms?
+
+Without viscosity - i.e. when ``μ = 0`` - the incompressible (inviscid) Navier-Stokes equations can be formulated like so:
+
+```julia
+eq11_inviscid_poisson = @decapode begin
+  d𝐮::DualForm2
+  𝐮::DualForm1
+  ψ::Form0
+
+  ψ == Δ⁻¹(⋆(d𝐮))
+  𝐮 == ⋆(d(ψ))
+
+  ∂ₜ(d𝐮) ==  (-1) * ∘(♭♯, ⋆₁, d̃₁)(∧ᵈᵖ₁₀(𝐮, ⋆(d𝐮)))
+end
+```
+
+Adding a viscosity term can be accomplished by simply added the appropriate term, and declaring the ``μ`` constant:
+
+```julia
+eq11_viscid_poisson = @decapode begin
+  d𝐮::DualForm2
+  𝐮::DualForm1
+  ψ::Form0
+  μ::Constant
+
+  ψ == Δ⁻¹(⋆(d𝐮))
+  𝐮 == ⋆(d(ψ))
+
+  ∂ₜ(d𝐮) ==  μ * ∘(⋆, d, ⋆, d)(d𝐮) + (-1) * ∘(♭♯, ⋆₁, d̃₁)(∧ᵈᵖ₁₀(𝐮, ⋆(d𝐮)))
+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 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,10,10,Point3D)
+
+# Binary subdivide 4 times
+series = PrimalGeometricMapSeries(s, binary_subdivision_map, 4);
+
+# Retrieve highest resolution mesh
+sd = finest_mesh(series)
+
+  ...
+
+f_mg = sim_mg(series, generate);
+```
+
+## 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.
+
+Similar workflows may retain the same model but may iterate on the types of meshes/initial conditions used. An excellent example of this is found in the [Glacial Flow page](../ice_dynamics/ice_dynamics.md) where the model is first run in a [1D](../ice_dynamics/ice_dynamics.md#Define-a-mesh) setting and then quickly promoted to both [2D](../ice_dynamics/ice_dynamics.md#Define-our-mesh) and [3D](../ice_dynamics/ice_dynamics.md#2-Manifold-in-3D). This allows either running some dynamics in a more complicated setting, as just discussed, or allows for simplifying higher dimensional models by some sort of symmetry.