created utils file; rm two files
2211567
parent
f39c85a0ad
commit
c98b5f0fca
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
154
src/solver.jl
154
src/solver.jl
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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[])")
|
||||
|
|
|
|||
Loading…
Reference in New Issue