Smoke Simulation
A simulation running on the GPU using SchroedingersSmoke.jl:
using SchroedingersSmoke, CLArrays
using Colors, GPUArrays, GeometryTypes, GLAbstraction
# can be any (GPU) Array type supporting the GPU Array interface and the needed intrinsics
# like sin, etc
ArrayType = CLArray
vol_size = (4,2,2)# box size
dims = (64, 32, 32) .* 2 # volume resolution
hbar = 0.1f0 # Planck constant
dt = 1f0/48f0 # time step
jet_velocity = (1f0, 0f0, 0f0)
nozzle_cen = Float32.((2-1.7, 1-0.034, 1+0.066))
nozzle_len = 0.5f0
nozzle_rad = 0.5f0
n_particles = 500 # number of particles
isf = ISF{ArrayType, UInt32, Float32}(vol_size, dims, hbar, dt);
# initialize psi
psi = ArrayType([(one(Complex64), one(Complex64) * 0.01f0) for i=1:dims[1], j=1:dims[2], k=1:dims[3]]);
psi .= normalize_psi.(psi);
kvec = jet_velocity ./ hbar;
omega = sum(jet_velocity.^2f0) / (2f0*hbar);
isjetarr = isjet.(isf.positions, (nozzle_cen,), (nozzle_len,), (nozzle_rad,))
# constrain velocity
for iter = 1:10
restrict_velocity!(isf, psi, kvec, isjetarr, 0f0)
pressure_project!(isf, psi)
end
particles = ArrayType(map(x-> (0f0, 0f0, 0f0), 1:(10^6)))
add_particles!(particles, 1:n_particles, nozzle_cen, nozzle_rad)
using GLVisualize;
w = glscreen(color = RGBA(0f0, 0f0, 0f0, 0f0), resolution = (1920, 1080));
lincolor = RGBA{Float32}(0.0f0,0.74736935f0,1.0f0,0.1f0)
offset = translationmatrix(Vec3f0(0))
dircolor = map(x-> RGBA{Float32}(x[1], x[2], x[3], 0.1), isf.velocity);
dircolorcpu = Array{RGBA{Float32}}(size(isf.velocity))
copy!(dircolorcpu, dircolor)
dircolorviz = visualize(
dircolorcpu, :absorption, model = offset,
dimensions = Vec3f0(vol_size)
).children[]
# first view sets up camera
_view(dircolorviz, camera = :perspective, position = Vec3f0(4, -7, 2.1), lookat = Vec3f0(4, 0, 2))
_view(visualize(
AABB(Vec3f0(0), Vec3f0(vol_size)), :lines,
color = lincolor,
model = offset,
), camera = :perspective)
_view(visualize(
"Velocity as Colors",
start_position = Point3f0(0), billboard = false,
model = translationmatrix(Vec3f0(0.1, -0.1, -0.2)) * rotationmatrix_x(Float32(0.5pi)),
relative_scale = 0.2f0, color = RGBA(1f0,1f0,1f0,1f0)
), camera = :perspective)
offset = translationmatrix(Vec3f0(0, 0, 2))
particlescpu = Array(particles)
particle_vis = visualize(
(GeometryTypes.Circle(GeometryTypes.Point2f0(0), 0.006f0), reinterpret(Point3f0, particlescpu)),
boundingbox = nothing, # don't waste time on bb computation
model = offset,
color = fill(RGBA{Float32}(0, 0, 0, 0.09), length(particles)),
billboard = true
).children[]
_view(particle_vis, camera = :perspective)
_view(visualize(
AABB(Vec3f0(0), Vec3f0(vol_size)), :lines,
color = lincolor,
model = offset,
), camera = :perspective)
particle_vis[:color][1:n_particles] = map(1:n_particles) do i
xx = (i / n_particles) * 2pi
RGBA{Float32}((sin(xx) + 1) / 2, (cos(xx) + 1.0) / 2.0, 0.0, 0.1)
end
_view(visualize(
"Smoke Particle",
start_position = Point3f0(0), billboard = false,
model = translationmatrix(Vec3f0(0.1, -0.1, 4.1)) * rotationmatrix_x(Float32(0.5pi)),
relative_scale = 0.2f0, color = RGBA(1f0,1f0,1f0,1f0)
), camera = :perspective)
offset = translationmatrix(Vec3f0(4, 0, 0))
magnitude = map(x->sqrt(sum(x .* x)), isf.velocity);
magnitudecpu = Array(magnitude)
magnitudeviz = visualize(
magnitudecpu, :mip,
model = offset,
color_map = [RGBA(0f0, 0f0, 0f0, 0f0), RGBA(0.2f0, 0f0, 0.9f0, 0.9f0), RGBA(0.9f0, 0.2f0, 0f0, 0f0)],
dimensions = Vec3f0(vol_size)
).children[]
_view(magnitudeviz, camera = :perspective)
_view(visualize(
AABB(Vec3f0(0), Vec3f0(vol_size)), :lines,
model = offset,
color = lincolor,
), camera = :perspective)
_view(visualize(
"Velocity Magnitude",
start_position = Point3f0(0), billboard = false,
model = translationmatrix(Vec3f0(4.1, -0.1, -0.2)) * rotationmatrix_x(Float32(0.5pi)),
relative_scale = 0.2f0, color = RGBA(1f0,1f0,1f0,1f0)
), camera = :perspective)
offset = translationmatrix(Vec3f0(4, 0, 2))
directions = copy(isf.velocity)
directionscpu = Array(directions)
directionsviz = visualize(
reinterpret(Vec3f0, directionscpu),
model = offset,
color_map = [RGBA(0f0, 0f0, 0f0, 0f0), RGBA(0.2f0, 0f0, 0.9f0, 0.9f0), RGBA(0.9f0, 0.2f0, 0f0, 0f0)],
ranges = map(x-> 0:x, vol_size)
).children[]
_view(directionsviz, camera = :perspective)
_view(visualize(
"Velocity as Vectorfield",
start_position = Point3f0(0), billboard = false,
model = translationmatrix(Vec3f0(4.1, -0.1, 4.1)) * rotationmatrix_x(Float32(0.5pi)),
relative_scale = 0.2f0, color = RGBA(1f0,1f0,1f0,1f0)
), camera = :perspective)
_view(visualize(
AABB(Vec3f0(0), Vec3f0(vol_size)), :lines,
model = offset,
color = lincolor,
), camera = :perspective)
dt = isf.dt; d = isf.d
iter = 1
io, buffer = GLVisualize.create_video_stream("test.mkv", w)
# main simulation loop
while isopen(w)
t = iter * dt
# incompressible Schroedinger flow
schroedinger_flow!(isf, psi)
psi .= normalize_psi.(psi)
pressure_project!(isf, psi)
start = mod((iter - 1) * n_particles + 1, length(particles))
stop = start + n_particles - 1
add_particles!(particles, start:stop, nozzle_cen, nozzle_rad)
particle_vis[:color][start:stop] = map(1:n_particles) do i
xx = (i / n_particles) * 2pi
RGBA{Float32}((sin(xx) + 1) / 2, (cos(xx) + 1.0) / 2.0, mod(iter, 100) / 100, 0.1)
end
#constrain velocity
restrict_velocity!(isf, psi, kvec, isjetarr, omega*t)
pressure_project!(isf, psi)
velocity_one_form!(isf, psi, isf.hbar)
# inplace StaggeredSharp
dinv = inv.(isf.d)
broadcast!((x, y)-> x .* y, isf.velocity, isf.velocity, (dinv,))
staggered_advect!(particles, isf)
copy!(particlescpu, particles)
GLAbstraction.set_arg!(particle_vis, :position, reinterpret(Point3f0, particlescpu))
dircolor .= (x-> RGBA{Float32}(x[1] * 10f0, x[2] * 10f0, x[3] * 10f0, 1.0)).(isf.velocity)
copy!(dircolorcpu, dircolor)
GLAbstraction.set_arg!(dircolorviz, :volumedata, dircolorcpu)
magnitude .= (x->sqrt(sum(x .* x))).(isf.velocity)
copy!(magnitudecpu, magnitude)
GLAbstraction.set_arg!(magnitudeviz, :volumedata, magnitudecpu)
copy!(directionscpu, isf.velocity)
GLAbstraction.set_arg!(directionsviz, :rotation, vec(reinterpret(Vec3f0, directionscpu)))
GLWindow.poll_glfw()
GLWindow.reactive_run_till_now()
GLWindow.render_frame(w)
GLWindow.swapbuffers(w)
GLVisualize.add_frame!(io, w, buffer)
iter += 1
end
close(io)
GLWindow.destroy!(w)
pwd()