From a2b92baffb3797edab16eb63c6d7e89d434e25c2 Mon Sep 17 00:00:00 2001
From: Niki Laptop <2212719@stud.hs-mannheim.de>
Date: Sat, 7 Jun 2025 12:34:38 +0200
Subject: [PATCH] add support for gray scott

---
 scripts/run_simulation.jl | 27 ++++++++++++++----
 src/constants.jl          | 11 +++++++-
 src/solver.jl             | 58 +++++++++++++++++++++++++++++++++++++--
 3 files changed, 87 insertions(+), 9 deletions(-)

diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl
index 3e43bb4..b2a08a8 100644
--- a/scripts/run_simulation.jl
+++ b/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)
\ No newline at end of file
+# 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)
diff --git a/src/constants.jl b/src/constants.jl
index d1734e3..3f6e39a 100644
--- a/src/constants.jl
+++ b/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
diff --git a/src/solver.jl b/src/solver.jl
index 8999319..510cc69 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -17,7 +17,7 @@ using .Constants
     # Returns
     - `du`: calculated derivatives put back into the du array
 """
-function fhn!(du, u, p::FHNParams, t = 0)
+function fhn!(du, u, p::FHNParams, t=0)
     u_mat = reshape(u[1:p.N^2], p.N, p.N) # activation variable
     v_mat = reshape(u[p.N^2+1:end], p.N, p.N) # deactivation variable
 
@@ -44,7 +44,7 @@ end
 """
 function run_simulation(tspan::Tuple{Float64,Float64}, N::Int)
     # Turing-spot parameters
-    p = FHNParams(N = N, dx = 1.0, Du = 1e-5, Dv = 1e-3, ϵ = 0.01, a = 0.1, b = 0.5)
+    p = FHNParams(N=N, dx=1.0, Du=1e-5, Dv=1e-3, ϵ=0.01, a=0.1, b=0.5)
 
     # Initial conditions (random noise)
     Random.seed!(1234)
@@ -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