From 408c66dc6d5fb3512ee53a40d21fcd4e6983ebda Mon Sep 17 00:00:00 2001
From: GeorgeR227 <78235421+GeorgeR227@users.noreply.github.com>
Date: Thu, 19 Dec 2024 17:33:17 -0500
Subject: [PATCH 1/2] Update gray_scott.jl

Updates the gray scott model with boundary conditions and constants to produce cool looking behavior. This should probably be turned into a doc page.
---
 examples/chemistry/gray_scott.jl | 232 +++++++++----------------------
 1 file changed, 69 insertions(+), 163 deletions(-)

diff --git a/examples/chemistry/gray_scott.jl b/examples/chemistry/gray_scott.jl
index bdafd2ba..8add5b8b 100644
--- a/examples/chemistry/gray_scott.jl
+++ b/examples/chemistry/gray_scott.jl
@@ -1,18 +1,11 @@
-using Catlab
-using Catlab.Graphics
 using CombinatorialSpaces
-using CombinatorialSpaces.ExteriorCalculus
 using DiagrammaticEquations
-using DiagrammaticEquations.Deca
 using Decapodes
 using MLStyle
 using OrdinaryDiffEq
 using LinearAlgebra
 using CairoMakie
 import CairoMakie: wireframe, mesh, Figure, Axis
-using Logging
-using JLD2
-using Printf
 using ComponentArrays
 
 using GeometryBasics: Point2, Point3
@@ -27,189 +20,102 @@ GrayScott = @decapode begin
   (UV2)::Form0
   (U̇, V̇)::Form0
   (f, k, rᵤ, rᵥ)::Constant
+  B::Constant
+
   UV2 == (U .* (V .* V))
-  U̇ == rᵤ * Δ(U) - UV2 + f * (1 .- U)
-  V̇ == rᵥ * Δ(V) + UV2 - (f + k) .* V
+  lap_U == mask(Δ(U), B)
+  lap_V == mask(Δ(V), B)
+
+  U̇ == rᵤ * lap_U - UV2 + f * (1 .- U)
+  V̇ == rᵥ * lap_V + UV2 - (f + k) .* V
   ∂ₜ(U) == U̇
   ∂ₜ(V) == V̇
 end
 
-# Visualize. You must have graphviz installed.
-to_graphviz(GrayScott)
-
-# We resolve types of intermediate variables using sets of rules.
-infer_types!(GrayScott)
-to_graphviz(GrayScott)
-
-# Resolve overloads. i.e. ~dispatch
-resolve_overloads!(GrayScott)
-to_graphviz(GrayScott)
-
-s = loadmesh(Rectangle_30x10())
-scaling_mat = Diagonal([1/maximum(x->x[1], s[:point]),
-                        1/maximum(x->x[2], s[:point]),
-                        1.0])
-s[:point] = map(x -> scaling_mat*x, s[:point])
-s[:edge_orientation] = false
-orient!(s)
-# Visualize the mesh.
-wireframe(s)
-sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s)
-subdivide_duals!(sd, Circumcenter())
-
-# Define how operations map to Julia functions.
-function generate(sd, my_symbol; hodge=GeometricHodge()) end
-
-# Create initial data.
-@assert all(map(sd[:point]) do (x,y)
-  0.0 ≤ x ≤ 1.0 && 0.0 ≤ y ≤ 1.0
-end)
-
-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
-
-constants_and_parameters = (
-  f = 0.024,
-  k = 0.055,
-  rᵤ = 0.01,
-  rᵥ = 0.005)
-
-# Generate the simulation.
-gensim(expand_operators(GrayScott))
-sim = eval(gensim(expand_operators(GrayScott)))
-fₘ = sim(sd, generate)
-
-# Create problem and run sim for t ∈ [0,tₑ).
-# Map symbols to data.
-u₀ = ComponentArray(U=U,V=V)
+n = 100
 
-# Visualize the initial conditions.
-# If GLMakie throws errors, then update your graphics drivers,
-# or use an alternative Makie backend like CairoMakie.
-fig_ic = Figure()
-p1 = mesh(fig_ic[1,2], s, color=u₀.U, colormap=:jet)
-p2 = mesh(fig_ic[1,3], s, color=u₀.V, colormap=:jet)
-display(fig_ic)
+s = triangulated_grid(n,n,0.8,0.8,Point3D);
+sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s);
+subdivide_duals!(sd, Circumcenter());
 
-tₑ = 11.5
+sim = eval(gensim(GrayScott))
 
