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 01/23] 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 -- 2.43.0 From 12d2e249acf21ae4d1846caa8fe4ff0652decffb Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Sat, 7 Jun 2025 12:38:31 +0200 Subject: [PATCH 02/23] add gif folder to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index c487889..7a4896f 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,4 @@ docs/site/ # environment. Manifest.toml +gif -- 2.43.0 From b675e17c8824be454f3a21b8f4ab8a8c95e010b0 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Sat, 7 Jun 2025 12:38:44 +0200 Subject: [PATCH 03/23] change save location for gifs --- scripts/run_simulation.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index b2a08a8..274501d 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -21,8 +21,9 @@ function animate_gray_scott(sol, N; var=:U) 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) + gif(anim, "gif/gray_scott.gif", fps=15) end +# set time to more than 1000 sol = AnimalFurFHN.run_simulationG((0.0, 10000.0), 256) animate_gray_scott(sol, 256, var=:U) -- 2.43.0 From 1d3540bfbefc0dddb61255df47317b033d92dd4f Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Sat, 7 Jun 2025 13:33:07 +0200 Subject: [PATCH 04/23] enable stepping through solution with slider in makie --- scripts/run_simulation.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 274501d..bb923e3 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -1,10 +1,10 @@ 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 +N = 256 # tspan = (0.0, 1000.0) # sol = AnimalFurFHN.run_simulation(tspan, N) @@ -26,4 +26,6 @@ end # set time to more than 1000 sol = AnimalFurFHN.run_simulationG((0.0, 10000.0), 256) -animate_gray_scott(sol, 256, var=:U) + +Visualization.step_through_solution(sol, N) +# animate_gray_scott(sol, 256, var=:U) -- 2.43.0 From 3e283fb94e0b04f23f1c309aff8591e70c383488 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Sun, 8 Jun 2025 18:56:06 +0200 Subject: [PATCH 05/23] add function without ode and give GSParams into function --- src/solver.jl | 44 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index 510cc69..259b5d8 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -92,9 +92,9 @@ function grayscott!(du, u, p::GSParams, t=0) end # Run simulation -function run_simulationG(tspan::Tuple{Float64,Float64}, N::Int) +function run_simulationG(tspan::Tuple{Float64,Float64}, N::Int, GSP::GSParams) # Replace this in run_simulation - p = GSParams(N, 1.0, 0.16, 0.08, 0.060, 0.062) + p = GSP # Initial conditions: U = 1, V = 0 everywhere U = ones(N, N) @@ -109,7 +109,45 @@ function run_simulationG(tspan::Tuple{Float64,Float64}, N::Int) y0 = vcat(vec(U), vec(V)) prob = ODEProblem(grayscott!, y0, tspan, p) - sol = solve(prob, BS3(), saveat=100.0) # or Rosenbrock23() + sol = solve(prob, BS3(), saveat=10.0) # or Rosenbrock23() return sol end + +function run_simulationG_no_ode(tspan::Tuple{Float64,Float64}, N::Int, GSP::GSParams) + p = GSP + dx2 = p.dx^2 + dt = 1.0 # fixed timestep + + # Initialize U and V + U = ones(N, N) + V = zeros(N, N) + + # Seed pattern 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 + + # Store history for visualization (optional) + snapshots = [] + + for i in 1:tspan[2] + print("Running $i/$(tspan[2])\r") + + lap_U = laplacianA(U, N, p.dx) + lap_V = laplacianA(V, N, p.dx) + + UVV = U .* V .* V + dU = p.Du * lap_U .- UVV .+ p.F * (1 .- U) + dV = p.Dv * lap_V .+ UVV .- (p.F + p.k) .* V + + U .+= dU * dt + V .+= dV * dt + + if i % 100 == 0 # save every 100 steps + push!(snapshots, (copy(U), copy(V))) + end + end + return snapshots +end -- 2.43.0 From 24d32768c55d0f4f9ca9f1ea1cc96439783750d8 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Sun, 8 Jun 2025 18:56:49 +0200 Subject: [PATCH 06/23] rework run script to enable stepping and simulating while watching --- scripts/run_simulation.jl | 148 +++++++++++++++++++++++++++++++------- 1 file changed, 124 insertions(+), 24 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index bb923e3..6eb584e 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -1,31 +1,131 @@ -include("../src/AnimalFurFHN.jl") # this loads the module code +using GLMakie, Observables using .AnimalFurFHN -# include optional visualizer only if needed: -include("../src/visualization.jl") -using .Visualization - -N = 256 -# tspan = (0.0, 1000.0) - -# sol = AnimalFurFHN.run_simulation(tspan, N) -# Visualization.step_through_solution(sol, N) +# Parameters and initial conditions +N = 128 +dx = 1.0 +Du, Dv = Observable(0.16), Observable(0.08) +F, k = Observable(0.060), Observable(0.062) +dt = 1.0 +params_obs = Observable(AnimalFurFHN.GSParams(N, dx, 0.16, 0.08, 0.08, 0.06)) -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, "gif/gray_scott.gif", fps=15) +function update_params!(params_obs, u, v, feed, kill) + old = params_obs[] + params_obs[] = AnimalFurFHN.GSParams(old.N, old.dx, u, v, feed, kill) end -# set time to more than 1000 -sol = AnimalFurFHN.run_simulationG((0.0, 10000.0), 256) +lift(Du, Dv, F, k) do u, v, F, k + update_params!(params_obs, u, v, F, k) +end +U = ones(N, N) +V = zeros(N, N) +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 -Visualization.step_through_solution(sol, N) -# animate_gray_scott(sol, 256, var=:U) +# Observable holding current U for heatmap +heat_obs = Observable(U) + +function laplacian5(f) + h2 = dx^2 + left = f[2:end-1, 1:end-2] + right = f[2:end-1, 3:end] + down = f[3:end, 2:end-1] + up = f[1:end-2, 2:end-1] + center = f[2:end-1, 2:end-1] + return (left .+ right .+ down .+ up .- 4 .* center) ./ h2 +end + +function step!(U, V) + lap_u = laplacian5(U) + lap_v = laplacian5(V) + + u = U[2:end-1, 2:end-1] + v = V[2:end-1, 2:end-1] + + uvv = u .* v .* v + u_new = u .+ (Du[] .* lap_u .- uvv .+ F[] .* (1 .- u)) .* dt + v_new = v .+ (Dv[] .* lap_v .+ uvv .- (F[] .+ k[]) .* v) .* dt + + # Update with new values + U[2:end-1, 2:end-1] .= u_new + V[2:end-1, 2:end-1] .= v_new + + # Periodic boundary conditions + U[1, :] .= U[end-1, :] + U[end, :] .= U[2, :] + U[:, 1] .= U[:, end-1] + U[:, end] .= U[:, 2] + + V[1, :] .= V[end-1, :] + V[end, :] .= V[2, :] + V[:, 1] .= V[:, end-1] + V[:, end] .= V[:, 2] + + # Update heatmap observable + heat_obs[] = copy(U) +end +function multi_step!(state, n_steps) + for _ in 1:n_steps + step!(state[1], state[2]) + end +end + +# Build GUI +fig = Figure(size=(600, 600)) +ax = Axis(fig[1, 1]) + + + +hm = Makie.heatmap!(ax, heat_obs, colormap=:magma) +ax.aspect = DataAspect() + +# # Controls + +fig[2, 1] = buttongrid = GridLayout(tellwidth=false) +btn_step = Button(buttongrid[1, 1], label="Step") +btn_start = Button(buttongrid[1, 2], label="Start") +btn_stop = Button(buttongrid[1, 3], label="Stop") + + +fig[3, 1] = textboxgrid = GridLayout(tellwidth=false) +box_u = Textbox(textboxgrid[1, 1], validator=Float64, placeholder="word") +# box_v = Textbox(textboxgrid[1, 2], validator=Float64, placeholder="word") +# box_feed = Textbox(textboxgrid[1, 3], validator=Float64, placeholder="word") +# box_kill = Textbox(textboxgrid[1, 4], validator=Float64, placeholder="word") +# Timer and state for animation +running = Observable(false) + +function animation_loop() + while running[] + # step!(U, V) + multi_step!((U, V), 30) + sleep(0.05) # ~20 FPS + end +end +on(box_u.stored_string) do s + try + Du[] = parse(Float64, s) + catch + @warn "Invalid input for Du: $s" + end +end + +on(btn_step.clicks) do _ + multi_step!((U, V), 30) +end + +on(btn_start.clicks) do _ + if !running[] + running[] = true + @async animation_loop() + end +end + +on(btn_stop.clicks) do _ + running[] = false +end + +display(fig) -- 2.43.0 From 6e5cf4b9df06c4b9b96bca387d921856fe7b1ce3 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Sun, 8 Jun 2025 18:57:17 +0200 Subject: [PATCH 07/23] add makie and observables as dependencies --- Manifest.toml | 6 +++--- Project.toml | 2 ++ 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 23b9f57..eceefdd 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.11.5" +julia_version = "1.11.4" manifest_format = "2.0" -project_hash = "0d147eaf78630b05e95b5317275a880443440272" +project_hash = "b5f5f0b50b1e0b7dd05bb93e0b1bb03fea235d53" [[deps.ADTypes]] git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32" @@ -1639,7 +1639,7 @@ version = "3.2.4+0" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.5+0" +version = "0.8.1+4" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] diff --git a/Project.toml b/Project.toml index 19c0af6..0672cef 100644 --- a/Project.toml +++ b/Project.toml @@ -1,4 +1,6 @@ [deps] DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +Observables = "510215fc-4207-5dde-b226-833fc4488ee2" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -- 2.43.0 From ccc364e64469239418464ce87958a8c6232c41c6 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Sun, 8 Jun 2025 18:57:39 +0200 Subject: [PATCH 08/23] add constructor and add more exports --- src/AnimalFurFHN.jl | 2 +- src/constants.jl | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/AnimalFurFHN.jl b/src/AnimalFurFHN.jl index 060abd7..1cde029 100644 --- a/src/AnimalFurFHN.jl +++ b/src/AnimalFurFHN.jl @@ -5,5 +5,5 @@ include("constants.jl") include("laplacian.jl") include("solver.jl") -export run_simulation # Make sure this is here! +export run_simulation, run_simulationG, run_simulationG_no_ode # Make sure this is here! end diff --git a/src/constants.jl b/src/constants.jl index 3f6e39a..4325f91 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -22,6 +22,8 @@ struct GSParams Dv::Float64 # diffusion rate V F::Float64 # feed rate k::Float64 # kill rate + GSParams(; N::Int, dx::Float64, Du::Float64, Dv::Float64, F::Float64, k::Float64) = + new(N, dx, Du, Dv, ϵ, a) end export FHNParams, GSParams -- 2.43.0 From 6bc2c3f06d032b3f616a204bc6d75c75d0184e51 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Sun, 8 Jun 2025 19:08:56 +0200 Subject: [PATCH 09/23] add reset functionality --- scripts/run_simulation.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 6eb584e..943c4cf 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -88,6 +88,7 @@ fig[2, 1] = buttongrid = GridLayout(tellwidth=false) btn_step = Button(buttongrid[1, 1], label="Step") btn_start = Button(buttongrid[1, 2], label="Start") btn_stop = Button(buttongrid[1, 3], label="Stop") +btn_reset = Button(buttongrid[1, 4], label="Reset") fig[3, 1] = textboxgrid = GridLayout(tellwidth=false) @@ -97,6 +98,15 @@ box_u = Textbox(textboxgrid[1, 1], validator=Float64, placeholder="word") # box_kill = Textbox(textboxgrid[1, 4], validator=Float64, placeholder="word") # Timer and state for animation running = Observable(false) +function reset!(U, V, heat_obs) + U .= 1.0 + V .= 0.0 + center = size(U, 1) ÷ 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 + heat_obs[] = copy(U) +end function animation_loop() while running[] @@ -108,6 +118,7 @@ end on(box_u.stored_string) do s try Du[] = parse(Float64, s) + println("Change Du ", Du[]) catch @warn "Invalid input for Du: $s" end @@ -128,4 +139,9 @@ on(btn_stop.clicks) do _ running[] = false end +on(btn_reset.clicks) do _ + running[] = false + reset!(U, V, heat_obs) + +end display(fig) -- 2.43.0 From 34da199ff25a4168ce5760d92c06789e078fd5d0 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Mon, 9 Jun 2025 11:52:59 +0200 Subject: [PATCH 10/23] create value boxes by method, fix heatmap size bug, remove dt from step method --- scripts/run_simulation.jl | 44 ++++++++++++++++++++++++++------------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 943c4cf..e5955a9 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -1,5 +1,6 @@ using GLMakie, Observables -using .AnimalFurFHN +include("../src/constants.jl") +using .Constants # Parameters and initial conditions @@ -8,11 +9,11 @@ dx = 1.0 Du, Dv = Observable(0.16), Observable(0.08) F, k = Observable(0.060), Observable(0.062) dt = 1.0 -params_obs = Observable(AnimalFurFHN.GSParams(N, dx, 0.16, 0.08, 0.08, 0.06)) +params_obs = Observable(GSParams(N, dx, 0.16, 0.08, 0.08, 0.06)) function update_params!(params_obs, u, v, feed, kill) old = params_obs[] - params_obs[] = AnimalFurFHN.GSParams(old.N, old.dx, u, v, feed, kill) + params_obs[] = GSParams(old.N, old.dx, u, v, feed, kill) end lift(Du, Dv, F, k) do u, v, F, k @@ -46,8 +47,8 @@ function step!(U, V) v = V[2:end-1, 2:end-1] uvv = u .* v .* v - u_new = u .+ (Du[] .* lap_u .- uvv .+ F[] .* (1 .- u)) .* dt - v_new = v .+ (Dv[] .* lap_v .+ uvv .- (F[] .+ k[]) .* v) .* dt + u_new = u .+ (Du[] .* lap_u .- uvv .+ F[] .* (1 .- u)) + v_new = v .+ (Dv[] .* lap_v .+ uvv .- (F[] .+ k[]) .* v) # Update with new values U[2:end-1, 2:end-1] .= u_new @@ -65,21 +66,22 @@ function step!(U, V) V[:, end] .= V[:, 2] # Update heatmap observable - heat_obs[] = copy(U) + return U end function multi_step!(state, n_steps) for _ in 1:n_steps - step!(state[1], state[2]) + heat_obs[] = step!(state[1], state[2]) end end # Build GUI fig = Figure(size=(600, 600)) -ax = Axis(fig[1, 1]) +gh = GridLayout(fig[1, 1]) +ax = Axis(gh[1, 1]) -hm = Makie.heatmap!(ax, heat_obs, colormap=:magma) +hm = Makie.heatmap!(ax, heat_obs, colormap=:viridis) ax.aspect = DataAspect() # # Controls @@ -90,12 +92,26 @@ btn_start = Button(buttongrid[1, 2], label="Start") btn_stop = Button(buttongrid[1, 3], label="Stop") btn_reset = Button(buttongrid[1, 4], label="Reset") +gh[1, 2] = textboxgrid = GridLayout(tellwidth=false) +rowsize!(gh, 1, Relative(1.0)) +function param_box!(row, labeltxt, observable) + Label(textboxgrid[row, 1], labeltxt) + box = Textbox(textboxgrid[row, 2], validator=Float64, width=50, placeholder=labeltxt, stored_string="$(observable[])") + on(box.stored_string) do s + try + observable[] = parse(Float64, s) + println("changed $labeltxt to $s") + catch + @warn "Invalid input for $labeltxt: $s" + end + end +end +param_box!(1, "Du", Du) +param_box!(2, "Dv", Dv) +param_box!(3, "Feed", F) +param_box!(4, "kill", k) + -fig[3, 1] = textboxgrid = GridLayout(tellwidth=false) -box_u = Textbox(textboxgrid[1, 1], validator=Float64, placeholder="word") -# box_v = Textbox(textboxgrid[1, 2], validator=Float64, placeholder="word") -# box_feed = Textbox(textboxgrid[1, 3], validator=Float64, placeholder="word") -# box_kill = Textbox(textboxgrid[1, 4], validator=Float64, placeholder="word") # Timer and state for animation running = Observable(false) function reset!(U, V, heat_obs) -- 2.43.0 From 1e9885bd4b85ccfd9479bea906a1ca76529afee4 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Mon, 9 Jun 2025 12:49:03 +0200 Subject: [PATCH 11/23] use params object in step method, preparing for seperating into multiple files --- scripts/run_simulation.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index e5955a9..8bd6c69 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -9,7 +9,7 @@ dx = 1.0 Du, Dv = Observable(0.16), Observable(0.08) F, k = Observable(0.060), Observable(0.062) dt = 1.0 -params_obs = Observable(GSParams(N, dx, 0.16, 0.08, 0.08, 0.06)) +params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[])) function update_params!(params_obs, u, v, feed, kill) old = params_obs[] @@ -39,16 +39,19 @@ function laplacian5(f) return (left .+ right .+ down .+ up .- 4 .* center) ./ h2 end -function step!(U, V) +function step!(U, V, params_obs::Observable) lap_u = laplacian5(U) lap_v = laplacian5(V) - + diff_u = params_obs[].Du + diff_v = params_obs[].Dv + feed_u = params_obs[].F + kill_v = params_obs[].k u = U[2:end-1, 2:end-1] v = V[2:end-1, 2:end-1] uvv = u .* v .* v - u_new = u .+ (Du[] .* lap_u .- uvv .+ F[] .* (1 .- u)) - v_new = v .+ (Dv[] .* lap_v .+ uvv .- (F[] .+ k[]) .* v) + u_new = u .+ (diff_u .* lap_u .- uvv .+ feed_u .* (1 .- u)) + v_new = v .+ (diff_v .* lap_v .+ uvv .- (feed_u .+ kill_v) .* v) # Update with new values U[2:end-1, 2:end-1] .= u_new @@ -70,7 +73,7 @@ function step!(U, V) end function multi_step!(state, n_steps) for _ in 1:n_steps - heat_obs[] = step!(state[1], state[2]) + heat_obs[] = step!(state[1], state[2], params_obs) end end -- 2.43.0 From 951f9281a54d98926aedfb52f94803a7506385b9 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Mon, 9 Jun 2025 12:49:50 +0200 Subject: [PATCH 12/23] remove stop button, and replace with start/pause combined into one button --- scripts/run_simulation.jl | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 8bd6c69..af80c61 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -88,12 +88,11 @@ hm = Makie.heatmap!(ax, heat_obs, colormap=:viridis) ax.aspect = DataAspect() # # Controls - +run_label = Observable("Run") fig[2, 1] = buttongrid = GridLayout(tellwidth=false) btn_step = Button(buttongrid[1, 1], label="Step") -btn_start = Button(buttongrid[1, 2], label="Start") -btn_stop = Button(buttongrid[1, 3], label="Stop") -btn_reset = Button(buttongrid[1, 4], label="Reset") +btn_start = Button(buttongrid[1, 2], width=50, label=run_label) +btn_reset = Button(buttongrid[1, 3], label="Reset") gh[1, 2] = textboxgrid = GridLayout(tellwidth=false) rowsize!(gh, 1, Relative(1.0)) @@ -134,29 +133,24 @@ function animation_loop() sleep(0.05) # ~20 FPS end end -on(box_u.stored_string) do s - try - Du[] = parse(Float64, s) - println("Change Du ", Du[]) - catch - @warn "Invalid input for Du: $s" + +on(running) do r + run_label[] = r ? "Pause" : "Run" end -end + on(btn_step.clicks) do _ multi_step!((U, V), 30) end - on(btn_start.clicks) do _ - if !running[] - running[] = true - @async animation_loop() + running[] = !running[] +end +on(btn_start.clicks) do _ + @async while running[] + animation_loop() end end -on(btn_stop.clicks) do _ - running[] = false -end on(btn_reset.clicks) do _ running[] = false -- 2.43.0 From 50d5766fef6ff041e7e1c4c2e273ad482cff96d3 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Mon, 9 Jun 2025 13:05:25 +0200 Subject: [PATCH 13/23] add documentation and simplify running loop --- scripts/run_simulation.jl | 43 ++++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index af80c61..84e01c2 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -6,7 +6,9 @@ using .Constants # Parameters and initial conditions N = 128 dx = 1.0 +# diffusion rates of substance 'u' and 'v' Du, Dv = Observable(0.16), Observable(0.08) +# feed rate of 'u' and kill rate of 'v' F, k = Observable(0.060), Observable(0.062) dt = 1.0 params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[])) @@ -16,6 +18,7 @@ function update_params!(params_obs, u, v, feed, kill) params_obs[] = GSParams(old.N, old.dx, u, v, feed, kill) end +# Whenever a value gets changed via the textbox, update params object to reflect changes in simulation lift(Du, Dv, F, k) do u, v, F, k update_params!(params_obs, u, v, F, k) end @@ -23,6 +26,7 @@ U = ones(N, N) V = zeros(N, N) center = N ÷ 2 radius = 10 +# set a cube in the center with starting concentrations for 'u' and 'v' U[center-radius:center+radius, center-radius:center+radius] .= 0.50 V[center-radius:center+radius, center-radius:center+radius] .= 0.25 @@ -85,8 +89,11 @@ gh = GridLayout(fig[1, 1]) ax = Axis(gh[1, 1]) hm = Makie.heatmap!(ax, heat_obs, colormap=:viridis) +Makie.deactivate_interaction!(ax, :rectanglezoom) ax.aspect = DataAspect() +spoint = select_point(ax.scene) + # # Controls run_label = Observable("Run") fig[2, 1] = buttongrid = GridLayout(tellwidth=false) @@ -126,32 +133,26 @@ function reset!(U, V, heat_obs) heat_obs[] = copy(U) end -function animation_loop() - while running[] - # step!(U, V) +on(running) do r + run_label[] = r ? "Pause" : "Run" +end + +# Button Listeners +on(btn_step.clicks) do _ + multi_step!((U, V), 30) +end + +on(btn_start.clicks) do _ + running[] = !running[] +end + +on(btn_start.clicks) do _ + @async while running[] multi_step!((U, V), 30) sleep(0.05) # ~20 FPS end end -on(running) do r - run_label[] = r ? "Pause" : "Run" - end - - -on(btn_step.clicks) do _ - multi_step!((U, V), 30) -end -on(btn_start.clicks) do _ - running[] = !running[] -end -on(btn_start.clicks) do _ - @async while running[] - animation_loop() - end -end - - on(btn_reset.clicks) do _ running[] = false reset!(U, V, heat_obs) -- 2.43.0 From 6ec7b9009fd67e91560de8484d093e2d322070ec Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Mon, 9 Jun 2025 13:07:00 +0200 Subject: [PATCH 14/23] remove constructor --- src/constants.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/constants.jl b/src/constants.jl index 4325f91..c06f46c 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -22,8 +22,7 @@ struct GSParams Dv::Float64 # diffusion rate V F::Float64 # feed rate k::Float64 # kill rate - GSParams(; N::Int, dx::Float64, Du::Float64, Dv::Float64, F::Float64, k::Float64) = - new(N, dx, Du, Dv, ϵ, a) + end export FHNParams, GSParams -- 2.43.0 From 5140a2f32c426c6dfe0d336286a4dd37d8b2520a Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Tue, 10 Jun 2025 13:30:19 +0200 Subject: [PATCH 15/23] add feature that you can drop disbalanced concentration at click point, to kickstart reaction there --- scripts/run_simulation.jl | 37 ++++++++++++++++++++++++++++++++----- 1 file changed, 32 insertions(+), 5 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 84e01c2..99c13a4 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -94,14 +94,41 @@ ax.aspect = DataAspect() spoint = select_point(ax.scene) +function coord_to_index(x, y, N) + ix = clamp(round(Int, x), 1, N) + iy = clamp(round(Int, y), 1, N) + return ix, iy +end +r = 5 + +on(spoint) do pt + if pt === nothing + return + end + x, y = pt + i, j = coord_to_index(x, y, N) + + # get corners of square that will get filled with concentration + imin = max(i - r, 1) + imax = min(i + r, N) + jmin = max(j - r, 1) + jmax = min(j + r, N) + + # set disbalanced concentration of U and V + U[imin:imax, jmin:jmax] .= 0.5 + V[imin:imax, jmin:jmax] .= 0.25 + + heat_obs[] = copy(U) +end + # # Controls run_label = Observable("Run") -fig[2, 1] = buttongrid = GridLayout(tellwidth=false) -btn_step = Button(buttongrid[1, 1], label="Step") +fig[2, 1] = buttongrid = GridLayout(ax.scene, tellwidth=false) +btn_step = Button(buttongrid[1, 1], width=50, label="Step") btn_start = Button(buttongrid[1, 2], width=50, label=run_label) -btn_reset = Button(buttongrid[1, 3], label="Reset") +btn_reset = Button(buttongrid[1, 3], width=50, label="Reset") -gh[1, 2] = textboxgrid = GridLayout(tellwidth=false) +gh[1, 2] = textboxgrid = GridLayout(ax.scene, tellwidth=false) rowsize!(gh, 1, Relative(1.0)) function param_box!(row, labeltxt, observable) Label(textboxgrid[row, 1], labeltxt) @@ -148,7 +175,7 @@ end on(btn_start.clicks) do _ @async while running[] - multi_step!((U, V), 30) + multi_step!((U, V), 60) sleep(0.05) # ~20 FPS end end -- 2.43.0 From d50c8437095797e094f4fcd4c401b3c6fc76c6da Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Wed, 11 Jun 2025 11:19:59 +0200 Subject: [PATCH 16/23] add textbox for setting stepsize and some documentation --- scripts/run_simulation.jl | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 99c13a4..effc9a6 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -2,7 +2,7 @@ using GLMakie, Observables include("../src/constants.jl") using .Constants - +# # Gray Scott Model Parameters # Parameters and initial conditions N = 128 dx = 1.0 @@ -13,6 +13,12 @@ F, k = Observable(0.060), Observable(0.062) dt = 1.0 params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[])) +# # GUI Parameters +run_label = Observable{String}("Run") +stepsize = Observable{Int}(30) + + + function update_params!(params_obs, u, v, feed, kill) old = params_obs[] params_obs[] = GSParams(old.N, old.dx, u, v, feed, kill) @@ -122,11 +128,11 @@ on(spoint) do pt end # # Controls -run_label = Observable("Run") fig[2, 1] = buttongrid = GridLayout(ax.scene, tellwidth=false) btn_step = Button(buttongrid[1, 1], width=50, label="Step") btn_start = Button(buttongrid[1, 2], width=50, label=run_label) btn_reset = Button(buttongrid[1, 3], width=50, label="Reset") +step_box = Textbox(buttongrid[1, 4], width=50, validator=Int, placeholder="Stepsize", stored_string="$(stepsize[])") gh[1, 2] = textboxgrid = GridLayout(ax.scene, tellwidth=false) rowsize!(gh, 1, Relative(1.0)) @@ -164,9 +170,18 @@ on(running) do r run_label[] = r ? "Pause" : "Run" end +on(step_box.stored_string) do s + try + stepsize[] = parse(Int, s) + println("Changed stepsize to $s") + catch + @warn "Invalid input for $labeltxt: $s" + end + +end # Button Listeners on(btn_step.clicks) do _ - multi_step!((U, V), 30) + multi_step!((U, V), stepsize[]) end on(btn_start.clicks) do _ @@ -175,7 +190,7 @@ end on(btn_start.clicks) do _ @async while running[] - multi_step!((U, V), 60) + multi_step!((U, V), stepsize[]) sleep(0.05) # ~20 FPS end end -- 2.43.0 From 2819c48076180571004b52f18183f3a7d29758bc Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Wed, 11 Jun 2025 11:44:36 +0200 Subject: [PATCH 17/23] replace step textbox with slider --- scripts/run_simulation.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index effc9a6..4949395 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -132,7 +132,17 @@ fig[2, 1] = buttongrid = GridLayout(ax.scene, tellwidth=false) btn_step = Button(buttongrid[1, 1], width=50, label="Step") btn_start = Button(buttongrid[1, 2], width=50, label=run_label) btn_reset = Button(buttongrid[1, 3], width=50, label="Reset") -step_box = Textbox(buttongrid[1, 4], width=50, validator=Int, placeholder="Stepsize", stored_string="$(stepsize[])") +slidergrid = SliderGrid(fig[3, 1], (label="Speed", range=1:1:100, format="{}x", width=350, startvalue=stepsize[])) + +speed_slider = slidergrid.sliders[1].value +on(speed_slider) do s + try + stepsize[] = s + println("Changed stepsize to $s") + catch + @warn "Invalid input for $s" + end +end gh[1, 2] = textboxgrid = GridLayout(ax.scene, tellwidth=false) rowsize!(gh, 1, Relative(1.0)) @@ -170,15 +180,7 @@ on(running) do r run_label[] = r ? "Pause" : "Run" end -on(step_box.stored_string) do s - try - stepsize[] = parse(Int, s) - println("Changed stepsize to $s") - catch - @warn "Invalid input for $labeltxt: $s" - end -end # Button Listeners on(btn_step.clicks) do _ multi_step!((U, V), stepsize[]) -- 2.43.0 From 32f3885e86959efb360b2a6ef363710b24e87f57 Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Wed, 11 Jun 2025 14:48:24 +0200 Subject: [PATCH 18/23] split gray scott simulation into multiple files --- scripts/main.jl | 28 ++++++++++ src/gray_scott_solver.jl | 55 ++++++++++++++++++++ src/visualization.jl | 109 +++++++++++++++++++++++++++++++++++++-- 3 files changed, 189 insertions(+), 3 deletions(-) create mode 100644 scripts/main.jl create mode 100644 src/gray_scott_solver.jl diff --git a/scripts/main.jl b/scripts/main.jl new file mode 100644 index 0000000..e5566c6 --- /dev/null +++ b/scripts/main.jl @@ -0,0 +1,28 @@ +include("../src/gray_scott_solver.jl") +include("../src/visualization.jl") +include("../src/constants.jl") + + +using Observables +using GLMakie +using .GrayScottSolver +using .Visualization +using .Constants + +const N = 128 +const dx = 1.0 +Du, Dv = Observable(0.16), Observable(0.08) +F, k = Observable(0.060), Observable(0.062) + +params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[])) +lift(Du, Dv, F, k) do u, v, f, ki + params_obs[] = GSParams(N, dx, u, v, f, ki) +end + + +U = ones(N, N) +V = zeros(N, N) +heat_obs = Observable(U) + +fig = build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs) +display(fig) diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl new file mode 100644 index 0000000..ca1840b --- /dev/null +++ b/src/gray_scott_solver.jl @@ -0,0 +1,55 @@ +module GrayScottSolver +using Observables +include("constants.jl") +using .Constants +export laplacian5, step!, multi_step! + + +function laplacian5(f, dx) + h2 = dx^2 + left = f[2:end-1, 1:end-2] + right = f[2:end-1, 3:end] + down = f[3:end, 2:end-1] + up = f[1:end-2, 2:end-1] + center = f[2:end-1, 2:end-1] + return (left .+ right .+ down .+ up .- 4 .* center) ./ h2 +end + +function step!(U, V, params_obs::Observable; dx=1) + lap_u = laplacian5(U, dx) + lap_v = laplacian5(V, dx) + diff_u = params_obs[].Du + diff_v = params_obs[].Dv + feed_u = params_obs[].F + kill_v = params_obs[].k + u = U[2:end-1, 2:end-1] + v = V[2:end-1, 2:end-1] + + uvv = u .* v .* v + u_new = u .+ (diff_u .* lap_u .- uvv .+ feed_u .* (1 .- u)) + v_new = v .+ (diff_v .* lap_v .+ uvv .- (feed_u .+ kill_v) .* v) + + # Update with new values + U[2:end-1, 2:end-1] .= u_new + V[2:end-1, 2:end-1] .= v_new + + # Periodic boundary conditions + U[1, :] .= U[end-1, :] + U[end, :] .= U[2, :] + U[:, 1] .= U[:, end-1] + U[:, end] .= U[:, 2] + + V[1, :] .= V[end-1, :] + V[end, :] .= V[2, :] + V[:, 1] .= V[:, end-1] + V[:, end] .= V[:, 2] + + # Update heatmap observable + return U +end +function multi_step!(state, n_steps, heat_obs::Observable, params_obs::Observable; dx=1) + for _ in 1:n_steps + heat_obs[] = step!(state[1], state[2], params_obs; dx=1) + end +end +end \ No newline at end of file diff --git a/src/visualization.jl b/src/visualization.jl index 7f1514a..ec4e778 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -1,6 +1,7 @@ module Visualization - -using GLMakie +include("gray_scott_solver.jl") +using GLMakie, Observables, Makie +using .GrayScottSolver """ step_through_solution(sol::SolutionType, N::Int) @@ -18,18 +19,120 @@ function step_through_solution(sol, N::Int) ax = Axis(fig[1, 1]) slider = Slider(fig[2, 1], range=1:length(sol), startvalue=1) + textbox = Textbox(fig[1, 2], placeholder="Feed Rate", validator=Float64) + F = Observable(0.060) # Initialize heatmap with first time step u0 = reshape(sol[1][1:N^2], N, N) heat_obs = Observable(u0) hmap = heatmap!(ax, heat_obs, colormap=:magma) + on(textbox.stored_string) do s + F[] = parse(Float64, s) + end # Update heatmap on slider movement on(slider.value) do i u = reshape(sol[i][1:N^2], N, N) heat_obs[] = u end + display(fig) end - +function coord_to_index(x, y, N) + ix = clamp(round(Int, x), 1, N) + iy = clamp(round(Int, y), 1, N) + return ix, iy +end +function reset!(U, V, heat_obs) + U .= 1.0 + V .= 0.0 + center = size(U, 1) ÷ 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 + heat_obs[] = copy(U) +end +function param_box!(grid, row, labeltxt, observable::Observable) + Label(grid[row, 1], labeltxt) + box = Textbox(grid[row, 2], validator=Float64, width=50, placeholder=labeltxt, stored_string="$(observable[])") + on(box.stored_string) do s + try + observable[] = parse(Float64, s) + println("changed $labeltxt to $s") + catch + @warn "Invalid input for $labeltxt: $s" + end + end +end + +function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs) + reset!(U, V, heat_obs) + fig = Figure(size=(800, 800)) + gh = GridLayout(fig[1, 1]) + ax = Axis(gh[1, 1]) + hm = Makie.heatmap!(ax, heat_obs, colormap=:viridis) + deactivate_interaction!(ax, :rectanglezoom) + ax.aspect = DataAspect() + + run_label = Observable{String}("Run") + stepsize = Observable(30) + # # Controls + fig[2, 1] = buttongrid = GridLayout(ax.scene, tellwidth=false) + btn_step = Button(buttongrid[1, 1], width=50, label="Step") + btn_start = Button(buttongrid[1, 2], width=50, label=run_label) + btn_reset = Button(buttongrid[1, 3], width=50, label="Reset") + slidergrid = SliderGrid(fig[3, 1], (label="Speed", range=1:1:100, format="{}x", width=350, startvalue=stepsize[])) + + speed_slider = slidergrid.sliders[1].value + + + gh[1, 2] = textboxgrid = GridLayout(ax.scene, tellwidth=false) + rowsize!(gh, 1, Relative(1.0)) + + + param_box!(textboxgrid, 1, "Du", Du) + param_box!(textboxgrid, 2, "Dv", Dv) + param_box!(textboxgrid, 3, "Feed", F) + param_box!(textboxgrid, 4, "kill", k) + + # Timer and state for animation + running = Observable(false) + + on(running) do r + run_label[] = r ? "Pause" : "Run" + end + on(speed_slider) do s + try + stepsize[] = s + println("Changed stepsize to $s") + catch + @warn "Invalid input for $s" + end + end + # Button Listeners + on(btn_step.clicks) do _ + multi_step!((U, V), stepsize[], heat_obs, params_obs) + end + + on(btn_start.clicks) do _ + running[] = !running[] + end + + on(btn_start.clicks) do _ + @async while running[] + multi_step!((U, V), stepsize[], heat_obs, params_obs) + sleep(0.0015) # ~20 FPS + end + end + + on(btn_reset.clicks) do _ + running[] = false + reset!(U, V, heat_obs) + + end + return fig +end + + +export step_through_solution, build_ui end -- 2.43.0 From 249715cbf08a763212a6ff3463ac938269f6e6fe Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Wed, 11 Jun 2025 14:48:42 +0200 Subject: [PATCH 19/23] delete unnecessary module --- src/AnimalFurFHN.jl | 9 --------- 1 file changed, 9 deletions(-) delete mode 100644 src/AnimalFurFHN.jl diff --git a/src/AnimalFurFHN.jl b/src/AnimalFurFHN.jl deleted file mode 100644 index 1cde029..0000000 --- a/src/AnimalFurFHN.jl +++ /dev/null @@ -1,9 +0,0 @@ -# src/AnimalFurFHN.jl -module AnimalFurFHN - -include("constants.jl") -include("laplacian.jl") -include("solver.jl") - -export run_simulation, run_simulationG, run_simulationG_no_ode # Make sure this is here! -end -- 2.43.0 From 33a092034e3beaf934d3d3316bbe04ab9fb7a31f Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Wed, 11 Jun 2025 14:49:00 +0200 Subject: [PATCH 20/23] add type annotation --- scripts/run_simulation.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 4949395..226ccf6 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -1,6 +1,8 @@ using GLMakie, Observables include("../src/constants.jl") using .Constants +include("../src/visualization.jl") +using .Visualization # # Gray Scott Model Parameters # Parameters and initial conditions @@ -19,7 +21,7 @@ stepsize = Observable{Int}(30) -function update_params!(params_obs, u, v, feed, kill) +function update_params!(params_obs::Observable, u, v, feed, kill) old = params_obs[] params_obs[] = GSParams(old.N, old.dx, u, v, feed, kill) end @@ -39,6 +41,7 @@ V[center-radius:center+radius, center-radius:center+radius] .= 0.25 # Observable holding current U for heatmap heat_obs = Observable(U) + function laplacian5(f) h2 = dx^2 left = f[2:end-1, 1:end-2] -- 2.43.0 From 2fc80227ca219080b459010adf5f25a0dd036efb Mon Sep 17 00:00:00 2001 From: Niki Laptop <2212719@stud.hs-mannheim.de> Date: Wed, 11 Jun 2025 14:49:14 +0200 Subject: [PATCH 21/23] include file so it works --- src/solver.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/solver.jl b/src/solver.jl index 259b5d8..3910c8f 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -1,5 +1,6 @@ using DifferentialEquations using Random +include("constants.jl") using .Constants """ -- 2.43.0 From f39c85a0adc4eb3537217db7471ccfa909c45226 Mon Sep 17 00:00:00 2001 From: Nikola <2212719@stud.hs-mannheim.de> Date: Thu, 12 Jun 2025 18:22:55 +0200 Subject: [PATCH 22/23] add functionality to add concentration at clickend spots dynamically --- src/visualization.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/visualization.jl b/src/visualization.jl index ec4e778..7d1a066 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -76,6 +76,8 @@ function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs) run_label = Observable{String}("Run") stepsize = Observable(30) + spoint = select_point(ax.scene) + # # Controls fig[2, 1] = buttongrid = GridLayout(ax.scene, tellwidth=false) btn_step = Button(buttongrid[1, 1], width=50, label="Step") @@ -130,6 +132,27 @@ function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs) reset!(U, V, heat_obs) end + on(spoint) do pt + r = 5 + if pt === nothing + return + end + x, y = pt + println("params_obs[].N, ", params_obs[].N) + i, j = coord_to_index(x, y, params_obs[].N) + + # get corners of square that will get filled with concentration + imin = max(i - r, 1) + imax = min(i + r, params_obs[].N) + jmin = max(j - r, 1) + jmax = min(j + r, params_obs[].N) + + # set disbalanced concentration of U and V + U[imin:imax, jmin:jmax] .= 0.5 + V[imin:imax, jmin:jmax] .= 0.25 + + heat_obs[] = copy(U) + end return fig end -- 2.43.0 From c98b5f0fca81354a08bb6d2ee5891fe8d317689b Mon Sep 17 00:00:00 2001 From: 2211567 Date: Fri, 13 Jun 2025 14:08:05 +0200 Subject: [PATCH 23/23] created utils file; rm two files --- scripts/main.jl | 2 +- scripts/run_simulation.jl | 208 ----------------------------------- src/gray_scott_solver.jl | 17 +-- src/solver.jl | 154 -------------------------- src/{ => utils}/constants.jl | 0 src/{ => utils}/laplacian.jl | 16 +++ src/visualization.jl | 3 + 7 files changed, 24 insertions(+), 376 deletions(-) delete mode 100644 scripts/run_simulation.jl delete mode 100644 src/solver.jl rename src/{ => utils}/constants.jl (100%) rename src/{ => utils}/laplacian.jl (71%) diff --git a/scripts/main.jl b/scripts/main.jl index e5566c6..e38e373 100644 --- a/scripts/main.jl +++ b/scripts/main.jl @@ -1,6 +1,6 @@ include("../src/gray_scott_solver.jl") include("../src/visualization.jl") -include("../src/constants.jl") +include("../src/utils/constants.jl") using Observables diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl deleted file mode 100644 index 226ccf6..0000000 --- a/scripts/run_simulation.jl +++ /dev/null @@ -1,208 +0,0 @@ -using GLMakie, Observables -include("../src/constants.jl") -using .Constants -include("../src/visualization.jl") -using .Visualization - -# # Gray Scott Model Parameters -# Parameters and initial conditions -N = 128 -dx = 1.0 -# diffusion rates of substance 'u' and 'v' -Du, Dv = Observable(0.16), Observable(0.08) -# feed rate of 'u' and kill rate of 'v' -F, k = Observable(0.060), Observable(0.062) -dt = 1.0 -params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[])) - -# # GUI Parameters -run_label = Observable{String}("Run") -stepsize = Observable{Int}(30) - - - -function update_params!(params_obs::Observable, u, v, feed, kill) - old = params_obs[] - params_obs[] = GSParams(old.N, old.dx, u, v, feed, kill) -end - -# Whenever a value gets changed via the textbox, update params object to reflect changes in simulation -lift(Du, Dv, F, k) do u, v, F, k - update_params!(params_obs, u, v, F, k) -end -U = ones(N, N) -V = zeros(N, N) -center = N ÷ 2 -radius = 10 -# set a cube in the center with starting concentrations for 'u' and 'v' -U[center-radius:center+radius, center-radius:center+radius] .= 0.50 -V[center-radius:center+radius, center-radius:center+radius] .= 0.25 - -# Observable holding current U for heatmap -heat_obs = Observable(U) - - -function laplacian5(f) - h2 = dx^2 - left = f[2:end-1, 1:end-2] - right = f[2:end-1, 3:end] - down = f[3:end, 2:end-1] - up = f[1:end-2, 2:end-1] - center = f[2:end-1, 2:end-1] - return (left .+ right .+ down .+ up .- 4 .* center) ./ h2 -end - -function step!(U, V, params_obs::Observable) - lap_u = laplacian5(U) - lap_v = laplacian5(V) - diff_u = params_obs[].Du - diff_v = params_obs[].Dv - feed_u = params_obs[].F - kill_v = params_obs[].k - u = U[2:end-1, 2:end-1] - v = V[2:end-1, 2:end-1] - - uvv = u .* v .* v - u_new = u .+ (diff_u .* lap_u .- uvv .+ feed_u .* (1 .- u)) - v_new = v .+ (diff_v .* lap_v .+ uvv .- (feed_u .+ kill_v) .* v) - - # Update with new values - U[2:end-1, 2:end-1] .= u_new - V[2:end-1, 2:end-1] .= v_new - - # Periodic boundary conditions - U[1, :] .= U[end-1, :] - U[end, :] .= U[2, :] - U[:, 1] .= U[:, end-1] - U[:, end] .= U[:, 2] - - V[1, :] .= V[end-1, :] - V[end, :] .= V[2, :] - V[:, 1] .= V[:, end-1] - V[:, end] .= V[:, 2] - - # Update heatmap observable - return U -end -function multi_step!(state, n_steps) - for _ in 1:n_steps - heat_obs[] = step!(state[1], state[2], params_obs) - end -end - -# Build GUI -fig = Figure(size=(600, 600)) - - -gh = GridLayout(fig[1, 1]) -ax = Axis(gh[1, 1]) - -hm = Makie.heatmap!(ax, heat_obs, colormap=:viridis) -Makie.deactivate_interaction!(ax, :rectanglezoom) -ax.aspect = DataAspect() - -spoint = select_point(ax.scene) - -function coord_to_index(x, y, N) - ix = clamp(round(Int, x), 1, N) - iy = clamp(round(Int, y), 1, N) - return ix, iy -end -r = 5 - -on(spoint) do pt - if pt === nothing - return - end - x, y = pt - i, j = coord_to_index(x, y, N) - - # get corners of square that will get filled with concentration - imin = max(i - r, 1) - imax = min(i + r, N) - jmin = max(j - r, 1) - jmax = min(j + r, N) - - # set disbalanced concentration of U and V - U[imin:imax, jmin:jmax] .= 0.5 - V[imin:imax, jmin:jmax] .= 0.25 - - heat_obs[] = copy(U) -end - -# # Controls -fig[2, 1] = buttongrid = GridLayout(ax.scene, tellwidth=false) -btn_step = Button(buttongrid[1, 1], width=50, label="Step") -btn_start = Button(buttongrid[1, 2], width=50, label=run_label) -btn_reset = Button(buttongrid[1, 3], width=50, label="Reset") -slidergrid = SliderGrid(fig[3, 1], (label="Speed", range=1:1:100, format="{}x", width=350, startvalue=stepsize[])) - -speed_slider = slidergrid.sliders[1].value -on(speed_slider) do s - try - stepsize[] = s - println("Changed stepsize to $s") - catch - @warn "Invalid input for $s" - end -end - -gh[1, 2] = textboxgrid = GridLayout(ax.scene, tellwidth=false) -rowsize!(gh, 1, Relative(1.0)) -function param_box!(row, labeltxt, observable) - Label(textboxgrid[row, 1], labeltxt) - box = Textbox(textboxgrid[row, 2], validator=Float64, width=50, placeholder=labeltxt, stored_string="$(observable[])") - on(box.stored_string) do s - try - observable[] = parse(Float64, s) - println("changed $labeltxt to $s") - catch - @warn "Invalid input for $labeltxt: $s" - end - end -end -param_box!(1, "Du", Du) -param_box!(2, "Dv", Dv) -param_box!(3, "Feed", F) -param_box!(4, "kill", k) - - -# Timer and state for animation -running = Observable(false) -function reset!(U, V, heat_obs) - U .= 1.0 - V .= 0.0 - center = size(U, 1) ÷ 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 - heat_obs[] = copy(U) -end - -on(running) do r - run_label[] = r ? "Pause" : "Run" -end - - -# Button Listeners -on(btn_step.clicks) do _ - multi_step!((U, V), stepsize[]) -end - -on(btn_start.clicks) do _ - running[] = !running[] -end - -on(btn_start.clicks) do _ - @async while running[] - multi_step!((U, V), stepsize[]) - sleep(0.05) # ~20 FPS - end -end - -on(btn_reset.clicks) do _ - running[] = false - reset!(U, V, heat_obs) - -end -display(fig) diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl index ca1840b..3668d26 100644 --- a/src/gray_scott_solver.jl +++ b/src/gray_scott_solver.jl @@ -1,19 +1,10 @@ module GrayScottSolver using Observables -include("constants.jl") +include("utils/constants.jl") +include("utils/laplacian.jl") using .Constants -export laplacian5, step!, multi_step! - - -function laplacian5(f, dx) - h2 = dx^2 - left = f[2:end-1, 1:end-2] - right = f[2:end-1, 3:end] - down = f[3:end, 2:end-1] - up = f[1:end-2, 2:end-1] - center = f[2:end-1, 2:end-1] - return (left .+ right .+ down .+ up .- 4 .* center) ./ h2 -end +using .Laplacian +export step!, multi_step! function step!(U, V, params_obs::Observable; dx=1) lap_u = laplacian5(U, dx) diff --git a/src/solver.jl b/src/solver.jl deleted file mode 100644 index 3910c8f..0000000 --- a/src/solver.jl +++ /dev/null @@ -1,154 +0,0 @@ -using DifferentialEquations -using Random -include("constants.jl") -using .Constants - -""" - fhn(du, u, p:FHNParams, t:) - - Implements the spatial dynamics of FitzHugh-Nagumo (fhn). Designed to be - within a larger numerical solver of partial differential equations. - - # Arguments - - `du`: output argument which stores the calculated derivatives - - `u`: input vector containing the current state of the system at time t - - `p`: holds all the fixed parameters of the FHN model - - `t`: current time - - # Returns - - `du`: calculated derivatives put back into the du array -""" -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 - - Δu = reshape(laplacian(u_mat, p.N, p.dx), p.N, p.N) - Δv = reshape(laplacian(v_mat, p.N, p.dx), p.N, p.N) - - fu = p.Du * Δu .+ u_mat .- u_mat .^ 3 ./ 3 .- v_mat - fv = p.Dv * Δv .+ p.ϵ * (u_mat .+ p.a .- p.b .* v_mat) - - du .= vcat(vec(fu), vec(fv)) -end - -""" - run_simulation(tspan::Tuple{Float64,Float64}, N::Int) - - solving the ODE and modelling it after FHN - - # Arguments - - `tspan`: tuple of two Float64's representing start and end times for simulation - - `N`: size of the N×N grid - - # Returns - - `sol`: solved differential equation (ODE) -""" -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) - - # Initial conditions (random noise) - Random.seed!(1234) - # Create two vectors with length N*N with numbers between 0.1 and 0.11 - u0 = vec(0.1 .+ 0.01 .* rand(N, N)) - v0 = vec(0.1 .+ 0.01 .* rand(N, N)) - - y0 = vcat(u0, v0) - - prob = ODEProblem(fhn!, y0, tspan, p) - sol = solve(prob, BS3(), saveat=1.0) # You can try `Rosenbrock23()` too - - 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, GSP::GSParams) - # Replace this in run_simulation - p = GSP - - # 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=10.0) # or Rosenbrock23() - - return sol -end - -function run_simulationG_no_ode(tspan::Tuple{Float64,Float64}, N::Int, GSP::GSParams) - p = GSP - dx2 = p.dx^2 - dt = 1.0 # fixed timestep - - # Initialize U and V - U = ones(N, N) - V = zeros(N, N) - - # Seed pattern 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 - - # Store history for visualization (optional) - snapshots = [] - - for i in 1:tspan[2] - print("Running $i/$(tspan[2])\r") - - lap_U = laplacianA(U, N, p.dx) - lap_V = laplacianA(V, N, p.dx) - - UVV = U .* V .* V - dU = p.Du * lap_U .- UVV .+ p.F * (1 .- U) - dV = p.Dv * lap_V .+ UVV .- (p.F + p.k) .* V - - U .+= dU * dt - V .+= dV * dt - - if i % 100 == 0 # save every 100 steps - push!(snapshots, (copy(U), copy(V))) - end - end - return snapshots -end diff --git a/src/constants.jl b/src/utils/constants.jl similarity index 100% rename from src/constants.jl rename to src/utils/constants.jl diff --git a/src/laplacian.jl b/src/utils/laplacian.jl similarity index 71% rename from src/laplacian.jl rename to src/utils/laplacian.jl index f9f4965..e04d5a2 100644 --- a/src/laplacian.jl +++ b/src/utils/laplacian.jl @@ -1,3 +1,5 @@ +module Laplacian + """ laplacian(U::Matrix{Float64}, N::Int, h::Float64) @@ -18,3 +20,17 @@ function laplacian(U::Matrix{Float64}, N::Int, h::Float64) circshift(U, (0, -1)) .+ circshift(U, (0, 1)) .- 4 .* U return vec(padded) ./ h^2 end + +function laplacian5(f, dx) + h2 = dx^2 + left = f[2:end-1, 1:end-2] + right = f[2:end-1, 3:end] + down = f[3:end, 2:end-1] + up = f[1:end-2, 2:end-1] + center = f[2:end-1, 2:end-1] + return (left .+ right .+ down .+ up .- 4 .* center) ./ h2 +end + +export laplacian, laplacian5 + +end # module Laplacian \ No newline at end of file diff --git a/src/visualization.jl b/src/visualization.jl index 7d1a066..d2f8343 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -38,11 +38,13 @@ function step_through_solution(sol, N::Int) display(fig) end + function coord_to_index(x, y, N) ix = clamp(round(Int, x), 1, N) iy = clamp(round(Int, y), 1, N) return ix, iy end + function reset!(U, V, heat_obs) U .= 1.0 V .= 0.0 @@ -52,6 +54,7 @@ function reset!(U, V, heat_obs) V[center-radius:center+radius, center-radius:center+radius] .= 0.25 heat_obs[] = copy(U) end + function param_box!(grid, row, labeltxt, observable::Observable) Label(grid[row, 1], labeltxt) box = Textbox(grid[row, 2], validator=Float64, width=50, placeholder=labeltxt, stored_string="$(observable[])") -- 2.43.0