Skip to content

Commit

Permalink
Update gray_scott.jl
Browse files Browse the repository at this point in the history 
Updates the gray scott model with boundary conditions and constants to produce cool looking behavior. This should probably be turned into a doc page.
  • Loading branch information
@GeorgeR227 @jpfairbanks
GeorgeR227 authored and jpfairbanks committed Jan 10, 2025
1 parent 8cb43b9 commit 19473fe
Showing 1 changed file with 69 additions and 163 deletions.
232 changes: 69 additions & 163 deletions examples/chemistry/gray_scott.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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")

0 comments on commit 19473fe

Please sign in to comment.