Commit 378125f3 authored by s.zitz's avatar s.zitz
Browse files

More efficient calculation of equilibria, reworked tests

parent da1d3d4a
......@@ -25,8 +25,8 @@ as zeroth population and
`` f^{eq}_i = w_i h \\Bigg(\\frac{3}{2}g h + 3\\mathbf{c}\\cdot\\mathbf{u} + \\frac{9}{2}(\\mathbf{c}\\cdot\\mathbf{u})^2 - \\frac{3}{2}\\mathbf{u}^2\\Bigg), ``
for ``i \\neq 0``.
for ``i \\neq 0``. This also what is implemented below.
# Arguments
- `dists::distributions`: container for the three distribution functions.
- `mom::moments`: struct that contains a macroscopic state, e.g. height and velocity.
......@@ -82,112 +82,32 @@ Paul Dellar and Rick Salmon did a lot of work for the shallow water lattice Bolt
See also: [`moments`](@ref), [`velocity`](@ref), [`params`](@ref)
"""
function equilibrium!(dists::Distribution, mom::Moment, lat::Latticevel, para::params)
len, wid = size(mom.height)
usquare = 0.0
ξ = zeros(len, wid)
function equilibrium!(dists::distributions, mom::moments, lat::lattice, para::params)
# ξ = zeros(len, wid) instead of this I will just use the pressure array, and overwrite it later
dists.feq[:, :, 1] = (mom.height .-
5/6 .* para.gravity .* mom.height.^2 .-
2/3 .* mom.height .* velsquare(mom))
for k in 2:9
ξ = lat.ci[k,1] .* mom.velocity.x .+ lat.ci[k,2] .* mom.velocity.y
mom.pressure = lat.ci[k,1] .* mom.velocity.x .+ lat.ci[k,2] .* mom.velocity.y
@inbounds dists.feq[:, :, k] = lat.wi[k] .* mom.height .* (1.5 .* para.gravity .* mom.height .+
3.0 .* ξ .+
4.5 .* ξ .* ξ .-
3.0 .* mom.pressure .+
4.5 .* mom.pressure .* mom.pressure .-
1.5 .* velsquare(mom))
end
return nothing
end
"""
kernel_equilibrium!(feq, height, velocityx, velocityy, ci, wi, parameters)
Computes the equilibrium distribution `feq` of a given macroscopic state `mom` on the GPU.
The current marcoscopic state `mom` is a struct which contains three elements.
First is being `mom.height` the height field as a scalar function of space.
Second being the velocity vector `mom.velocity` which is two dimensional `mom.velocity.x` & `mom.velocity.y`.
Last and for convenience I save the pressure as well in `mom` which can be accessed with `mom.pressure`.
Two of those are needed for the equilibrium distribution, further the lattice velocities `lattice.ci` and weights `lattice.wi`.
Since the `parameters` struct contains the gravity `parameters.gravity` it is needed as well.
# Arguments
- `feq::CuArray{T, 3}`: the equilibrium distribution.
- `height::CuArray{T,2}`: height of the films surface at every lattice point.
- `velocityx::CuArray{T,2}`: x-component of the velocity vector.
- `velocityx::CuArray{T,2}`: x-component of the velocity vector.
- `ci::CuArray{T,2}`: the two dimensional lattice speeds vector.
- `velocityy::CuArray{T,2}`: the D2Q9 lattice weights.
- `parameters::param`: struct that stores lattice size and gravitational acceleration.
- `T::AbstractFloat`: Floating point accuarcy.
# Examples
```jldoctest
julia> using JuThinFilm, CUDA
julia> feq = CUDA.zeros(5, 5, 9);
julia> mom = JuThinFilm.default_cumoments(Float32, 5, 5);
julia> lat = culattice{Float32}(); parameters = params{Float32}(lx=5, ly=5);
julia> @cuda blocks=(4,4) threads=(4,4) kernel_equilibrium!(feq, mom.height, mom.velocity.x, mom.velocity.y, lat.ci, lat.wi, parameters, Float32)
julia> feq[:, :, 1] # Now it should be filled with mom.height
5×5 CuArray{Float32,2}:
1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0
julia> mom.velocity.x = CUDA.fill(0.1f0, (5,5)); # With some small velocity
julia> @cuda blocks=(4,4) threads=(4,4) kernel_equilibrium!(feq, mom.height, mom.velocity.x, mom.velocity.y, lat.ci, lat.wi, parameters, Float32)
julia> feq[:, :, 1]
5×5 CuArray{Float32,2}:
0.993333 0.993333 0.993333 0.993333 0.993333
0.993333 0.993333 0.993333 0.993333 0.993333
0.993333 0.993333 0.993333 0.993333 0.993333
0.993333 0.993333 0.993333 0.993333 0.993333
0.993333 0.993333 0.993333 0.993333 0.993333
julia> feq[:, :, 2] # Here with (ci . u)
5×5 CuArray{Float32,2}:
0.0366667 0.0366667 0.0366667 0.0366667 0.0366667
0.0366667 0.0366667 0.0366667 0.0366667 0.0366667
0.0366667 0.0366667 0.0366667 0.0366667 0.0366667
0.0366667 0.0366667 0.0366667 0.0366667 0.0366667
0.0366667 0.0366667 0.0366667 0.0366667 0.0366667
```
# References
See also: [`cumoments`](@ref), [`cuvelocity`](@ref), [`culattice`](@ref), [`params`](@ref), [`equilibrium`](@ref)
"""
function kernel_equilibrium!(feq, height, velocityx, velocityy, ci, wi, para::params, T)
len = size(height,1)
wid = size(height,2)
c, ξ, velsquare = T(1.0), T(0.0), T(0.0)
@inbounds for i in 1:len
for j in 1:wid
# u * u
vsquare = velocityx[i,j] * velocityx[i,j] + velocityy[i,j] * velocityy[i,j]
# The zeroth distribution function
feq[i, j, 1] = (height[i,j] -
T(5.0/6.0) * T(para.gravity) * height[i,j] * height[i,j]/ (c * c) -
T(2.0/3.0) * height[i,j] / (c * c) * vsquare)
for k in 2:9
# ci * u
ξ = ci[k,1] * velocityx[i,j] + ci[k,2] * velocityy[i,j]
# All other feq's
feq[i, j, k] = wi[k] * height[i,j] * (T(3.0/2.0) * T(para.gravity) * height[i,j] +
T(3.0) * ξ + T(9.0/2.0) * ξ * ξ -
T(3.0/2.0) * vsquare)
end
end
# This is waaaaaaaaayyyyy faster, really like ten times faster!
function equilibrium!(feq, height, velocityx, velocityy, ξ, ci, wi, para::params)
# ξ = zeros(len, wid) instead of this I will just use the pressure array, and overwrite it later
feq[:, :, 1] .= (height .-
5/6 .* para.gravity .* height.^2 .-
2/3 .* height .* (velocityx .* velocityx .+ velocityy .* velocityy))
for k in 2:9
ξ .= ci[k,1] .* velocityx .+ ci[k,2] .* velocityy
@inbounds feq[:, :, k] .= wi[k] .* height .* (1.5 .* para.gravity .* height .+
3.0 .* ξ .+
4.5 .* ξ .* ξ .-
1.5 .* (velocityx .* velocityx .+ velocityy .* velocityy))
end
return nothing
end
......@@ -237,8 +157,4 @@ See also: [`equilibrium!`](@ref)
"""
function velsquare(mom::moments)
res =@. mom.velocity.x * mom.velocity.x + mom.velocity.y * mom.velocity.y
end
function velsquare(mom::cumoments)
res =@. mom.velocity.x * mom.velocity.x + mom.velocity.y * mom.velocity.y
end
\ No newline at end of file
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment