diff --git a/.gitignore b/.gitignore index c487889..7a4896f 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,4 @@ docs/site/ # environment. Manifest.toml +gif diff --git a/Manifest.toml b/Manifest.toml index 6aaa7a8..df84095 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 = "0465e7f271b455d66e099c4ff05fca7f6e5b2b64" +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"] @@ -1665,29 +1665,7 @@ version = "1.12.0" [deps.Optim.weakdeps] MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -[[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.Opus_jll]]using .Laplacian [[deps.OrdinaryDiffEqBDF]] deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqSDIRK", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "StaticArrays", "TruncatedStacktraces"] git-tree-sha1 = "42755bd13fe56e9d9ce1bc005f8b206a6b56b731" diff --git a/Project.toml b/Project.toml index de40aee..0672cef 100644 --- a/Project.toml +++ b/Project.toml @@ -2,4 +2,5 @@ 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" diff --git a/scripts/main.jl b/scripts/main.jl new file mode 100644 index 0000000..e38e373 --- /dev/null +++ b/scripts/main.jl @@ -0,0 +1,28 @@ +include("../src/gray_scott_solver.jl") +include("../src/visualization.jl") +include("../src/utils/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..3668d26 --- /dev/null +++ b/src/gray_scott_solver.jl @@ -0,0 +1,46 @@ +module GrayScottSolver +using Observables +include("utils/constants.jl") +include("utils/laplacian.jl") +using .Constants +using .Laplacian +export step!, multi_step! + +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/constants.jl b/src/utils/constants.jl similarity index 65% rename from src/constants.jl rename to src/utils/constants.jl index d1734e3..c06f46c 100644 --- a/src/constants.jl +++ b/src/utils/constants.jl @@ -15,6 +15,16 @@ 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/laplacian.jl b/src/utils/laplacian.jl similarity index 70% rename from src/laplacian.jl rename to src/utils/laplacian.jl index 6a04d64..a9cebdb 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,4 +20,18 @@ function laplacian(U::Matrix{Float64}, N::Int, h::Float64) Δ[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 +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 diff --git a/src/visualization.jl b/src/visualization.jl index e73d9f1..51c8ccc 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,146 @@ 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) 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) + 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") + 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 + 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 + + +export step_through_solution, build_ui end