added Fokker-Planck; Gray-Scott to literate

examples/chemistry/gray_scott.jl → docs/literate/chemistry/gray_scott.jl

@@ -102,13 +102,13 @@ constants_and_parameters = (
 tₑ = 10_000
 
 @info("Solving")
-prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
-soln = solve(prob, Tsit5())
+problem = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
+solution = solve(problem, Tsit5())
 @info("Done")
 
 function save_dynamics(save_file_name)
   time = Observable(0.0)
-  u = @lift(soln($time).U)
+  u = @lift(solution($time).U)
   f = Figure()
   ax_U = CairoMakie.Axis(f[1,1], title = @lift("Concentration of U at Time $($time)"))
 
@@ -121,4 +121,4 @@ function save_dynamics(save_file_name)
   end
 end
 
-save_dynamics("gs_f=$(f)_k=$(k).mp4")
+save_dynamics("gs_f=$(f)_k=$(k).mp4")

docs/make.jl

@@ -60,7 +60,7 @@ makedocs(
     "Concepts" => Any[
         "Equations" => "concepts/equations.md",
         "Meshes" => "concepts/meshes.md",
-        "Composition" => "concepts/composition.md",
+        "Custom Operators" => "concepts/generate.md",
     ],
     "Zoo" => Any[
         "Vortices" => "navier_stokes/ns.md",
@@ -71,14 +71,18 @@ makedocs(
         "CISM v2.1" => "cism/cism.md",
         "Grigoriev Ice Cap" => "grigoriev/grigoriev.md", # Requires ice_dynamics
         "Budyko-Sellers-Halfar" => "bsh/budyko_sellers_halfar.md", # Requires ice_dynamics
-        "Halfar-EBM-Water" => "ebm_melt/ebm_melt.md",
+        # "Halfar-EBM-Water" => "ebm_melt/ebm_melt.md",
         "Halfar-NS" => "halmo/halmo.md", # Requires grigoriev
         "NHS" => "nhs/nhs_lite.md",
         "Pipe Flow" => "poiseuille/poiseuille.md",
+        "Fokker-Planck" => "fokker_planck/fokker_planck.md"
     ],
     "Examples" => Any[
-        "Oncology" => "examples/oncology/tumor_proliferation_invasion.md"
-        "MHD" => "examples/mhd.md" # TODO convert original file to a docs page
+        "Heat" => "examples/diff_adv/heat.md",
+        "Burger" => "examples/diff_adv/burger.md",
+        "Gray-Scott" => "examples/chemistry/gray_scott.md",
+        "Oncology" => "examples/oncology/tumor_proliferation_invasion.md",
+        "MHD" => "examples/mhd.md", # TODO convert original file to a docs page
     ],
     "Misc Features" => "bc/bc_debug.md", # Requires overview
     "FAQ" => "faq/faq.md",

docs/src/cism/cism.md

@@ -227,9 +227,9 @@ Julia is a "Just-In-Time" compiled language. That means that functions are compi
 # Julia will pre-compile the generated simulation the first time it is run.
 @info("Precompiling Solver")
 # We run for a short timespan to pre-compile.
-prob = ODEProblem(fₘ, u₀, (0, 1e-8), constants_and_parameters)
-soln = solve(prob, Tsit5())
-soln.retcode != :Unstable || error("Solver was not stable")
+problem = ODEProblem(fₘ, u₀, (0, 1e-8), constants_and_parameters)
+solution = solve(problem, Tsit5())
+solution.retcode != :Unstable || error("Solver was not stable")
 ```
 
 ```@example DEC
@@ -238,9 +238,9 @@ tₑ = 200
 
 # This next run should be fast.
 @info("Solving")
-prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
-soln = solve(prob, Tsit5(), saveat=0.1)
-@show soln.retcode
+problem = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
+solution = solve(problem, Tsit5(), saveat=0.1)
+@show solution.retcode
 @info("Done")
 ```
 
@@ -249,14 +249,14 @@ We can benchmark the compiled simulation with `@benchmarkable`. This macro runs
 ```@example DEC
 # Time the simulation
 
-b = @benchmarkable solve(prob, Tsit5(), saveat=0.1)
+b = @benchmarkable solve(problem, Tsit5(), saveat=0.1)
 c = run(b)
 ```
 
 Here we save the solution information to a [file](ice_dynamics2D.jld2).
 
 ```@example DEC
-@save "ice_dynamics2D.jld2" soln
+@save "ice_dynamics2D.jld2" solution
 ```
 
 ## Result Comparison and Analysis
@@ -272,7 +272,7 @@ function plot_final_conditions()
   ax = CairoMakie.Axis(fig[1,1],
     title="Modeled thickness (m) at time 200.0",
     aspect=0.6, xticks = [0, 3e4, 6e4])
-  msh = mesh!(ax, s, color=soln(200.0).dynamics_h, colormap=:jet)
+  msh = mesh!(ax, s, color=solution(200.0).dynamics_h, colormap=:jet)
   Colorbar(fig[1,2], msh)
   fig
 end
@@ -306,7 +306,7 @@ nothing # hide
 # Plot the error.
 function plot_error()
   hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s))
-  h_diff = soln(tₑ).dynamics_h - hₐ
+  h_diff = solution(tₑ).dynamics_h - hₐ
   extrema(h_diff)
   fig = Figure()
   ax = CairoMakie.Axis(fig[1,1],
@@ -328,7 +328,7 @@ We compute below that the maximum absolute error is approximately 89 meters. We
 ```@example DEC
 # Compute max absolute error:
 hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s))
-h_diff = soln(tₑ).dynamics_h - hₐ
+h_diff = solution(tₑ).dynamics_h - hₐ
 maximum(abs.(h_diff))
 ```
 
@@ -338,14 +338,14 @@ We compute root-mean-square (RMS) error as well, both over the entire domain, an
 # Compute RMS not considering the "outside".
 hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s))
 nonzeros = findall(!=(0), hₐ)
-h_diff = soln(tₑ).dynamics_h - hₐ
+h_diff = solution(tₑ).dynamics_h - hₐ
 rmse = sqrt(sum(map(x -> x*x, h_diff[nonzeros])) / length(nonzeros))
 ```
 
 ```@example DEC
 # Compute RMS of the entire domain.
 hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s))
-h_diff = soln(tₑ).dynamics_h - hₐ
+h_diff = solution(tₑ).dynamics_h - hₐ
 rmse = sqrt(sum(map(x -> x*x, h_diff)) / length(h_diff))
 ```
 
@@ -355,10 +355,10 @@ begin
   frames = 100
   fig = Figure()
   ax = CairoMakie.Axis(fig[1,1], aspect=0.6, xticks = [0, 3e4, 6e4])
-  msh = mesh!(ax, s, color=soln(0).dynamics_h, colormap=:jet, colorrange=extrema(soln(tₑ).dynamics_h))
+  msh = mesh!(ax, s, color=solution(0).dynamics_h, colormap=:jet, colorrange=extrema(solution(tₑ).dynamics_h))
   Colorbar(fig[1,2], msh)
   record(fig, "ice_dynamics_cism.gif", range(0.0, tₑ; length=frames); framerate = 30) do t
-    msh.color = soln(t).dynamics_h
+    msh.color = solution(t).dynamics_h
   end
 end
 ```

docs/src/concepts/generate.md

@@ -1,3 +1,7 @@
+``` @example DEC
+using Decapodes
+```
+
 ## Custom Operators
 
 Decapodes.jl already defines a suite of operators from the Discrete Exterior Calculus. However it is often the case that an implementation requires custom operators. Sometimes, this is just matter of building operators through composition. However Decapodes accepts a lookup table of functions which are included when parsing a Decapodes expression. This allows for new operators with their own symbols to be defined.
@@ -14,7 +18,7 @@ Decapodes uses the Discrete Exterior Calculus to discretize our differential ope
 
 If this code seems too low level, do not worry. Decapodes defines and caches for you many differential operators behind the scenes, so you do not have to worry about defining your own.
 
-```@example DEC
+```
 lap_mat = dec_hodge_star(1,dualmesh) * dec_differential(0,dualmesh) * dec_inv_hodge_star(0,dualmesh) * dec_dual_derivative(0,dualmesh)
 
 function generate(dualmesh, my_symbol; hodge=DiagonalHodge())

docs/src/ebm_melt/ebm_melt.md

@@ -15,6 +15,7 @@ using DiagrammaticEquations
 using Decapodes
 
 # External Dependencies
+using Artifacts
 using ComponentArrays
 using CoordRefSystems
 using CairoMakie
@@ -41,8 +42,7 @@ wireframe(s_plots)
 The data provided by the Polar Science Center is given as a NetCDF file. Ice thickness is a matrix with the same dimensions as a matrix provided Latitude and Longitude of the associated point on the Earth's surface. We need to convert between polar and Cartesian coordinates to use this data on our mesh.
 
 ``` @example DEC
-ice_thickness_file = "piomas20c.heff.1901.2010.v1.0.nc"
-run(`curl -o $ice_thickness_file https://pscfiles.apl.uw.edu/axel/piomas20c/v1.0/monthly/piomas20c.heff.1901.2010.v1.0.nc`)
+ice_thickness_file = rootpath"piomas20c.heff.1901.2010.v1.0.nc"
 
 # Use ncinfo(ice_thickness_file) to interactively get information on variables.
 # Sea ice thickness ("sit") has dimensions of [y, x, time].

docs/src/fokker_planck/fokker_planck.md

@@ -0,0 +1,88 @@
+# Fokker-Planck
+
+``` @example DEC
+using Catlab, CombinatorialSpaces, Decapodes, DiagrammaticEquations
+using CairoMakie, ComponentArrays, LinearAlgebra, MLStyle, ComponentArrays
+using OrdinaryDiffEq
+using GeometryBasics: Point3
+Point3D = Point3{Float64}
+using Arpack
+```
+
+Let's specify physics
+``` @example DEC
+Fokker_Planck = @decapode begin
+  (ρ,Ψ)::Form0
+  β⁻¹::Constant
+  ∂ₜ(ρ) == ∘(⋆,d,⋆)(d(Ψ)∧ρ) + β⁻¹*Δ(ρ)
+end
+```
+
+Specify the domain
+``` @example DEC
+spheremesh = loadmesh(Icosphere(6))
+dualmesh = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(spheremesh);
+subdivide_duals!(dualmesh, Barycenter())
+```
+
+Compile the simulation
+``` @example DEC
+simulation = eval(gensim(Fokker_Planck))
+f = simulation(dualmesh, nothing)
+```
+
+Specify initial conditions. Ψ must be a smooth function. Choose an interesting eigenfunction. We require that ρ integrated over the surface is 1, since it is a PDF. On a sphere where ρ(x,y,z) is proportional to the x-coordinate, that means divide by 2π.
+``` @example DEC
+Δ0 = Δ(0, dualmesh)
+Ψ = real.(eigs(Δ0, nev=32, which=:LR)[2][:,32])
+ρ = map(point(dualmesh)) do (x,y,z)
+  abs(x)
+end / 2π
+```
+
+Let's define the structures which hold the constants and state variables for the
+simulation, respectively.
+``` @example DEC
+constants_and_parameters = (β⁻¹ = 1e-2,)
+u0 = ComponentArray(Ψ=Ψ, ρ=ρ)
+```
+
+Run the simulation.
+``` @example DEC
+tₑ= 20.0
+problem = ODEProblem(f, u0, (0, tₑ), constants_and_parameters);
+solution = solve(problem, Tsit5(), progress=true, progress_steps=1);
+```
+
+Verify that the probability distribbution function is still a probability distribution. We'll show that the sum of the values on the
+dual 2-form integrate (sum to) unity,
+``` @example DEC
+s0 = dec_hodge_star(0, dualmesh)
+@info sum(s0 * solution(tₑ).ρ)
+@info any(solution(tₑ).ρ .≤ 0)
+```
+
+Now we will create a GIF.
+
+``` @example DEC
+function save_gif(file_name, soln)
+  time = Observable(0.0)
+  fig = Figure()
+  Label(fig[1, 1, Top()], @lift("ρ at $($time)"), padding = (0, 0, 5, 0))
+  ax = LScene(fig[1,1], scenekw=(lights=[],))
+  msh = CairoMakie.mesh!(ax, spheremesh,
+    color=@lift(soln($time).ρ),
+    colorrange=(0,1),
+    colormap=:jet)
+
+  Colorbar(fig[1,2], msh)
+  frames = range(0.0, tₑ; length=21)
+  record(fig, file_name, frames; framerate = 10) do t
+    time[] = t
+  end
+end
+gif = save_gif("fokker_planck.gif", solution)
+```
+
+!["FokkerPlanck"](fokker_planck.gif)
+

docs/src/klausmeier/klausmeier.md

@@ -127,6 +127,19 @@ Let's pass our mesh and methods of generating operators to our simulation code.
 
 ```@example DEC
 # Instantiate Simulation
+
+lap_mat = dec_hodge_star(1,dualmesh) * dec_differential(0,dualmesh) * dec_inv_hodge_star(0,dualmesh) * dec_dual_derivative(0,dualmesh)
+
+function generate(sd, my_symbol; hodge=DiagonalHodge())
+  op = @match my_symbol begin
+    :Δ => x -> begin
+      lap_mat * x
+    end
+  end
+  return (args...) -> op(args...)
+end
+
+
 fₘ = sim(dualmesh, generate, DiagonalHodge())
 ```
 
@@ -156,31 +169,31 @@ Let's execute our simulation.
 ```@example DEC
 # Run Simulation
 tₑ = 300.0
-prob = ODEProblem(fₘ, u₀, (0.0, tₑ), cs_ps)
-soln = solve(prob, Tsit5(), saveat=0.1, save_idxs=[:N, :W])
-soln.retcode
+problem = ODEProblem(fₘ, u₀, (0.0, tₑ), cs_ps)
+solution = solve(problem, Tsit5(), saveat=0.1, save_idxs=[:N, :W])
+solution.retcode
 ```
 
 ## Animation
 
 Let's perform some basic visualization and analysis of our results to verify our dynamics.
 
 ```@example DEC
-n = soln(0).N
-nₑ = soln(tₑ).N
-w = soln(0).W
-wₑ = soln(tₑ).W
+n = solution(0).N
+nₑ = solution(tₑ).N
+w = solution(0).W
+wₑ = solution(tₑ).W
 nothing # hide
 ```
 
 ```@example DEC
 # Animate dynamics
 function save_dynamics(form_name, framerate, filename)
   time = Observable(0.0)
-  ys = @lift(getproperty(soln($time), form_name))
+  ys = @lift(getproperty(solution($time), form_name))
   xcoords = [0, accumulate(+, dualmesh[:length])[1:end-1]...]
   fig = lines(xcoords, ys, color=:green, linewidth=4.0,
-    colorrange=extrema(getproperty(soln(0), form_name));
+    colorrange=extrema(getproperty(solution(0), form_name));
     axis = (; title = @lift("Klausmeier $(String(form_name)) at $($time)")))
   timestamps = range(0, tₑ, step=1)
   record(fig, filename, timestamps; framerate=framerate) do t