-@info("Precompiling Solver")
-prob = ODEProblem(fₘ, u₀, (0, 1e-4), constants_and_parameters)
-soln = solve(prob, Tsit5())
-soln.retcode != :Unstable || error("Solver was not stable")
-@info("Solving")
-prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
-soln = solve(prob, Tsit5())
-@info("Done")
+left_wall_idxs = findall(x -> x[1] <= 1.0, s[:point])
+right_wall_idxs = findall(x -> x[1] >= n - 1.0, s[:point])
+top_wall_idxs = findall(y -> y[2] == 0.0, s[:point])
+bot_wall_idxs = findall(y -> y[2] == n, s[:point])
 
-@save "gray_scott.jld2" soln
-
-# Visualize the final conditions.
-mesh(s, color=soln(tₑ).U, colormap=:jet)
-
-begin # BEGIN Gif creation
-frames = 100
-## Initial frame
-fig = Figure(resolution = (1200, 800))
-p1 = mesh(fig[1,2], s, color=soln(0).U, colormap=:jet, colorrange=extrema(soln(0).U))
-p2 = mesh(fig[1,4], s, color=soln(0).V, colormap=:jet, colorrange=extrema(soln(0).V))
-ax1 = Axis(fig[1,2], width = 400, height = 400)
-ax2 = Axis(fig[1,4], width = 400, height = 400)
-hidedecorations!(ax1)
-hidedecorations!(ax2)
-hidespines!(ax1)
-hidespines!(ax2)
-Colorbar(fig[1,1], colormap=:jet, colorrange=extrema(soln(0).U))
-Colorbar(fig[1,5], colormap=:jet, colorrange=extrema(soln(0).V))
-Label(fig[1,2,Top()], "U")
-Label(fig[1,4,Top()], "V")
-lab1 = Label(fig[1,3], "")
-
-## Animation
-record(fig, "gray_scott.gif", range(0.0, tₑ; length=frames); framerate = 15) do t
-    p1.plot.color = soln(t).U
-    p2.plot.color = soln(t).V
-    lab1.text = @sprintf("%.2f",t)
+wall_idxs = unique(vcat(left_wall_idxs, right_wall_idxs, top_wall_idxs, bot_wall_idxs))
+function generate(sd, my_symbol; hodge=GeometricHodge())
+  op = @match my_symbol begin
+    :mask => (x,y) -> begin
+      x[wall_idxs] .= y
+      x
+    end
+    _ => error("Unmatched operator $my_symbol")
+  end
 end
 
-end # END Gif creation
+fₘ = sim(sd, generate, DiagonalHodge())
 
-# Run on the sphere.
-# You can use lower resolution meshes, such as Icosphere(3).
-s = loadmesh(Icosphere(5))
-orient!(s)
-# Visualize the mesh.
-wireframe(s)
-sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s)
-subdivide_duals!(sd, Circumcenter())
+init_multi = 0.5
 
-# Define how operations map to Julia functions.
-function generate(sd, my_symbol; hodge=GeometricHodge()) end
+U = rand(0.0:0.001:0.1, nv(sd))
+V = zeros(nv(sd))
 
-# Create initial data.
-U = map(sd[:point]) do (_,y,_)
-  abs(y)
-end
+mid = div(n, 2)
 
-V = map(sd[:point]) do (x,_,_)
-  abs(x)
-end
+mid_p = Point2D(mid, mid)
 
-constants_and_parameters = (
-  f = 0.024,
-  k = 0.055,
-  rᵤ = 0.01,
-  rᵥ = 0.005)
+init = map(p -> if norm(p - mid_p, Inf) <= 5; 1.0 .* init_multi; else 0.0; end, sd[:point])
 
-# Generate the simulation.
-fₘ = sim(sd, generate)
+# Set up an initial small disturbance
+U .+= init
+V .+= 0.5 * init
 
-# Create problem and run sim for t ∈ [0,tₑ).
-# Map symbols to data.
 u₀ = ComponentArray(U=U,V=V)
 
-# Visualize the initial conditions.
-# If GLMakie throws errors, then update your graphics drivers,
-# or use an alternative Makie backend like CairoMakie.
-fig_ic = Figure()
-p1 = mesh(fig_ic[1,2], s, color=u₀.U, colormap=:jet)
-p2 = mesh(fig_ic[1,3], s, color=u₀.V, colormap=:jet)
-display(fig_ic)
-
-tₑ = 11.5
+f = 0.055
+k = 0.062
+constants_and_parameters = (
+  rᵤ = 0.16,
+  rᵥ = 0.08,
+  f = f,
+  k = k,
+  B = 0)
+
+# fig = Figure();
+# ax = CairoMakie.Axis(fig[1,1], aspect=1, title = "Initial value of U")
+# msh = CairoMakie.mesh!(ax, s, color=U, colormap=:jet, colorrange=(extrema(U)))
+# Colorbar(fig[1,2], msh)
+# display(fig)
+
+# fig = Figure()
+# 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)
+# fig
+
+tₑ = 10_000
 
