add support for gray scott

scripts/run_simulation.jl

@@ -1,13 +1,28 @@
 include("../src/AnimalFurFHN.jl")  # this loads the module code
 using .AnimalFurFHN
 # include optional visualizer only if needed:
-include("../src/visualization.jl")
-using .Visualization
+# include("../src/visualization.jl")
+# using .Visualization
 
-N = 100
-tspan = (0.0, 1000.0)
+# N = 100
+# tspan = (0.0, 1000.0)
 
-sol = AnimalFurFHN.run_simulation(tspan, N)
+# sol = AnimalFurFHN.run_simulation(tspan, N)
 
 
-Visualization.step_through_solution(sol, N)
+# Visualization.step_through_solution(sol, N)
+
+using Plots
+
+function animate_gray_scott(sol, N; var=:U)
+    anim = @animate for t in 1:length(sol.t)
+        u = reshape(sol[t][1:N^2], N, N)
+        v = reshape(sol[t][N^2+1:end], N, N)
+        data = var == :U ? u : v
+        Plots.heatmap(data, c=:magma, title="$var at t=$(Int(sol.t[t]))", clims=(0, 1))
+    end
+    gif(anim, "gray_scott.gif", fps=15)
+end
+
+sol = AnimalFurFHN.run_simulationG((0.0, 10000.0), 256)
+animate_gray_scott(sol, 256, var=:U)

src/constants.jl

@@ -15,6 +15,15 @@ struct FHNParams
         new(N, dx, Du, Dv, ϵ, a, b)
 end
 
-export FHNParams
+struct GSParams
+    N::Int       # grid size
+    dx::Float64  # grid spacing
+    Du::Float64  # diffusion rate U
+    Dv::Float64  # diffusion rate V
+    F::Float64   # feed rate
+    k::Float64   # kill rate
+end
+
+export FHNParams, GSParams
 
 end # module Constants

src/solver.jl

@@ -59,3 +59,57 @@ function run_simulation(tspan::Tuple{Float64,Float64}, N::Int)
 
     return sol
 end
+
+# Laplacian with 5-point stencil
+function laplacianA(A::Matrix{Float64}, N::Int, dx::Float64)
+    dx2 = dx^2
+    lap = zeros(N, N)
+    lap[2:end-1, 2:end-1] .= (
+        A[1:end-2, 2:end-1] .+ A[3:end, 2:end-1]  # up/down
+        .+
+        A[2:end-1, 1:end-2] .+ A[2:end-1, 3:end]  # left/right
+        .-
+        4 .* A[2:end-1, 2:end-1]
+    ) ./ dx2
+    return lap
+end
+
+# Gray-Scott model: reaction + diffusion
+function grayscott!(du, u, p::GSParams, t=0)
+    N = p.N
+    U = reshape(u[1:N^2], N, N)
+    V = reshape(u[N^2+1:end], N, N)
+
+    ΔU = laplacianA(U, N, p.dx)
+    ΔV = laplacianA(V, N, p.dx)
+
+    UVV = U .* V .* V
+
+    dU = p.Du * ΔU .- UVV .+ p.F .* (1 .- U)
+    dV = p.Dv * ΔV .+ UVV .- (p.F .+ p.k) .* V
+
+    du .= vcat(vec(dU), vec(dV))
+end
+
+# Run simulation
+function run_simulationG(tspan::Tuple{Float64,Float64}, N::Int)
+    # Replace this in run_simulation
+    p = GSParams(N, 1.0, 0.16, 0.08, 0.060, 0.062)
+
+    # Initial conditions: U = 1, V = 0 everywhere
+    U = ones(N, N)
+    V = zeros(N, N)
+
+    # Seed a square in the center
+    center = N ÷ 2
+    radius = 10
+    U[center-radius:center+radius, center-radius:center+radius] .= 0.50
+    V[center-radius:center+radius, center-radius:center+radius] .= 0.25
+
+    y0 = vcat(vec(U), vec(V))
+
+    prob = ODEProblem(grayscott!, y0, tspan, p)
+    sol = solve(prob, BS3(), saveat=100.0)  # or Rosenbrock23()
+
+    return sol
+end