From 7b761e7ff994e31627b6c1131a700a4ad8c0eb03 Mon Sep 17 00:00:00 2001 From: 2211567 Date: Sun, 8 Jun 2025 17:48:09 +0200 Subject: [PATCH 01/11] FHN model WIP --- Manifest.toml | 2 +- Project.toml | 1 + scripts/run_simulation.jl | 4 +- src/laplacian.jl | 11 +++--- src/solver.jl | 79 ++++++++++++++++++++++++++++++++++++--- src/visualization.jl | 2 +- 6 files changed, 84 insertions(+), 15 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 23b9f57..6aaa7a8 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.11.5" manifest_format = "2.0" -project_hash = "0d147eaf78630b05e95b5317275a880443440272" +project_hash = "0465e7f271b455d66e099c4ff05fca7f6e5b2b64" [[deps.ADTypes]] git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32" diff --git a/Project.toml b/Project.toml index 19c0af6..de40aee 100644 --- a/Project.toml +++ b/Project.toml @@ -1,4 +1,5 @@ [deps] DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 3e43bb4..6146ea8 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -4,8 +4,8 @@ using .AnimalFurFHN include("../src/visualization.jl") using .Visualization -N = 100 -tspan = (0.0, 1000.0) +N = 256 +tspan = (0.0, 5000.0) sol = AnimalFurFHN.run_simulation(tspan, N) diff --git a/src/laplacian.jl b/src/laplacian.jl index f9f4965..6a04d64 100644 --- a/src/laplacian.jl +++ b/src/laplacian.jl @@ -13,8 +13,9 @@ - `Vector{Float64}`: A flattened (vectorized) representation of the approximated Laplacian values for each element in `U`. The boundary conditions are handled circularly. """ function laplacian(U::Matrix{Float64}, N::Int, h::Float64) - # shifts matrices and sums them up - padded = circshift(U, (-1, 0)) .+ circshift(U, (1, 0)) .+ - circshift(U, (0, -1)) .+ circshift(U, (0, 1)) .- 4 .* U - return vec(padded) ./ h^2 -end + Δ = zeros(N, N) + for i in 2:N-1, j in 2:N-1 + Δ[i, j] = (U[i+1, j] + U[i-1, j] + U[i, j+1] + U[i, j-1] - 4 * U[i, j]) / h^2 + end + return vec(Δ) +end \ No newline at end of file diff --git a/src/solver.jl b/src/solver.jl index 8999319..3ddd9f0 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 @@ -30,6 +30,70 @@ function fhn!(du, u, p::FHNParams, t = 0) du .= vcat(vec(fu), vec(fv)) end +# helper functions for filling cells in specific places of the matrix +function blocks_ic(N) + u = fill(1.0, N, N) + v = fill(0.0, N, N) + p = div(N, 2) + + function safe_block!(u, row_center, col_center) + row_start = max(row_center - 8, 1) + row_end = min(row_center + 7, N) + col_start = max(col_center - 8, 1) + col_end = min(col_center + 7, N) + u[row_start:row_end, col_start:col_end] .= -0.534522 + end + + safe_block!(u, p, p) + + return vec(u), vec(v) +end + +function column_ic(N) + u = fill(0.534522, N, N) + v = fill(0.381802, N, N) + + col_center = div(N, 2) + col_width = 8 # You can adjust this + + col_start = max(col_center - div(col_width, 2), 1) + col_end = min(col_center + div(col_width, 2) - 1, N) + + u[col_start:col_end, :] .= -0.534522 + + return vec(u), vec(v) +end + +function center_band_ic(N) + u = fill(0.0, N, N) + v = fill(0.0, N, N) + + band_width = div(N, 8) + row_start = div(N, 2) - div(band_width, 2) + row_end = div(N, 2) + div(band_width, 2) + + u[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N) + v[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N) + + return vec(u), vec(v) +end + +function circle_ic(N) + u = fill(0.534522, N, N) + v = fill(0.381802, N, N) + cx, cy = div(N, 2), div(N, 2) # center of matrix + radius = 0.250 * N # circle radius = 3/4 of N divided by 2 + + for i in 1:N, j in 1:N + if (i - cx)^2 + (j - cy)^2 ≤ radius^2 + u[i, j] = -0.534522 + end + end + + return vec(u), vec(v) +end + + """ run_simulation(tspan::Tuple{Float64,Float64}, N::Int) @@ -44,18 +108,21 @@ 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=0.00016, Dv=0.001, ϵ=0.1, a=0.5, b=0.9) # 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)) + #u0 = vec(0.4 .+ 0.01 .* rand(N, N)) + #v0 = vec(0.4 .+ 0.01 .* rand(N, N)) - y0 = vcat(u0, v0) + # Or use this + u0, v0 = center_band_ic(N) + + y0 = vcat(vcat(u0, v0)) prob = ODEProblem(fhn!, y0, tspan, p) - sol = solve(prob, BS3(), saveat=1.0) # You can try `Rosenbrock23()` too + sol = solve(prob, BS3(), saveat=3.0) # You can try `Rosenbrock23()` too return sol end diff --git a/src/visualization.jl b/src/visualization.jl index 7f1514a..00a0b71 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -21,7 +21,7 @@ function step_through_solution(sol, N::Int) # 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) + hmap = heatmap!(ax, heat_obs, colormap=:berlin) # Update heatmap on slider movement on(slider.value) do i From 2f498367e8dfa7797a8a0192794fd38041cb0559 Mon Sep 17 00:00:00 2001 From: 2211567 Date: Tue, 10 Jun 2025 12:07:13 +0200 Subject: [PATCH 02/11] adding a squiggle on the matrix now possible --- src/solver.jl | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index 3ddd9f0..16287aa 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -93,6 +93,33 @@ function circle_ic(N) return vec(u), vec(v) end +function squiggle_ic(N, Lx=400.0, Ly=400.0) + uplus = 0.534522 + vplus = 0.381802 + uminus = -uplus + + # Create coordinate grids + x = LinRange(0, Lx, N) + y = LinRange(0, Ly, N) + X = repeat(x', N, 1) # Transposed to align with meshgrid X + Y = repeat(y, 1, N) # Broadcasted to align with meshgrid Y + + # Squiggle pattern + cos_term = 0.05 * Lx .* sin.(10 * 2π .* Y ./ Ly .+ π * 0.3) + exp_term = exp.(-((Y .- Ly / 2) ./ (0.1 * Ly)) .^ 2) + width = 0.05 * Lx + Z = exp.(-((X .- Lx / 2 .+ cos_term .* exp_term) ./ width) .^ 2) + + + u = fill(uplus, N, N) + v = fill(vplus, N, N) + + # Apply squiggle + u[Z .> 0.8] .= uminus + + return vec(u), vec(v) +end + """ run_simulation(tspan::Tuple{Float64,Float64}, N::Int) @@ -117,12 +144,12 @@ function run_simulation(tspan::Tuple{Float64,Float64}, N::Int) #v0 = vec(0.4 .+ 0.01 .* rand(N, N)) # Or use this - u0, v0 = center_band_ic(N) + u0, v0 = squiggle_ic(N) y0 = vcat(vcat(u0, v0)) prob = ODEProblem(fhn!, y0, tspan, p) - sol = solve(prob, BS3(), saveat=3.0) # You can try `Rosenbrock23()` too + sol = solve(prob, BS3(), saveat=10.0) # You can try `Rosenbrock23()` too return sol end From ece09a29f5b12cad19196db9da461d58242385c4 Mon Sep 17 00:00:00 2001 From: 2211567 Date: Tue, 10 Jun 2025 13:44:43 +0200 Subject: [PATCH 03/11] Zebra provisorisch erstellt lol --- scripts/run_simulation.jl | 2 +- src/solver.jl | 91 +++++++++++++++++++++++++++++++++------ 2 files changed, 79 insertions(+), 14 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 6146ea8..0e253a3 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -5,7 +5,7 @@ include("../src/visualization.jl") using .Visualization N = 256 -tspan = (0.0, 5000.0) +tspan = (0.0, 1500.0) sol = AnimalFurFHN.run_simulation(tspan, N) diff --git a/src/solver.jl b/src/solver.jl index 16287aa..a2be899 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -41,7 +41,7 @@ function blocks_ic(N) row_end = min(row_center + 7, N) col_start = max(col_center - 8, 1) col_end = min(col_center + 7, N) - u[row_start:row_end, col_start:col_end] .= -0.534522 + u[row_start:row_end, col_start:col_end] .= -0.01 end safe_block!(u, p, p) @@ -50,8 +50,8 @@ function blocks_ic(N) end function column_ic(N) - u = fill(0.534522, N, N) - v = fill(0.381802, N, N) + u = fill(0.01, N, N) + v = fill(0.99, N, N) col_center = div(N, 2) col_width = 8 # You can adjust this @@ -59,7 +59,43 @@ function column_ic(N) col_start = max(col_center - div(col_width, 2), 1) col_end = min(col_center + div(col_width, 2) - 1, N) - u[col_start:col_end, :] .= -0.534522 + u[col_start:col_end, :] .= -0.01 + + return vec(u), vec(v) +end + +function two_rows_edge_distance_ic(N) + row_width = 8 + distance_from_edge = 50 + u = fill(0.01, N, N) + v = fill(0.99, N, N) + + # --- Input Validation --- + if row_width <= 0 || distance_from_edge < 0 + error("row_width must be positive and distance_from_edge must be non-negative.") + end + + # Calculate column 1 (from the left edge) + col1_start = distance_from_edge + 1 + col1_end = col1_start + row_width - 1 + + # Calculate column 2 (from the right edge) + col2_end = N - distance_from_edge + col2_start = col2_end - row_width + 1 + + # --- Further Validation for placement --- + if col1_end > N || col2_start < 1 + error("Columns go out of bounds. Adjust N, row_width, or distance_from_edge.") + end + if col1_end >= col2_start # Check for overlap or touching + error("Columns overlap or touch. Adjust N, row_width, or distance_from_edge such that 2 * (distance_from_edge + row_width) <= N.") + end + + # Apply the first column + u[:, col1_start:col1_end] .= -0.01 + + # Apply the second column + u[:, col2_start:col2_end] .= -0.01 return vec(u), vec(v) end @@ -79,14 +115,43 @@ function center_band_ic(N) end function circle_ic(N) - u = fill(0.534522, N, N) - v = fill(0.381802, N, N) + u = fill(0.01, N, N) + v = fill(0.99, N, N) cx, cy = div(N, 2), div(N, 2) # center of matrix - radius = 0.250 * N # circle radius = 3/4 of N divided by 2 + radius = 0.125 * N # circle radius = 3/4 of N divided by 2 for i in 1:N, j in 1:N if (i - cx)^2 + (j - cy)^2 ≤ radius^2 - u[i, j] = -0.534522 + u[i, j] = -0.01 + end + end + + return vec(u), vec(v) +end + +function three_circles_random_ic(N) + u = fill(0.01, N, N) + v = fill(0.99, N, N) + radius = 0.125 * N + + # Define the bounds for random centers to ensure the circle stays within the matrix + min_coord = ceil(Int, radius) + 1 + max_coord = floor(Int, N - radius) + + if min_coord > max_coord + error("Matrix size N is too small to place circles of this radius without overlap or going out of bounds.") + end + + for _ in 1:5 # Place 3 circles + # Generate random center coordinates + cx = rand(min_coord:max_coord) + cy = rand(min_coord:max_coord) + + # Apply the circle to the matrix + for i in 1:N, j in 1:N + if (i - cx)^2 + (j - cy)^2 ≤ radius^2 + u[i, j] = -0.01 + end end end @@ -94,8 +159,8 @@ function circle_ic(N) end function squiggle_ic(N, Lx=400.0, Ly=400.0) - uplus = 0.534522 - vplus = 0.381802 + uplus = 0.01 + vplus = 0.99 uminus = -uplus # Create coordinate grids @@ -135,16 +200,16 @@ end """ function run_simulation(tspan::Tuple{Float64,Float64}, N::Int) # Turing-spot parameters - p = FHNParams(N=N, dx=1.0, Du=0.00016, Dv=0.001, ϵ=0.1, a=0.5, b=0.9) + p = FHNParams(N=N, dx=1.0, Du=0.016, Dv=0.1, ϵ=0.1, a=0.5, b=0.9) # Initial conditions (random noise) - Random.seed!(1234) + #Random.seed!(4321) # Create two vectors with length N*N with numbers between 0.1 and 0.11 #u0 = vec(0.4 .+ 0.01 .* rand(N, N)) #v0 = vec(0.4 .+ 0.01 .* rand(N, N)) # Or use this - u0, v0 = squiggle_ic(N) + u0, v0 = two_rows_edge_distance_ic(N) y0 = vcat(vcat(u0, v0)) From f5d3ebd0c2cf84ef52218904deb7b3e25bbb1638 Mon Sep 17 00:00:00 2001 From: 2211567 Date: Tue, 10 Jun 2025 16:43:10 +0200 Subject: [PATCH 04/11] cheetah conditions wip --- scripts/run_simulation.jl | 2 +- src/AnimalFurFHN.jl | 1 + src/solver.jl | 174 +---------------------------------- src/utils.jl | 184 ++++++++++++++++++++++++++++++++++++++ src/visualization.jl | 2 +- 5 files changed, 191 insertions(+), 172 deletions(-) create mode 100644 src/utils.jl diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 0e253a3..b724908 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -5,7 +5,7 @@ include("../src/visualization.jl") using .Visualization N = 256 -tspan = (0.0, 1500.0) +tspan = (1.0, 3000.0) sol = AnimalFurFHN.run_simulation(tspan, N) diff --git a/src/AnimalFurFHN.jl b/src/AnimalFurFHN.jl index 060abd7..f3f77c3 100644 --- a/src/AnimalFurFHN.jl +++ b/src/AnimalFurFHN.jl @@ -3,6 +3,7 @@ module AnimalFurFHN include("constants.jl") include("laplacian.jl") +include("utils.jl") include("solver.jl") export run_simulation # Make sure this is here! diff --git a/src/solver.jl b/src/solver.jl index a2be899..07df5c5 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -30,162 +30,6 @@ function fhn!(du, u, p::FHNParams, t=0) du .= vcat(vec(fu), vec(fv)) end -# helper functions for filling cells in specific places of the matrix -function blocks_ic(N) - u = fill(1.0, N, N) - v = fill(0.0, N, N) - p = div(N, 2) - - function safe_block!(u, row_center, col_center) - row_start = max(row_center - 8, 1) - row_end = min(row_center + 7, N) - col_start = max(col_center - 8, 1) - col_end = min(col_center + 7, N) - u[row_start:row_end, col_start:col_end] .= -0.01 - end - - safe_block!(u, p, p) - - return vec(u), vec(v) -end - -function column_ic(N) - u = fill(0.01, N, N) - v = fill(0.99, N, N) - - col_center = div(N, 2) - col_width = 8 # You can adjust this - - col_start = max(col_center - div(col_width, 2), 1) - col_end = min(col_center + div(col_width, 2) - 1, N) - - u[col_start:col_end, :] .= -0.01 - - return vec(u), vec(v) -end - -function two_rows_edge_distance_ic(N) - row_width = 8 - distance_from_edge = 50 - u = fill(0.01, N, N) - v = fill(0.99, N, N) - - # --- Input Validation --- - if row_width <= 0 || distance_from_edge < 0 - error("row_width must be positive and distance_from_edge must be non-negative.") - end - - # Calculate column 1 (from the left edge) - col1_start = distance_from_edge + 1 - col1_end = col1_start + row_width - 1 - - # Calculate column 2 (from the right edge) - col2_end = N - distance_from_edge - col2_start = col2_end - row_width + 1 - - # --- Further Validation for placement --- - if col1_end > N || col2_start < 1 - error("Columns go out of bounds. Adjust N, row_width, or distance_from_edge.") - end - if col1_end >= col2_start # Check for overlap or touching - error("Columns overlap or touch. Adjust N, row_width, or distance_from_edge such that 2 * (distance_from_edge + row_width) <= N.") - end - - # Apply the first column - u[:, col1_start:col1_end] .= -0.01 - - # Apply the second column - u[:, col2_start:col2_end] .= -0.01 - - return vec(u), vec(v) -end - -function center_band_ic(N) - u = fill(0.0, N, N) - v = fill(0.0, N, N) - - band_width = div(N, 8) - row_start = div(N, 2) - div(band_width, 2) - row_end = div(N, 2) + div(band_width, 2) - - u[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N) - v[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N) - - return vec(u), vec(v) -end - -function circle_ic(N) - u = fill(0.01, N, N) - v = fill(0.99, N, N) - cx, cy = div(N, 2), div(N, 2) # center of matrix - radius = 0.125 * N # circle radius = 3/4 of N divided by 2 - - for i in 1:N, j in 1:N - if (i - cx)^2 + (j - cy)^2 ≤ radius^2 - u[i, j] = -0.01 - end - end - - return vec(u), vec(v) -end - -function three_circles_random_ic(N) - u = fill(0.01, N, N) - v = fill(0.99, N, N) - radius = 0.125 * N - - # Define the bounds for random centers to ensure the circle stays within the matrix - min_coord = ceil(Int, radius) + 1 - max_coord = floor(Int, N - radius) - - if min_coord > max_coord - error("Matrix size N is too small to place circles of this radius without overlap or going out of bounds.") - end - - for _ in 1:5 # Place 3 circles - # Generate random center coordinates - cx = rand(min_coord:max_coord) - cy = rand(min_coord:max_coord) - - # Apply the circle to the matrix - for i in 1:N, j in 1:N - if (i - cx)^2 + (j - cy)^2 ≤ radius^2 - u[i, j] = -0.01 - end - end - end - - return vec(u), vec(v) -end - -function squiggle_ic(N, Lx=400.0, Ly=400.0) - uplus = 0.01 - vplus = 0.99 - uminus = -uplus - - # Create coordinate grids - x = LinRange(0, Lx, N) - y = LinRange(0, Ly, N) - X = repeat(x', N, 1) # Transposed to align with meshgrid X - Y = repeat(y, 1, N) # Broadcasted to align with meshgrid Y - - # Squiggle pattern - cos_term = 0.05 * Lx .* sin.(10 * 2π .* Y ./ Ly .+ π * 0.3) - exp_term = exp.(-((Y .- Ly / 2) ./ (0.1 * Ly)) .^ 2) - width = 0.05 * Lx - Z = exp.(-((X .- Lx / 2 .+ cos_term .* exp_term) ./ width) .^ 2) - - - u = fill(uplus, N, N) - v = fill(vplus, N, N) - - # Apply squiggle - u[Z .> 0.8] .= uminus - - return vec(u), vec(v) -end - - """ run_simulation(tspan::Tuple{Float64,Float64}, N::Int) @@ -199,21 +43,11 @@ end - `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=0.016, Dv=0.1, ϵ=0.1, a=0.5, b=0.9) + # Initial conditions + # params, y0 = zebra_conditions(N) # tspan of ~3500 enough + params, y0 = cheetah_conditions(N) - # Initial conditions (random noise) - #Random.seed!(4321) - # Create two vectors with length N*N with numbers between 0.1 and 0.11 - #u0 = vec(0.4 .+ 0.01 .* rand(N, N)) - #v0 = vec(0.4 .+ 0.01 .* rand(N, N)) - - # Or use this - u0, v0 = two_rows_edge_distance_ic(N) - - y0 = vcat(vcat(u0, v0)) - - prob = ODEProblem(fhn!, y0, tspan, p) + prob = ODEProblem(fhn!, y0, tspan, params) sol = solve(prob, BS3(), saveat=10.0) # You can try `Rosenbrock23()` too return sol diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..cf54588 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,184 @@ + +# initial conditions for different patterns +function zebra_conditions(N) + # Turing-spot parameters + params = FHNParams(N=N, dx=1.0, Du=0.016, Dv=0.1, ϵ=0.1, a=0.5, b=0.9) + + # Or use this + u0, v0 = two_rows_edge_distance_ic(N) + + return params, vcat(u0, v0) +end + +function cheetah_conditions(N) + # Turing-spot parameters + params = FHNParams(N=N, dx=1.0, Du=0.0025, Dv=0.6, ϵ=0.05, a=0.7, b=0.8) + + # Approximate Homogenous Steady State (HSS) for a=0.7, b=0.8 + u_hss = -1.2 + v_hss = -0.625 + + # Small perturbations around the HSS + # rand(N,N) generates values between 0 and 1. + # (2 .* rand(N,N) .- 1) generates values between -1 and 1 symmetrically around 0. + perturbation_amplitude = 0.01 + u0 = vec(u_hss .+ perturbation_amplitude .* (2 .* rand(N, N) .- 1)) + v0 = vec(v_hss .+ perturbation_amplitude .* (2 .* rand(N, N) .- 1)) + + return params, vcat(u0, v0) +end + +# helper functions for filling cells in specific places of the matrix +function blocks_ic(N) + u = fill(1.0, N, N) + v = fill(0.0, N, N) + p = div(N, 2) + + function safe_block!(u, row_center, col_center) + row_start = max(row_center - 8, 1) + row_end = min(row_center + 7, N) + col_start = max(col_center - 8, 1) + col_end = min(col_center + 7, N) + u[row_start:row_end, col_start:col_end] .= -0.01 + end + + safe_block!(u, p, p) + + return vec(u), vec(v) +end + +function column_ic(N) + u = fill(0.01, N, N) + v = fill(0.99, N, N) + + col_center = div(N, 2) + col_width = 8 # You can adjust this + + col_start = max(col_center - div(col_width, 2), 1) + col_end = min(col_center + div(col_width, 2) - 1, N) + + u[col_start:col_end, :] .= -0.01 + + return vec(u), vec(v) +end + +function two_rows_edge_distance_ic(N) + row_width = 8 + distance_from_edge = 50 + u = fill(0.01, N, N) + v = fill(0.99, N, N) + + # --- Input Validation --- + if row_width <= 0 || distance_from_edge < 0 + error("row_width must be positive and distance_from_edge must be non-negative.") + end + + # Calculate column 1 (from the left edge) + col1_start = distance_from_edge + 1 + col1_end = col1_start + row_width - 1 + + # Calculate column 2 (from the right edge) + col2_end = N - distance_from_edge + col2_start = col2_end - row_width + 1 + + # --- Further Validation for placement --- + if col1_end > N || col2_start < 1 + error("Columns go out of bounds. Adjust N, row_width, or distance_from_edge.") + end + if col1_end >= col2_start # Check for overlap or touching + error("Columns overlap or touch. Adjust N, row_width, or distance_from_edge such that 2 * (distance_from_edge + row_width) <= N.") + end + + # Apply the first column + u[:, col1_start:col1_end] .= -0.01 + + # Apply the second column + u[:, col2_start:col2_end] .= -0.01 + + return vec(u), vec(v) +end + +function center_band_ic(N) + u = fill(0.0, N, N) + v = fill(0.0, N, N) + + band_width = div(N, 8) + row_start = div(N, 2) - div(band_width, 2) + row_end = div(N, 2) + div(band_width, 2) + + u[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N) + v[row_start:row_end, :] .= 0.1 .+ 0.01 .* randn(band_width + 1, N) + + return vec(u), vec(v) +end + +function circle_ic(N) + u = fill(0.01, N, N) + v = fill(0.99, N, N) + cx, cy = div(N, 2), div(N, 2) # center of matrix + radius = 0.125 * N # circle radius = 3/4 of N divided by 2 + + for i in 1:N, j in 1:N + if (i - cx)^2 + (j - cy)^2 ≤ radius^2 + u[i, j] = -0.01 + end + end + + return vec(u), vec(v) +end + +function three_circles_random_ic(N) + u = fill(0.01, N, N) + v = fill(0.99, N, N) + radius = 0.125 * N + + # Define the bounds for random centers to ensure the circle stays within the matrix + min_coord = ceil(Int, radius) + 1 + max_coord = floor(Int, N - radius) + + if min_coord > max_coord + error("Matrix size N is too small to place circles of this radius without overlap or going out of bounds.") + end + + for _ in 1:5 # Place 3 circles + # Generate random center coordinates + cx = rand(min_coord:max_coord) + cy = rand(min_coord:max_coord) + + # Apply the circle to the matrix + for i in 1:N, j in 1:N + if (i - cx)^2 + (j - cy)^2 ≤ radius^2 + u[i, j] = -0.01 + end + end + end + + return vec(u), vec(v) +end + +function squiggle_ic(N, Lx=400.0, Ly=400.0) + uplus = 0.01 + vplus = 0.99 + uminus = -uplus + + # Create coordinate grids + x = LinRange(0, Lx, N) + y = LinRange(0, Ly, N) + X = repeat(x', N, 1) # Transposed to align with meshgrid X + Y = repeat(y, 1, N) # Broadcasted to align with meshgrid Y + + # Squiggle pattern + cos_term = 0.05 * Lx .* sin.(10 * 2π .* Y ./ Ly .+ π * 0.3) + exp_term = exp.(-((Y .- Ly / 2) ./ (0.1 * Ly)) .^ 2) + width = 0.05 * Lx + Z = exp.(-((X .- Lx / 2 .+ cos_term .* exp_term) ./ width) .^ 2) + + + u = fill(uplus, N, N) + v = fill(vplus, N, N) + + # Apply squiggle + u[Z .> 0.8] .= uminus + + return vec(u), vec(v) +end diff --git a/src/visualization.jl b/src/visualization.jl index 00a0b71..e73d9f1 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -21,7 +21,7 @@ function step_through_solution(sol, N::Int) # 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=:berlin) + hmap = heatmap!(ax, heat_obs, colormap=:RdGy) # Update heatmap on slider movement on(slider.value) do i From 9223916f326d621b4dc46dbdeeb28950f4917a31 Mon Sep 17 00:00:00 2001 From: 2211567 Date: Wed, 11 Jun 2025 15:31:55 +0200 Subject: [PATCH 05/11] cheetah design more or less done --- scripts/run_simulation.jl | 4 ++-- src/solver.jl | 2 +- src/utils.jl | 25 ++++++++++++++----------- 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index b724908..87de992 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -4,8 +4,8 @@ using .AnimalFurFHN include("../src/visualization.jl") using .Visualization -N = 256 -tspan = (1.0, 3000.0) +N = 100 +tspan = (1.0, 1000.0) sol = AnimalFurFHN.run_simulation(tspan, N) diff --git a/src/solver.jl b/src/solver.jl index 07df5c5..91bacdd 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -45,7 +45,7 @@ end function run_simulation(tspan::Tuple{Float64,Float64}, N::Int) # Initial conditions # params, y0 = zebra_conditions(N) # tspan of ~3500 enough - params, y0 = cheetah_conditions(N) + params, y0 = cell_conditions(N) prob = ODEProblem(fhn!, y0, tspan, params) sol = solve(prob, BS3(), saveat=10.0) # You can try `Rosenbrock23()` too diff --git a/src/utils.jl b/src/utils.jl index cf54588..f79072f 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -12,20 +12,23 @@ end function cheetah_conditions(N) # Turing-spot parameters - params = FHNParams(N=N, dx=1.0, Du=0.0025, Dv=0.6, ϵ=0.05, a=0.7, b=0.8) + #params = FHNParams(N=N, dx=1.0, Du=5e-6, Dv=2e-3, ϵ=0.025, a=0.6, b=0.15) + params = FHNParams(N=N, dx=1.0, Du=0.182, Dv=0.5, ϵ=0.05, a=0.2, b=2.0) - # Approximate Homogenous Steady State (HSS) for a=0.7, b=0.8 - u_hss = -1.2 - v_hss = -0.625 + u0 = 0.05 .* (2 .* rand(N, N) .- 1) # noise in [-0.05, 0.05] + v0 = 0.05 .* (2 .* rand(N, N) .- 1) - # Small perturbations around the HSS - # rand(N,N) generates values between 0 and 1. - # (2 .* rand(N,N) .- 1) generates values between -1 and 1 symmetrically around 0. - perturbation_amplitude = 0.01 - u0 = vec(u_hss .+ perturbation_amplitude .* (2 .* rand(N, N) .- 1)) - v0 = vec(v_hss .+ perturbation_amplitude .* (2 .* rand(N, N) .- 1)) + return params, vcat(vec(u0), vec(v0)) +end - return params, vcat(u0, v0) +function cell_conditions(N) + # Turing-spot parameters + params = FHNParams(N=N, dx=1.0, Du=5e-5, Dv=6e-3, ϵ=0.05, a=0.6, b=0.1) + + u0 = 0.534522 .* (2 .* rand(N, N) .- 1) # noise in [-0.05, 0.05] + v0 = 0.381802 .* (2 .* rand(N, N) .- 1) + + return params, vcat(vec(u0), vec(v0)) end # helper functions for filling cells in specific places of the matrix From 0dd6062c5d18008ed58dca8e5dacb6153c8baaef Mon Sep 17 00:00:00 2001 From: 2211567 Date: Wed, 11 Jun 2025 18:13:12 +0200 Subject: [PATCH 06/11] coral pattern kinda done --- scripts/run_simulation.jl | 4 ++-- src/solver.jl | 4 ++-- src/utils.jl | 21 ++++++++++++++++----- 3 files changed, 20 insertions(+), 9 deletions(-) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 87de992..2d5ba7e 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -4,8 +4,8 @@ using .AnimalFurFHN include("../src/visualization.jl") using .Visualization -N = 100 -tspan = (1.0, 1000.0) +N = 128 +tspan = (0.0, 2000.0) sol = AnimalFurFHN.run_simulation(tspan, N) diff --git a/src/solver.jl b/src/solver.jl index 91bacdd..b99226c 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -45,10 +45,10 @@ end function run_simulation(tspan::Tuple{Float64,Float64}, N::Int) # Initial conditions # params, y0 = zebra_conditions(N) # tspan of ~3500 enough - params, y0 = cell_conditions(N) + params, y0 = coral_conditions(N) prob = ODEProblem(fhn!, y0, tspan, params) - sol = solve(prob, BS3(), saveat=10.0) # You can try `Rosenbrock23()` too + sol = solve(prob, BS3(), saveat=50.0) # You can try `Rosenbrock23()` too return sol end diff --git a/src/utils.jl b/src/utils.jl index f79072f..7298bf4 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -21,14 +21,13 @@ function cheetah_conditions(N) return params, vcat(vec(u0), vec(v0)) end -function cell_conditions(N) +function coral_conditions(N) # Turing-spot parameters - params = FHNParams(N=N, dx=1.0, Du=5e-5, Dv=6e-3, ϵ=0.05, a=0.6, b=0.1) + params = FHNParams(N=N, dx=1.0, Du=0.001, Dv=0.06, ϵ=0.05, a=0.0, b=1.2) - u0 = 0.534522 .* (2 .* rand(N, N) .- 1) # noise in [-0.05, 0.05] - v0 = 0.381802 .* (2 .* rand(N, N) .- 1) + u0, v0 = coral_ic(N) - return params, vcat(vec(u0), vec(v0)) + return params, vcat(u0, v0) end # helper functions for filling cells in specific places of the matrix @@ -185,3 +184,15 @@ function squiggle_ic(N, Lx=400.0, Ly=400.0) return vec(u), vec(v) end + +function coral_ic(N) + u = fill(0.534522, N, N) + v = fill(0.381802, N, N) + + for _ in 1:40 # place 15 noisy seeds + i, j = rand(10:N-10), rand(10:N-10) + u[i-2:i+2, j-2:j+2] .= -0.534522 .+ 0.2 * rand(5, 5) + end + + return vec(u), vec(v) +end From d5ea440656b2ec2ec767b29643ba88436890a7c8 Mon Sep 17 00:00:00 2001 From: 2211567 Date: Sat, 14 Jun 2025 13:23:14 +0200 Subject: [PATCH 07/11] build_ui also more modular; adjusted includes and namespaces --- Manifest.toml | 180 +++++++++++++++------------ Project.toml | 1 + scripts/main.jl | 25 ++-- scripts/run_simulation.jl | 2 +- scripts/show_frame.jl | 21 ---- src/AnimalFurFHN.jl | 10 -- src/{solver.jl => fhn_solver.jl} | 3 +- src/gray_scott_solver.jl | 5 +- src/utils/constants.jl | 8 +- src/{utils.jl => utils/templates.jl} | 0 src/visualization.jl | 26 ++-- 11 files changed, 141 insertions(+), 140 deletions(-) delete mode 100644 scripts/show_frame.jl delete mode 100644 src/AnimalFurFHN.jl rename src/{solver.jl => fhn_solver.jl} (96%) rename src/{utils.jl => utils/templates.jl} (100%) diff --git a/Manifest.toml b/Manifest.toml index df84095..6c62bc0 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.4" +julia_version = "1.11.5" manifest_format = "2.0" -project_hash = "b5f5f0b50b1e0b7dd05bb93e0b1bb03fea235d53" +project_hash = "5fd84347cd356de5b819aa3ea793fa63c45180e4" [[deps.ADTypes]] git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32" @@ -194,9 +194,9 @@ version = "0.1.6" [[deps.BoundaryValueDiffEq]] deps = ["ADTypes", "BoundaryValueDiffEqAscher", "BoundaryValueDiffEqCore", "BoundaryValueDiffEqFIRK", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqMIRKN", "BoundaryValueDiffEqShooting", "DiffEqBase", "FastClosures", "ForwardDiff", "LinearAlgebra", "Reexport", "SciMLBase"] -git-tree-sha1 = "6606f3f7d43b038a8c34ee4cdc31bbd93b447407" +git-tree-sha1 = "d6ec33e4516b2e790a64128afdb54f3b536667a7" uuid = "764a87c0-6b3e-53db-9096-fe964310641d" -version = "5.17.0" +version = "5.18.0" [deps.BoundaryValueDiffEq.extensions] BoundaryValueDiffEqODEInterfaceExt = "ODEInterface" @@ -206,39 +206,39 @@ version = "5.17.0" [[deps.BoundaryValueDiffEqAscher]] deps = ["ADTypes", "AlmostBlockDiagonals", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield"] -git-tree-sha1 = "add7560694794f30eec7803b9c4036efb96bf493" +git-tree-sha1 = "64a777e06d995f677c86c7ddbb85f393074a0877" uuid = "7227322d-7511-4e07-9247-ad6ff830280e" -version = "1.6.0" +version = "1.7.0" [[deps.BoundaryValueDiffEqCore]] deps = ["ADTypes", "Adapt", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "ForwardDiff", "LineSearch", "LinearAlgebra", "Logging", "NonlinearSolveFirstOrder", "PreallocationTools", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays", "SparseConnectivityTracer", "SparseMatrixColorings"] -git-tree-sha1 = "1e97d81b6ea5e8e938c57a3c640a5cff1e181a2f" +git-tree-sha1 = "c83bf97da90dd379b1e3f4d9c6f3d0ae48eb0b29" uuid = "56b672f2-a5fe-4263-ab2d-da677488eb3a" -version = "1.9.0" +version = "1.10.0" [[deps.BoundaryValueDiffEqFIRK]] deps = ["ADTypes", "ArrayInterface", "BandedMatrices", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastAlmostBandedMatrices", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays"] -git-tree-sha1 = "525c65da4b46271b4b5fc2178ccde16c99c21a41" +git-tree-sha1 = "e58ee9acfc6dce6dcc368fc0483952bec5625513" uuid = "85d9eb09-370e-4000-bb32-543851f73618" -version = "1.7.0" +version = "1.8.1" [[deps.BoundaryValueDiffEqMIRK]] deps = ["ADTypes", "ArrayInterface", "BandedMatrices", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastAlmostBandedMatrices", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays"] -git-tree-sha1 = "f729bdedaedb537bf7883c9f71e7577e1c7a07f6" +git-tree-sha1 = "43debeee94167e2dc744f4a385213c4f0d16b4c3" uuid = "1a22d4ce-7765-49ea-b6f2-13c8438986a6" -version = "1.7.0" +version = "1.8.1" [[deps.BoundaryValueDiffEqMIRKN]] deps = ["ADTypes", "ArrayInterface", "BandedMatrices", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastAlmostBandedMatrices", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays"] -git-tree-sha1 = "dfde7afe34137182a4bb8cd98f01d870772be999" +git-tree-sha1 = "1d92c9f7567b627514e143a3caf93af6d235c2db" uuid = "9255f1d6-53bf-473e-b6bd-23f1ff009da4" -version = "1.6.0" +version = "1.7.1" [[deps.BoundaryValueDiffEqShooting]] deps = ["ADTypes", "ArrayInterface", "BandedMatrices", "BoundaryValueDiffEqCore", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastAlmostBandedMatrices", "FastClosures", "ForwardDiff", "LinearAlgebra", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "SparseArrays"] -git-tree-sha1 = "255c15c2d262ea8d3c5db1e0f130f829362c0fd4" +git-tree-sha1 = "1460c449c91c9bc4c8fb08fb5e878811ab38d667" uuid = "ed55bfe0-3725-4db6-871e-a1dc9f42a757" -version = "1.7.0" +version = "1.8.1" [[deps.BracketingNonlinearSolve]] deps = ["CommonSolve", "ConcreteStructs", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase"] @@ -390,9 +390,9 @@ uuid = "2569d6c7-a4a2-43d3-a901-331e8e4be471" version = "0.2.3" [[deps.ConstructionBase]] -git-tree-sha1 = "76219f1ed5771adbb096743bff43fb5fdd4c1157" +git-tree-sha1 = "b4b092499347b18a015186eae3042f72267106cb" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" -version = "1.5.8" +version = "1.6.0" weakdeps = ["IntervalSets", "LinearAlgebra", "StaticArrays"] [deps.ConstructionBase.extensions] @@ -457,9 +457,9 @@ version = "5.53.1" [[deps.DiffEqBase]] deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "TruncatedStacktraces"] -git-tree-sha1 = "a0e5b5669df9465bc3dd32ea4a8ddeefbc0f7b5c" +git-tree-sha1 = "2d87d7bd165c1ca0d11923a9fabe90a9d71e88a6" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.175.0" +version = "6.176.0" [deps.DiffEqBase.extensions] DiffEqBaseCUDAExt = "CUDA" @@ -533,9 +533,9 @@ version = "7.16.1" [[deps.DifferentiationInterface]] deps = ["ADTypes", "LinearAlgebra"] -git-tree-sha1 = "c8d85ecfcbaef899308706bebdd8b00107f3fb43" +git-tree-sha1 = "210933c93f39f832d92f9efbbe69a49c453db36d" uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" -version = "0.6.54" +version = "0.7.1" [deps.DifferentiationInterface.extensions] DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore" @@ -614,9 +614,9 @@ version = "0.25.120" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[deps.DocStringExtensions]] -git-tree-sha1 = "e7b7e6f178525d17c720ab9c081e4ef04429f860" +git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.4" +version = "0.9.5" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] @@ -635,9 +635,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" version = "1.0.5" [[deps.EnzymeCore]] -git-tree-sha1 = "1eb59f40a772d0fbd4cb75e00b3fa7f5f79c975a" +git-tree-sha1 = "7d7822a643c33bbff4eab9c87ca8459d7c688db0" uuid = "f151be2c-9106-41f4-ab19-57ee4f262869" -version = "0.8.9" +version = "0.8.11" weakdeps = ["Adapt"] [deps.EnzymeCore.extensions] @@ -694,9 +694,9 @@ version = "6.1.2+0" [[deps.FFTW]] deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] -git-tree-sha1 = "7de7c78d681078f027389e067864a8d53bd7c3c9" +git-tree-sha1 = "797762812ed063b9b94f6cc7742bc8883bb5e69e" uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "1.8.1" +version = "1.9.0" [[deps.FFTW_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -852,9 +852,9 @@ version = "2.13.4+0" [[deps.FreeTypeAbstraction]] deps = ["BaseDirs", "ColorVectorSpace", "Colors", "FreeType", "GeometryBasics", "Mmap"] -git-tree-sha1 = "eaca92bac73aa42f68c57d1b8df1b746eeb2bdaa" +git-tree-sha1 = "4ebb930ef4a43817991ba35db6317a05e59abd11" uuid = "663a7486-cb36-511b-a19d-713bb74d65c9" -version = "0.10.7" +version = "0.10.8" [[deps.FriBidi_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -898,9 +898,9 @@ version = "3.4.0+2" [[deps.GLMakie]] deps = ["ColorTypes", "Colors", "FileIO", "FixedPointNumbers", "FreeTypeAbstraction", "GLFW", "GeometryBasics", "LinearAlgebra", "Makie", "Markdown", "MeshIO", "ModernGL", "Observables", "PrecompileTools", "Printf", "ShaderAbstractions", "StaticArrays"] -git-tree-sha1 = "0650df73822ce8808dc473e3e1c7581bdb367083" +git-tree-sha1 = "aacae8d91f7da30979e189e8046c5bd655421687" uuid = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" -version = "0.11.8" +version = "0.12.0" [[deps.GPUArraysCore]] deps = ["Adapt"] @@ -956,10 +956,10 @@ uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" version = "1.3.15+0" [[deps.Graphs]] -deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] -git-tree-sha1 = "3169fd3440a02f35e549728b0890904cfd4ae58a" +deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "c5abfa0ae0aaee162a3fbb053c13ecda39be545b" uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" -version = "1.12.1" +version = "1.13.0" [[deps.GridLayoutBase]] deps = ["GeometryBasics", "InteractiveUtils", "Observables"] @@ -1145,9 +1145,9 @@ version = "3.1.1+0" [[deps.JumpProcesses]] deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FunctionWrappers", "Graphs", "LinearAlgebra", "Markdown", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "SymbolicIndexingInterface", "UnPack"] -git-tree-sha1 = "216c196df09c8b80a40a2befcb95760eb979bcfd" +git-tree-sha1 = "fb7fd516de38db80f50fe15e57d44da2836365e7" uuid = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" -version = "9.15.0" +version = "9.16.0" weakdeps = ["FastBroadcast"] [[deps.KernelDensity]] @@ -1308,9 +1308,9 @@ weakdeps = ["LineSearches"] [[deps.LineSearches]] deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] -git-tree-sha1 = "e4c3be53733db1051cc15ecf573b1042b3a712a1" +git-tree-sha1 = "4adee99b7262ad2a1a4bbbc59d993d24e55ea96f" uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" -version = "7.3.0" +version = "7.4.0" [[deps.LinearAlgebra]] deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] @@ -1319,9 +1319,9 @@ version = "1.11.0" [[deps.LinearSolve]] deps = ["ArrayInterface", "ChainRulesCore", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "LazyArrays", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "StaticArraysCore", "UnPack"] -git-tree-sha1 = "c618a6a774d5712c6bf02dbcceb51b6dc6b9bb89" +git-tree-sha1 = "c0d1a91a50af6778863d320761f807f641f74935" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -version = "3.16.0" +version = "3.17.0" [deps.LinearSolve.extensions] LinearSolveBandedMatricesExt = "BandedMatrices" @@ -1392,15 +1392,15 @@ version = "0.5.16" [[deps.Makie]] deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Dates", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageBase", "ImageIO", "InteractiveUtils", "Interpolations", "IntervalSets", "InverseFunctions", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "PNGFiles", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun", "Unitful"] -git-tree-sha1 = "968f03dc65c8144a728f988051a88c78fe625e26" +git-tree-sha1 = "29c5cf33287afa8c06711fec1b33d074f88dd354" uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" -version = "0.22.7" +version = "0.23.0" [[deps.MakieCore]] deps = ["ColorTypes", "GeometryBasics", "IntervalSets", "Observables"] -git-tree-sha1 = "733d910c70805e7114c82508bae99c6cdf004466" +git-tree-sha1 = "6ec4aba99ce1a8c247822ce346b273d0d686e141" uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b" -version = "0.9.3" +version = "0.10.0" [[deps.ManualMemory]] git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd" @@ -1493,9 +1493,9 @@ version = "0.2.4" [[deps.NLSolversBase]] deps = ["ADTypes", "DifferentiationInterface", "Distributed", "FiniteDiff", "ForwardDiff"] -git-tree-sha1 = "b14c7be6046e7d48e9063a0053f95ee0fc954176" +git-tree-sha1 = "25a6638571a902ecfb1ae2a18fc1575f86b1d4df" uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" -version = "7.9.1" +version = "7.10.0" [[deps.NLsolve]] deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"] @@ -1639,7 +1639,7 @@ version = "3.2.4+0" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.1+4" +version = "0.8.5+0" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -1655,9 +1655,9 @@ version = "0.5.6+0" [[deps.Optim]] deps = ["Compat", "EnumX", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] -git-tree-sha1 = "31b3b1b8e83ef9f1d50d74f1dd5f19a37a304a1f" +git-tree-sha1 = "61942645c38dd2b5b78e2082c9b51ab315315d10" uuid = "429524aa-4258-5aef-a3af-852621145aeb" -version = "1.12.0" +version = "1.13.2" [deps.Optim.extensions] OptimMOIExt = "MathOptInterface" @@ -1665,18 +1665,40 @@ version = "1.12.0" [deps.Optim.weakdeps] MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -[[deps.Opus_jll]]using .Laplacian +[[deps.Opus_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6703a85cb3781bd5909d48730a67205f3f31a575" +uuid = "91d4177d-7536-5919-b921-800302f37372" +version = "1.3.3+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "05868e21324cede2207c6f0f466b4bfef6d5e7ee" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.8.1" + +[[deps.OrdinaryDiffEq]] +deps = ["ADTypes", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "ExponentialUtilities", "FastBroadcast", "FastClosures", "FillArrays", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "InteractiveUtils", "LineSearches", "LinearAlgebra", "LinearSolve", "Logging", "MacroTools", "MuladdMacro", "NonlinearSolve", "OrdinaryDiffEqAdamsBashforthMoulton", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqExplicitRK", "OrdinaryDiffEqExponentialRK", "OrdinaryDiffEqExtrapolation", "OrdinaryDiffEqFIRK", "OrdinaryDiffEqFeagin", "OrdinaryDiffEqFunctionMap", "OrdinaryDiffEqHighOrderRK", "OrdinaryDiffEqIMEXMultistep", "OrdinaryDiffEqLinear", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqLowStorageRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqNordsieck", "OrdinaryDiffEqPDIRK", "OrdinaryDiffEqPRK", "OrdinaryDiffEqQPRK", "OrdinaryDiffEqRKN", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqStabilizedIRK", "OrdinaryDiffEqStabilizedRK", "OrdinaryDiffEqSymplecticRK", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "Polyester", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleNonlinearSolve", "SimpleUnPack", "SparseArrays", "Static", "StaticArrayInterface", "StaticArrays", "TruncatedStacktraces"] +git-tree-sha1 = "1c2b2df870944e0dc01454fd87479847c55fa26c" +uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +version = "6.98.0" + +[[deps.OrdinaryDiffEqAdamsBashforthMoulton]] +deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "Polyester", "RecursiveArrayTools", "Reexport", "Static"] +git-tree-sha1 = "82f78099ecf4e0fa53545811318520d87e7fe0b8" +uuid = "89bda076-bce5-4f1c-845f-551c83cdda9a" +version = "1.2.0" + [[deps.OrdinaryDiffEqBDF]] deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqSDIRK", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "StaticArrays", "TruncatedStacktraces"] -git-tree-sha1 = "42755bd13fe56e9d9ce1bc005f8b206a6b56b731" +git-tree-sha1 = "9124a686af119063bb4d3a8f87044a8f312fcad9" uuid = "6ad6398a-0878-4a85-9266-38940aa047c8" -version = "1.5.1" +version = "1.6.0" [[deps.OrdinaryDiffEqCore]] deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "FastBroadcast", "FastClosures", "FastPower", "FillArrays", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleUnPack", "Static", "StaticArrayInterface", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"] -git-tree-sha1 = "d29adfeb720dd7c251b216d91c4bd4fe67c087df" +git-tree-sha1 = "08dac9c6672a4548439048089bac293759a897fd" uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" -version = "1.26.0" +version = "1.26.1" weakdeps = ["EnzymeCore"] [deps.OrdinaryDiffEqCore.extensions] @@ -1690,9 +1712,9 @@ version = "1.4.0" [[deps.OrdinaryDiffEqDifferentiation]] deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseMatrixColorings", "StaticArrayInterface", "StaticArrays"] -git-tree-sha1 = "c78060115fa4ea9d70ac47fa49496acbc630aefa" +git-tree-sha1 = "efecf0c4cc44e16251b0e718f08b0876b2a82b80" uuid = "4302a76b-040a-498a-8c04-15b101fed76b" -version = "1.9.1" +version = "1.10.0" [[deps.OrdinaryDiffEqExplicitRK]] deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "TruncatedStacktraces"] @@ -1798,9 +1820,9 @@ version = "1.1.0" [[deps.OrdinaryDiffEqRosenbrock]] deps = ["ADTypes", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "Static"] -git-tree-sha1 = "063e5ff1447b3869856ed264b6dcbb21e6e8bdb0" +git-tree-sha1 = "1ce0096d920e95773220e818f29bf4b37ea2bb78" uuid = "43230ef6-c299-4910-a778-202eb28ce4ce" -version = "1.10.1" +version = "1.11.0" [[deps.OrdinaryDiffEqSDIRK]] deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"] @@ -2147,9 +2169,9 @@ version = "0.1.0" [[deps.SciMLBase]] deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "Moshi", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"] -git-tree-sha1 = "ce947672206f6a3a2bee1017c690cfd5fd82d897" +git-tree-sha1 = "04bbcdc8d1f7d6f667f75fbcc68728231e21fabe" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -version = "2.96.0" +version = "2.101.0" [deps.SciMLBase.extensions] SciMLBaseChainRulesCoreExt = "ChainRulesCore" @@ -2180,9 +2202,9 @@ version = "0.1.6" [[deps.SciMLOperators]] deps = ["Accessors", "ArrayInterface", "DocStringExtensions", "LinearAlgebra", "MacroTools"] -git-tree-sha1 = "85608e4aaf758547ffc4030c908318b432114ec9" +git-tree-sha1 = "3249fe77f322fe539e935ecb388c8290cd38a3fc" uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" -version = "1.3.0" +version = "1.3.1" weakdeps = ["SparseArrays", "StaticArraysCore"] [deps.SciMLOperators.extensions] @@ -2286,9 +2308,9 @@ version = "1.11.0" [[deps.SparseConnectivityTracer]] deps = ["ADTypes", "DocStringExtensions", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "fadb2d7010dd92912e5eb31a493613ad4b8c9583" +git-tree-sha1 = "affde0bfd920cfcaa0944d3c0eb3a573fa9c4d1e" uuid = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" -version = "0.6.18" +version = "0.6.20" [deps.SparseConnectivityTracer.extensions] SparseConnectivityTracerDataInterpolationsExt = "DataInterpolations" @@ -2495,9 +2517,9 @@ version = "1.0.1" [[deps.Tables]] deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"] -git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297" +git-tree-sha1 = "f2c1efbc8f3a609aadf318094f8fc5204bdaf344" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.12.0" +version = "1.12.1" [[deps.Tar]] deps = ["ArgTools", "SHA"] @@ -2517,15 +2539,15 @@ version = "1.11.0" [[deps.ThreadingUtilities]] deps = ["ManualMemory"] -git-tree-sha1 = "2d529b6b22791f3e22e7ec5c60b9016e78f5f6bf" +git-tree-sha1 = "d969183d3d244b6c33796b5ed01ab97328f2db85" uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5" -version = "0.5.4" +version = "0.5.5" [[deps.TiffImages]] -deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "ProgressMeter", "SIMD", "UUIDs"] -git-tree-sha1 = "f21231b166166bebc73b99cea236071eb047525b" +deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "PrecompileTools", "ProgressMeter", "SIMD", "UUIDs"] +git-tree-sha1 = "02aca429c9885d1109e58f400c333521c13d48a0" uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69" -version = "0.11.3" +version = "0.11.4" [[deps.TimerOutputs]] deps = ["ExprTools", "Printf"] @@ -2577,14 +2599,16 @@ version = "0.4.1" [[deps.Unitful]] deps = ["Dates", "LinearAlgebra", "Random"] -git-tree-sha1 = "d62610ec45e4efeabf7032d67de2ffdea8344bed" +git-tree-sha1 = "d2282232f8a4d71f79e85dc4dd45e5b12a6297fb" uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" -version = "1.22.1" -weakdeps = ["ConstructionBase", "InverseFunctions"] +version = "1.23.1" +weakdeps = ["ConstructionBase", "ForwardDiff", "InverseFunctions", "Printf"] [deps.Unitful.extensions] ConstructionBaseUnitfulExt = "ConstructionBase" + ForwardDiffExt = "ForwardDiff" InverseFunctionsUnitfulExt = "InverseFunctions" + PrintfExt = "Printf" [[deps.Wayland_jll]] deps = ["Artifacts", "EpollShim_jll", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "XML2_jll"] @@ -2593,10 +2617,10 @@ uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89" version = "1.23.1+0" [[deps.Wayland_protocols_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "5db3e9d307d32baba7067b13fc7b5aa6edd4a19a" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "54b8a029ac145ebe8299463447fd1590b2b1d92f" uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91" -version = "1.36.0+0" +version = "1.44.0+0" [[deps.WebP]] deps = ["CEnum", "ColorTypes", "FileIO", "FixedPointNumbers", "ImageCore", "libwebp_jll"] diff --git a/Project.toml b/Project.toml index 0672cef..a3758f0 100644 --- a/Project.toml +++ b/Project.toml @@ -4,3 +4,4 @@ GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Observables = "510215fc-4207-5dde-b226-833fc4488ee2" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" diff --git a/scripts/main.jl b/scripts/main.jl index e38e373..69ce53f 100644 --- a/scripts/main.jl +++ b/scripts/main.jl @@ -1,28 +1,33 @@ -include("../src/gray_scott_solver.jl") include("../src/visualization.jl") -include("../src/utils/constants.jl") - +include("../src/gray_scott_solver.jl") using Observables using GLMakie -using .GrayScottSolver -using .Visualization -using .Constants -const N = 128 -const dx = 1.0 +using .Constants +using .Visualization + +N = 128 +dx = 1.0 Du, Dv = Observable(0.16), Observable(0.08) F, k = Observable(0.060), Observable(0.062) +param_observables = ( + Du=Du, + Dv=Dv, + F=F, + k=k, +) + 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) +fig = Visualization.build_ui(U, V, param_observables, params_obs, heat_obs) display(fig) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl index 2d5ba7e..3a5a5e5 100644 --- a/scripts/run_simulation.jl +++ b/scripts/run_simulation.jl @@ -5,7 +5,7 @@ include("../src/visualization.jl") using .Visualization N = 128 -tspan = (0.0, 2000.0) +tspan = (0.0, 1000.0) sol = AnimalFurFHN.run_simulation(tspan, N) diff --git a/scripts/show_frame.jl b/scripts/show_frame.jl deleted file mode 100644 index e645396..0000000 --- a/scripts/show_frame.jl +++ /dev/null @@ -1,21 +0,0 @@ -using GLMakie -include("../src/AnimalFurFHN.jl") -using .AnimalFurFHN - -# === Parameters === -N = 100 -tspan = (0.0, 100.0) -frame_index = 1 # final time step - -# === Run simulation === -sol = run_simulation(N, tspan) - -# === Extract activator u === -u = sol.u[frame_index][1:N^2] -u_mat = reshape(u, N, N) - -# === Plot === -fig = Figure() -ax = Axis(fig[1, 1]) -heatmap!(ax, u_mat, colormap=:viridis) -fig diff --git a/src/AnimalFurFHN.jl b/src/AnimalFurFHN.jl deleted file mode 100644 index f3f77c3..0000000 --- a/src/AnimalFurFHN.jl +++ /dev/null @@ -1,10 +0,0 @@ -# src/AnimalFurFHN.jl -module AnimalFurFHN - -include("constants.jl") -include("laplacian.jl") -include("utils.jl") -include("solver.jl") - -export run_simulation # Make sure this is here! -end diff --git a/src/solver.jl b/src/fhn_solver.jl similarity index 96% rename from src/solver.jl rename to src/fhn_solver.jl index b99226c..41e2524 100644 --- a/src/solver.jl +++ b/src/fhn_solver.jl @@ -1,5 +1,6 @@ using DifferentialEquations using Random +using .Laplacian using .Constants """ @@ -45,7 +46,7 @@ end function run_simulation(tspan::Tuple{Float64,Float64}, N::Int) # Initial conditions # params, y0 = zebra_conditions(N) # tspan of ~3500 enough - params, y0 = coral_conditions(N) + params, y0 = zebra_conditions(N) prob = ODEProblem(fhn!, y0, tspan, params) sol = solve(prob, BS3(), saveat=50.0) # You can try `Rosenbrock23()` too diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl index 3668d26..2c213f0 100644 --- a/src/gray_scott_solver.jl +++ b/src/gray_scott_solver.jl @@ -1,10 +1,8 @@ -module GrayScottSolver -using Observables include("utils/constants.jl") include("utils/laplacian.jl") +using Observables using .Constants using .Laplacian -export step!, multi_step! function step!(U, V, params_obs::Observable; dx=1) lap_u = laplacian5(U, dx) @@ -43,4 +41,3 @@ function multi_step!(state, n_steps, heat_obs::Observable, params_obs::Observabl heat_obs[] = step!(state[1], state[2], params_obs; dx=1) end end -end \ No newline at end of file diff --git a/src/utils/constants.jl b/src/utils/constants.jl index c06f46c..a07f4f0 100644 --- a/src/utils/constants.jl +++ b/src/utils/constants.jl @@ -1,6 +1,8 @@ module Constants -struct FHNParams +abstract type PDEParams end + +struct FHNParams <: PDEParams N::Int dx::Float64 # grid spacing Du::Float64 @@ -15,7 +17,7 @@ struct FHNParams new(N, dx, Du, Dv, ϵ, a, b) end -struct GSParams +struct GSParams <: PDEParams N::Int # grid size dx::Float64 # grid spacing Du::Float64 # diffusion rate U @@ -25,6 +27,6 @@ struct GSParams end -export FHNParams, GSParams +export PDEParams, FHNParams, GSParams end # module Constants diff --git a/src/utils.jl b/src/utils/templates.jl similarity index 100% rename from src/utils.jl rename to src/utils/templates.jl diff --git a/src/visualization.jl b/src/visualization.jl index 51c8ccc..e50ceac 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -1,7 +1,8 @@ module Visualization include("gray_scott_solver.jl") + using GLMakie, Observables, Makie -using .GrayScottSolver + """ step_through_solution(sol::SolutionType, N::Int) @@ -19,16 +20,11 @@ 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=:RdGy) - 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) @@ -68,7 +64,7 @@ function param_box!(grid, row, labeltxt, observable::Observable) end end -function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs) +function build_ui(U, V, param_obs_map::NamedTuple, params_obs, heat_obs) reset!(U, V, heat_obs) fig = Figure(size=(800, 800)) gh = GridLayout(fig[1, 1]) @@ -95,10 +91,17 @@ function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs) 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) + param_box!(textboxgrid, 1, "Du", param_obs_map.Du) + param_box!(textboxgrid, 2, "Dv", param_obs_map.Dv) + + if haskey(param_obs_map, :F) + param_box!(textboxgrid, 3, "Feed", param_obs_map.F) + param_box!(textboxgrid, 4, "kill", param_obs_map.k) + elseif haskey(param_obs_map, :ϵ) + param_box!(textboxgrid, 3, "ϵ", param_obs_map.ϵ) + param_box!(textboxgrid, 4, "a", param_obs_map.a) + param_box!(textboxgrid, 5, "b", param_obs_map.b) + end # Timer and state for animation running = Observable(false) @@ -159,6 +162,5 @@ function build_ui(U, V, Du, Dv, F, k, params_obs, heat_obs) return fig end - export step_through_solution, build_ui end From 884e87383cf50fa2acaf988462a635a3e4890b59 Mon Sep 17 00:00:00 2001 From: 2211567 Date: Sat, 14 Jun 2025 14:36:31 +0200 Subject: [PATCH 08/11] compatibility of new UI with FHN implemented --- scripts/main.jl | 31 ++++++++++++++++++++--- scripts/run_simulation.jl | 13 ---------- src/fhn_solver.jl | 46 ++++++++++++++++++++------------- src/gray_scott_solver.jl | 2 ++ src/visualization.jl | 53 +++++++++++++-------------------------- 5 files changed, 76 insertions(+), 69 deletions(-) delete mode 100644 scripts/run_simulation.jl diff --git a/scripts/main.jl b/scripts/main.jl index 69ce53f..a95c179 100644 --- a/scripts/main.jl +++ b/scripts/main.jl @@ -1,5 +1,7 @@ +include("../src/utils/constants.jl") +include("../src/fhn_solver.jl") +#include("../src/gray_scott_solver.jl") include("../src/visualization.jl") -include("../src/gray_scott_solver.jl") using Observables using GLMakie @@ -7,6 +9,8 @@ using GLMakie using .Constants using .Visualization +""" +# GSParams N = 128 dx = 1.0 Du, Dv = Observable(0.16), Observable(0.08) @@ -18,11 +22,32 @@ param_observables = ( F=F, k=k, ) +""" -params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[])) +# FHNParams +N = 128 +dx = 1.0 +Du, Dv = Observable(0.016), Observable(0.1) +ϵ, a, b = Observable(0.1), Observable(0.5), Observable(0.9) +param_observables = ( + Du=Du, + Dv=Dv, + ϵ=ϵ, + a=a, + b=b +) +#params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[])) +params_obs = Observable(Constants.FHNParams(N=N, dx=dx, Du=Du[], Dv=Dv[], ϵ=ϵ[], a=a[], b=b[])) + +""" lift(Du, Dv, F, k) do u, v, f, ki - params_obs[] = GSParams(N, dx, u, v, f, ki) + params_obs[] = GSParams(N, dx, u, v, f, ki) +end +""" + +lift(Du, Dv, ϵ, a, b) do d_u, d_v, eps, aa, bb + params_obs[] = Constants.FHNParams(N=N, dx=dx, Du=d_u, Dv=d_v, ϵ=eps, a=aa, b=bb) end U = ones(N, N) diff --git a/scripts/run_simulation.jl b/scripts/run_simulation.jl deleted file mode 100644 index 3a5a5e5..0000000 --- a/scripts/run_simulation.jl +++ /dev/null @@ -1,13 +0,0 @@ -include("../src/AnimalFurFHN.jl") # this loads the module code -using .AnimalFurFHN -# include optional visualizer only if needed: -include("../src/visualization.jl") -using .Visualization - -N = 128 -tspan = (0.0, 1000.0) - -sol = AnimalFurFHN.run_simulation(tspan, N) - - -Visualization.step_through_solution(sol, N) \ No newline at end of file diff --git a/src/fhn_solver.jl b/src/fhn_solver.jl index 41e2524..4118c01 100644 --- a/src/fhn_solver.jl +++ b/src/fhn_solver.jl @@ -1,7 +1,11 @@ +include("utils/laplacian.jl") + using DifferentialEquations using Random +using Observables + +using ..Constants using .Laplacian -using .Constants """ fhn(du, u, p:FHNParams, t:) @@ -18,7 +22,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) 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 @@ -31,25 +35,31 @@ function fhn!(du, u, p::FHNParams, t=0) du .= vcat(vec(fu), vec(fv)) end -""" - run_simulation(tspan::Tuple{Float64,Float64}, N::Int) +function step!(U, V, params_obs::Observable; dx=1) + p = params_obs[] - solving the ODE and modelling it after FHN + # Flatten initial condition (activation u, recovery v) + u0 = vec(U) + v0 = vec(V) + u_init = vcat(u0, v0) - # Arguments - - `tspan`: tuple of two Float64's representing start and end times for simulation - - `N`: size of the N×N grid + # Define one integration step using your fhn! function + prob = ODEProblem((du, u, p, t) -> fhn!(du, u, p, t), u_init, (0.0, 0.1), p) + sol = solve(prob, Tsit5(); dt=0.1, save_start=false, saveat=10.0) - # Returns - - `sol`: solved differential equation (ODE) -""" -function run_simulation(tspan::Tuple{Float64,Float64}, N::Int) - # Initial conditions - # params, y0 = zebra_conditions(N) # tspan of ~3500 enough - params, y0 = zebra_conditions(N) + # Extract solution and reshape + u_new = reshape(sol[end][1:p.N^2], p.N, p.N) + v_new = reshape(sol[end][p.N^2+1:end], p.N, p.N) - prob = ODEProblem(fhn!, y0, tspan, params) - sol = solve(prob, BS3(), saveat=50.0) # You can try `Rosenbrock23()` too + # Update matrices in-place + U .= u_new + V .= v_new - return sol + 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 diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl index 2c213f0..eeb0ed8 100644 --- a/src/gray_scott_solver.jl +++ b/src/gray_scott_solver.jl @@ -1,6 +1,8 @@ include("utils/constants.jl") include("utils/laplacian.jl") + using Observables + using .Constants using .Laplacian diff --git a/src/visualization.jl b/src/visualization.jl index e50ceac..792a87f 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -1,46 +1,29 @@ module Visualization -include("gray_scott_solver.jl") +#include("gray_scott_solver.jl") +include("fhn_solver.jl") using GLMakie, Observables, Makie -""" - step_through_solution(sol::SolutionType, N::Int) - - Function for visualization for the output of run_simulation - - # Arguments - - `sol`: computed differential equation by run_simulation - - `N`: size of the N×N grid - - # Returns - - ``: Displays created figure -""" -function step_through_solution(sol, N::Int) - fig = Figure(resolution=(600, 600)) - ax = Axis(fig[1, 1]) - slider = Slider(fig[2, 1], range=1:length(sol), startvalue=1) - - # 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=:RdGy) - - # 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 +""" + reset(U, V, heat_obs) + + Resets heatmap to original state by replacing each cell. + Currently only places a square in the center + + # Arguments + - `U`: Matrix filled with ones + - `V`: Empty matrix filled with zeros + - `heat_obs`: Heatmap observable + + # Returns + - ``: resetted heatmap observable +""" function reset!(U, V, heat_obs) U .= 1.0 V .= 0.0 @@ -77,7 +60,7 @@ function build_ui(U, V, param_obs_map::NamedTuple, params_obs, heat_obs) stepsize = Observable(30) spoint = select_point(ax.scene) - # # Controls + # 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) @@ -162,5 +145,5 @@ function build_ui(U, V, param_obs_map::NamedTuple, params_obs, heat_obs) return fig end -export step_through_solution, build_ui +export build_ui end From abc2a4997d6aa6efeb82e854d11da65b3da5da8c Mon Sep 17 00:00:00 2001 From: 2211567 Date: Sat, 14 Jun 2025 14:42:01 +0200 Subject: [PATCH 09/11] adjusted project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index a3758f0..0672cef 100644 --- a/Project.toml +++ b/Project.toml @@ -4,4 +4,3 @@ GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Observables = "510215fc-4207-5dde-b226-833fc4488ee2" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" From a4f09dc79b97ce9b5aa94007abba2cf75e303442 Mon Sep 17 00:00:00 2001 From: 2211567 Date: Sat, 14 Jun 2025 14:48:48 +0200 Subject: [PATCH 10/11] Fixed bug in GS implentation --- scripts/main.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/main.jl b/scripts/main.jl index a95c179..8c35b49 100644 --- a/scripts/main.jl +++ b/scripts/main.jl @@ -37,12 +37,12 @@ param_observables = ( b=b ) -#params_obs = Observable(GSParams(N, dx, Du[], Dv[], F[], k[])) +#params_obs = Observable(Constants.GSParams(N, dx, Du[], Dv[], F[], k[])) params_obs = Observable(Constants.FHNParams(N=N, dx=dx, Du=Du[], Dv=Dv[], ϵ=ϵ[], a=a[], b=b[])) """ lift(Du, Dv, F, k) do u, v, f, ki - params_obs[] = GSParams(N, dx, u, v, f, ki) + params_obs[] = Constants.GSParams(N, dx, u, v, f, ki) end """ From e63fdfde655370d1a70d61ec8c526ceee212f9ab Mon Sep 17 00:00:00 2001 From: 2211567 Date: Sat, 14 Jun 2025 22:09:57 +0200 Subject: [PATCH 11/11] Made universal laplacian function; removed redundant code --- scripts/main.jl | 28 ++++++++++++++-------------- src/fhn_solver.jl | 24 +++++++++++++++++------- src/gray_scott_solver.jl | 28 ++++++++++++++++------------ src/utils/laplacian.jl | 37 ++++++++----------------------------- src/visualization.jl | 4 ++-- 5 files changed, 57 insertions(+), 64 deletions(-) diff --git a/scripts/main.jl b/scripts/main.jl index 8c35b49..9543a01 100644 --- a/scripts/main.jl +++ b/scripts/main.jl @@ -1,6 +1,6 @@ include("../src/utils/constants.jl") -include("../src/fhn_solver.jl") -#include("../src/gray_scott_solver.jl") +#include("../src/fhn_solver.jl") +include("../src/gray_scott_solver.jl") include("../src/visualization.jl") using Observables @@ -9,7 +9,6 @@ using GLMakie using .Constants using .Visualization -""" # GSParams N = 128 dx = 1.0 @@ -22,8 +21,8 @@ param_observables = ( F=F, k=k, ) -""" +""" # FHNParams N = 128 dx = 1.0 @@ -36,23 +35,24 @@ param_observables = ( a=a, b=b ) - -#params_obs = Observable(Constants.GSParams(N, dx, Du[], Dv[], F[], k[])) -params_obs = Observable(Constants.FHNParams(N=N, dx=dx, Du=Du[], Dv=Dv[], ϵ=ϵ[], a=a[], b=b[])) - """ +params_obs = Observable(Constants.GSParams(N, dx, Du[], Dv[], F[], k[])) +#params_obs = Observable(FHNParams(N=N, dx=dx, Du=Du[], Dv=Dv[], ϵ=ϵ[], a=a[], b=b[])) + + lift(Du, Dv, F, k) do u, v, f, ki - params_obs[] = Constants.GSParams(N, dx, u, v, f, ki) + params_obs[] = GSParams(N, dx, u, v, f, ki) +end + +""" +lift(Du, Dv, ϵ, a, b) do d_u, d_v, eps, aa, bb + params_obs[] = FHNParams(N=N, dx=dx, Du=d_u, Dv=d_v, ϵ=eps, a=aa, b=bb) end """ -lift(Du, Dv, ϵ, a, b) do d_u, d_v, eps, aa, bb - params_obs[] = Constants.FHNParams(N=N, dx=dx, Du=d_u, Dv=d_v, ϵ=eps, a=aa, b=bb) -end - U = ones(N, N) V = zeros(N, N) heat_obs = Observable(U) -fig = Visualization.build_ui(U, V, param_observables, params_obs, heat_obs) +fig = build_ui(U, V, param_observables, params_obs, heat_obs) display(fig) diff --git a/src/fhn_solver.jl b/src/fhn_solver.jl index 4118c01..75ca9b3 100644 --- a/src/fhn_solver.jl +++ b/src/fhn_solver.jl @@ -23,16 +23,26 @@ using .Laplacian - `du`: calculated derivatives put back into the du array """ function fhn!(du, u, p::FHNParams, t) - 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 + N = p.N + u_mat = reshape(u[1:N^2], N, N) + v_mat = reshape(u[N^2+1:end], N, N) - Δ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) + Δu = laplacian(u_mat, p.dx) + Δv = laplacian(v_mat, p.dx) - 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) + u_in = u_mat[2:end-1, 2:end-1] + v_in = v_mat[2:end-1, 2:end-1] - du .= vcat(vec(fu), vec(fv)) + fu = p.Du * Δu .+ u_in .- u_in .^ 3 ./ 3 .- v_in + fv = p.Dv * Δv .+ p.ϵ * (u_in .+ p.a .- p.b .* v_in) + + # Construct output with zero boundary (padding) + fu_full = zeros(N, N) + fv_full = zeros(N, N) + fu_full[2:end-1, 2:end-1] .= fu + fv_full[2:end-1, 2:end-1] .= fv + + du .= vcat(vec(fu_full), vec(fv_full)) end function step!(U, V, params_obs::Observable; dx=1) diff --git a/src/gray_scott_solver.jl b/src/gray_scott_solver.jl index eeb0ed8..e9e2c5f 100644 --- a/src/gray_scott_solver.jl +++ b/src/gray_scott_solver.jl @@ -7,24 +7,29 @@ using .Constants using .Laplacian 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 + # Extract parameters + p = params_obs[] + Du, Dv = p.Du, p.Dv + F, k = p.F, p.k + + # Compute Laplacians on the interior + lap_u = laplacian(U, dx) # shape (N-2, N-2) + lap_v = laplacian(V, dx) + + # Extract interior values u = U[2:end-1, 2:end-1] v = V[2:end-1, 2:end-1] + # Gray-Scott update equations (Euler-style) 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) + u_new = u .+ (Du .* lap_u .- uvv .+ F .* (1 .- u)) + v_new = v .+ (Dv .* lap_v .+ uvv .- (F .+ k) .* v) - # Update with new values + # Update interior U[2:end-1, 2:end-1] .= u_new V[2:end-1, 2:end-1] .= v_new - # Periodic boundary conditions + # Apply periodic boundary conditions U[1, :] .= U[end-1, :] U[end, :] .= U[2, :] U[:, 1] .= U[:, end-1] @@ -35,8 +40,7 @@ function step!(U, V, params_obs::Observable; dx=1) V[:, 1] .= V[:, end-1] V[:, end] .= V[:, 2] - # Update heatmap observable - return U + return U # for visualization end function multi_step!(state, n_steps, heat_obs::Observable, params_obs::Observable; dx=1) for _ in 1:n_steps diff --git a/src/utils/laplacian.jl b/src/utils/laplacian.jl index a9cebdb..814573a 100644 --- a/src/utils/laplacian.jl +++ b/src/utils/laplacian.jl @@ -1,35 +1,14 @@ module Laplacian -""" - laplacian(U::Matrix{Float64}, N::Int, h::Float64) - - Computes the discrete 2D Laplacian of a matrix `U` using a 5-point stencil - and circular boundary conditions. - - # Arguments - - `U::Matrix{Float64}`: The input 2D matrix representing the field or image. - - `N::Int`: Integer - - `h::Float64`: The spatial step size or grid spacing between points in the discretization. - - # Returns - - `Vector{Float64}`: A flattened (vectorized) representation of the approximated Laplacian values for each element in `U`. The boundary conditions are handled circularly. -""" -function laplacian(U::Matrix{Float64}, N::Int, h::Float64) - Δ = zeros(N, N) - for i in 2:N-1, j in 2:N-1 - Δ[i, j] = (U[i+1, j] + U[i-1, j] + U[i, j+1] + U[i, j-1] - 4 * U[i, j]) / h^2 - end - return vec(Δ) -end - -function laplacian5(f, dx) +function laplacian(U::AbstractMatrix{<:Real}, dx::Real) 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 + center = U[2:end-1, 2:end-1] + up = U[1:end-2, 2:end-1] + down = U[3:end, 2:end-1] + left = U[2:end-1, 1:end-2] + right = U[2:end-1, 3:end] + + return (up .+ down .+ left .+ right .- 4 .* center) ./ h2 end export laplacian, laplacian5 diff --git a/src/visualization.jl b/src/visualization.jl index 792a87f..c5a6fb0 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -1,6 +1,6 @@ module Visualization -#include("gray_scott_solver.jl") -include("fhn_solver.jl") +include("gray_scott_solver.jl") +#include("fhn_solver.jl") using GLMakie, Observables, Makie