-@info("Precompiling Solver")
-prob = ODEProblem(fₘ, u₀, (0, 1e-4), constants_and_parameters)
-soln = solve(prob, Tsit5())
-soln.retcode != :Unstable || error("Solver was not stable")
 @info("Solving")
 prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
 soln = solve(prob, Tsit5())
 @info("Done")
 
-@save "gray_scott_sphere.jld2" soln
-
-# Visualize the final conditions.
-mesh(s, color=soln(tₑ).U, colormap=:jet)
+function save_dynamics(save_file_name)
+  time = Observable(0.0)
+  u = @lift(soln($time).U)
+  f = Figure()
+  ax_U = CairoMakie.Axis(f[1,1], title = @lift("Concentration of U at Time $($time)"))
 
-begin # BEGIN Gif creation
-frames = 800
-## Initial frame
-fig = Figure(resolution = (1200, 1200))
-p1 = mesh(fig[1,1], s, color=soln(0).U, colormap=:jet, colorrange=extrema(soln(0).U))
-p2 = mesh(fig[2,1], s, color=soln(0).V, colormap=:jet, colorrange=extrema(soln(0).V))
-Colorbar(fig[1,2], colormap=:jet, colorrange=extrema(soln(0).U))
-Colorbar(fig[2,2], colormap=:jet, colorrange=extrema(soln(0).V))
+  msh_U = mesh!(ax_U, s, color=u, colormap=:jet, colorrange=(0, 1.1))
+  Colorbar(f[1,2], msh_U)
 
-## Animation
-record(fig, "gray_scott_sphere.gif", range(0.0, tₑ; length=frames); framerate = 30) do t
-    p1.plot.color = soln(t).U
-    p2.plot.color = soln(t).V
+  timestamps = range(0, tₑ, step=50)
+  record(f, save_file_name, timestamps; framerate = 30) do t
+    time[] = t
+  end
 end
 
-end # END Gif creation
+save_dynamics("gs_f=$(f)_k=$(k).mp4")
\ No newline at end of file

From 3fdd9dea7864bbe602a9818cdce605993a4ed424 Mon Sep 17 00:00:00 2001
From: GeorgeR227 <78235421+GeorgeR227@users.noreply.github.com>
Date: Thu, 19 Dec 2024 17:47:50 -0500
Subject: [PATCH 2/2] Add back Catlab plus updated sources

---
 examples/chemistry/gray_scott.jl | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/examples/chemistry/gray_scott.jl b/examples/chemistry/gray_scott.jl
index 8add5b8b..fbd52830 100644
--- a/examples/chemistry/gray_scott.jl
+++ b/examples/chemistry/gray_scott.jl
@@ -1,3 +1,4 @@
+using Catlab
 using CombinatorialSpaces
 using DiagrammaticEquations
 using Decapodes
@@ -12,9 +13,10 @@ using GeometryBasics: Point2, Point3
 Point2D = Point2{Float64}
 Point3D = Point3{Float64}
 
-# We use the model equations as stated here and use the initial conditions for
-# f, k, rᵤ, rᵥ as listed for experiment 4.
-# https://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/
+# We use the model equations as stated here:
+# https://github.com/JuliaParallel/julia-hpc-tutorial-sc24/blob/main/parts/gpu/gray-scott.ipynb
+# Initial conditions were based off those given here:
+# https://itp.uni-frankfurt.de/~gros/StudentProjects/Projects_2020/projekt_schulz_kaefer/#header
 GrayScott = @decapode begin
   (U, V)::Form0
   (UV2)::Form0
@@ -33,15 +35,16 @@ GrayScott = @decapode begin
 end
 
 n = 100
+h = 1
 
-s = triangulated_grid(n,n,0.8,0.8,Point3D);
+s = triangulated_grid(n,n,h,h,Point3D);
 sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s);
 subdivide_duals!(sd, Circumcenter());
 
 sim = eval(gensim(GrayScott))
 
-left_wall_idxs = findall(x -> x[1] <= 1.0, s[:point])
-right_wall_idxs = findall(x -> x[1] >= n - 1.0, s[:point])
+left_wall_idxs = findall(x -> x[1] <= h, s[:point])
+right_wall_idxs = findall(x -> x[1] >= n - h, s[:point])
 top_wall_idxs = findall(y -> y[2] == 0.0, s[:point])
 bot_wall_idxs = findall(y -> y[2] == n, s[:point])