Skip to content
GitLab
Commit 776fd029 authored by s.zitz's avatar s.zitz
Browse files

Seems like my idea of an kernel is rather wrong...

parent 0cbdb650
Pipeline #23087 passed with stage
in 50 minutes and 50 seconds
Hide whitespace changes
Inline Side-by-side
notebooks/droplets_and_patterns.jl
View file @ 776fd029
### A Pluto.jl notebook ###
# v0.11.8
# v0.11.14
using Markdown
using InteractiveUtils
......@@ -7,9 +7,6 @@ using InteractiveUtils
# ╔═╡ fbec2be8-e775-11ea-3b9d-a95eac0bb6b7
using DrWatson
# ╔═╡ da280ffa-e77e-11ea-30d5-65eb4b74d8c7
using JuThinFilm
# ╔═╡ 40854412-e76f-11ea-01aa-379d8e332c43
using DataFrames, CSV, Mmap, CodecZlib, Plots, LazySets # Needed to import the data to dataframes
......@@ -27,11 +24,16 @@ All the data can be found in the `data/sims/relaxation_different_shape_angle/` f
The relexation experiments generated those csv files. `measurements.csv` contains the length of the three phase contact line (**TPC**) and the differential surface area. On the other hand `final_hdvars.csv` contains an snapshot of the simulations state at the last time step."
# ╔═╡ 2aabe662-e771-11ea-37c9-0f1ca548971f
measurements = CSV.File(transcode(GzipDecompressor, Mmap.mmap("../data/sims/relaxation_different_shape_angle/measurements.csv.gz"))) |> DataFrame # three phase contact line lengths and surface areas
begin
measurements = CSV.File(transcode(GzipDecompressor, Mmap.mmap("../data/sims/relaxation_different_shape_angle/measurements.csv.gz"))) |> DataFrame # three phase contact line lengths and surface areas
measurements1 = CSV.File(transcode(GzipDecompressor, Mmap.mmap("../data/sims/relaxation_different_shape_angle/measurements_1.csv.gz"))) |> DataFrame
end
# ╔═╡ 377e4c0e-e771-11ea-31d1-8730208b9cad
HD_snaps = CSV.File(transcode(GzipDecompressor, Mmap.mmap("../data/sims/relaxation_different_shape_angle/final_hdvars.csv.gz"))) |> DataFrame # Snapshot of the final state of the simulation
begin
HD_snaps = CSV.File(transcode(GzipDecompressor, Mmap.mmap("../data/sims/relaxation_different_shape_angle/final_hdvars.csv.gz"))) |> DataFrame # Snapshot of the final state of the simulation
HD_snaps1 = CSV.File(transcode(GzipDecompressor, Mmap.mmap("../data/sims/relaxation_different_shape_angle/final_hdvars_1.csv.gz"))) |> DataFrame # Snapshot of the final state of the simulation
end
# ╔═╡ d719478e-e771-11ea-0a0c-95dacc886b4b
md"## Data Analysis
......@@ -81,9 +83,9 @@ end
# ╔═╡ 508d3c22-e772-11ea-3b00-03c64a90ecd2
begin
plot(polygon, fill=:blue)
height = reshape(HD_snaps[Symbol("$(hd_var)_shape_$(shape)_contrast_$(contrast)_vol_$(vol)")],100,100)
height = reshape(HD_snaps1[Symbol("$(hd_var)_shape_$(shape)_contrast_$(contrast)_vol_$(vol)")],256,256)
contourf!(height,aspect_ratio=1)
plot!(xlims = (0,100))
# plot!(xlims = (0,100))
end
# ╔═╡ 3b0a5dca-e7a0-11ea-318a-dd9d83591a23
......@@ -161,12 +163,11 @@ end
# ╔═╡ Cell order:
# ╠═fbec2be8-e775-11ea-3b9d-a95eac0bb6b7
# ╠═da280ffa-e77e-11ea-30d5-65eb4b74d8c7
# ╠═15c6d3ec-e776-11ea-271b-7d32a12e86d5
# ╟─e84a8518-e39e-11ea-1cd9-59bd95e672bf
# ╠═40854412-e76f-11ea-01aa-379d8e332c43
# ╟─2aabe662-e771-11ea-37c9-0f1ca548971f
# ╟─377e4c0e-e771-11ea-31d1-8730208b9cad
# ╠═2aabe662-e771-11ea-37c9-0f1ca548971f
# ╠═377e4c0e-e771-11ea-31d1-8730208b9cad
# ╟─d719478e-e771-11ea-0a0c-95dacc886b4b
# ╟─36aee6d2-e777-11ea-1ffe-4b6172bdffb6
# ╟─1886fee8-e77d-11ea-1e4e-037a987a2cd3
......
scripts/relax_analysis.jl
View file @ 776fd029
......@@ -27,10 +27,10 @@ center = (lx÷2, ly÷2)
measureframe = DataFrame()
dropletframe = DataFrame()
# Geometric parameters for the patterns
boxside = 124
triangleside = 188.5
circlerad = 70
largehalfax, smallhalfax = 50, 98
boxside = 62
triangleside = 94.3
circlerad = 35
largehalfax, smallhalfax = 25, 49
# Some structs that contain mostly default values
parameters = JuThinFilm.params{T}(lx=lx, ly=ly, mt=maxtime)
collisionparameters = JuThinFilm.collisionparams{T}(δ=δ₀)
......@@ -65,17 +65,19 @@ for R₀ in [20 25 30]
Fsum = default_force(T, lx, ly)
Fslip = default_force(T, lx, ly)
Fbody = default_force(T, lx, ly)
# Get the first equilibrium here
equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
copy!(dists.ftemp,dists.feq)
# Start of the main time loop
for t in 1:parameters.mt
# Measurements of surface area and three phase contact line lenght
JuThinFilm.measuresurfaceandenergy!(diffsurf, energy, mom, θₙ)
JuThinFilm.measurecontactlinelength!(contactline, mom)
# IO for progression
if t%500 == 0
if t%10 == 0
# Measurements of surface area and three phase contact line lenght
JuThinFilm.measuresurfaceandenergy!(diffsurf, energy, mom, θₙ)
JuThinFilm.measurecontactlinelength!(contactline, mom)
end
if t%1000 == 0
# IO for progression
mass = round(sum(mom.height),digits=3)
println("time step $t and mass $mass")
end
......@@ -94,8 +96,9 @@ for R₀ in [20 25 30]
end
# Save measurements and the last snapshot of the height field
angle_diff = rad2deg(δᵪ*π)
measureframe["Line_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = contactline ./ contactline[1]
measureframe["Surf_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = diffsurf ./ diffsurf[1]
measureframe["L_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = contactline # ./ contactline[1] but use it afterwards not already!
measureframe["S_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = diffsurf
measureframe["E_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = energy
dropletframe["height_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.height)
dropletframe["velx_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.velocity.x)
dropletframe["vely_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.velocity.y)
......@@ -104,12 +107,12 @@ for R₀ in [20 25 30]
end
end
savepath = "data/sims/relaxation_different_shape_angle/"
results_measures = "measurements_1.csv"
results_drop = "final_hdvars_1.csv"
results_measures = "measurements_2.csv"
results_drop = "final_hdvars_2.csv"
CSV.write(savepath * results_measures, measureframe)
CSV.write(savepath * results_drop, dropletframe)
zip1 = `gzip data/sims/relaxation_different_shape_angle/measurements_1.csv`
zip2 = `gzip data/sims/relaxation_different_shape_angle/final_hdvars_1.csv`
zip1 = `gzip data/sims/relaxation_different_shape_angle/measurements_2.csv`
zip2 = `gzip data/sims/relaxation_different_shape_angle/final_hdvars_2.csv`
run(zip1)
run(zip2)
println("Simulations done results are saved in $savepath")
\ No newline at end of file
src/calcmoments.jl
View file @ 776fd029
......@@ -90,8 +90,8 @@ end
function macroscopicmoments!(mom::cumoments, dists::Distribution)
sum!(mom.height, dists.fout)
mom.velocity.x .= dists.fout[:,:,2] .- dists.fout[:,:,4] .+ dists.fout[:,:,6] .- dists.fout[:,:,7] .- dists.fout[:,:,8] .+ dists.fout[:,:,9] # This is not as easy as thought
mom.velocity.y .= dists.fout[:,:,3] .- dists.fout[:,:,5] .+ dists.fout[:,:,6] .+ dists.fout[:,:,7] .- dists.fout[:,:,8] .- dists.fout[:,:,9]
mom.velocity.x .= view(dists.fout,:,:,2) .- view(dists.fout,:,:,4) .+ view(dists.fout,:,:,6) .- view(dists.fout,:,:,7) .- view(dists.fout,:,:,8) .+ view(dists.fout,:,:,9) # This is not as easy as thought
mom.velocity.y .= view(dists.fout,:,:,3) .- view(dists.fout,:,:,5) .+ view(dists.fout,:,:,6) .+ view(dists.fout,:,:,7) .- view(dists.fout,:,:,8) .- view(dists.fout,:,:,9)
# Normalize with height
mom.velocity.x ./= mom.height
......
src/collisions.jl
View file @ 776fd029
......@@ -424,3 +424,31 @@ function BGKwithGuoStream!(dists::Distribution, mom::Moment, F::Force, col::coll
return nothing
end
# ======================== Untested GPU code ============================== #
function kernel_BGKandStream!(fout, feq, ftemp, Fx, Fy, col::JuThinFilm.collisionparams)
# The lattice speeds
len, wid = size(Fx)
omeg = 1-col.ω
# Loop over all nine populations
@inbounds for i in 1:len
for j in 1:wid
ip = mod1(i+1+len,len)
im = mod1(i-1-len,len)
jp = mod1(j+1+wid,wid)
jm = mod1(j-1-wid,wid)
# This is three times faster than the above!
fout[i,j,1] = omeg * ftemp[i,j,1] + col.ω * feq[i,j,1]
fout[ip,j,2] = omeg * ftemp[i,j,2] + col.ω * feq[i,j,2] + 1/3 * Fx[i,j]
fout[i,jp,3] = omeg * ftemp[i,j,3] + col.ω * feq[i,j,3] + 1/3 * Fy[i,j]
fout[im,j,4] = omeg * ftemp[i,j,4] + col.ω * feq[i,j,4] - 1/3 * Fx[i,j]
fout[i,jm,5] = omeg * ftemp[i,j,5] + col.ω * feq[i,j,5] - 1/3 * Fy[i,j]
fout[ip,jp,6] = omeg * ftemp[i,j,6] + col.ω * feq[i,j,6] + 1/24 * (Fx[i,j] + Fy[i,j])
fout[im,jp,7] = omeg * ftemp[i,j,7] + col.ω * feq[i,j,7] + 1/24 * (Fy[i,j] - Fx[i,j])
fout[im,jm,8] = omeg * ftemp[i,j,8] + col.ω * feq[i,j,8] - 1/24 * (Fx[i,j] + Fy[i,j])
fout[ip,jm,9] = omeg * ftemp[i,j,9] + col.ω * feq[i,j,9] + 1/24 * (Fx[i,j] - Fy[i,j])
end
end
return nothing
end
\ No newline at end of file
src/equilibria.jl
View file @ 776fd029
......@@ -159,4 +159,25 @@ See also: [`equilibrium!`](@ref)
"""
function velsquare(mom::moments)
res =@. mom.velocity.x * mom.velocity.x + mom.velocity.y * mom.velocity.y
end
# Untested CUDA stuff here
function kernel_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
len, wid = size(height)
@inbounds for i in 1:len
for j in 1:wid
feq[i, j, 1] = (height[i,j] -
5/6 * para.gravity * height[i,j]^2 -
2/3 * height[i,j] * (velocityx[i,j] * velocityx[i,j] + velocityy[i,j] * velocityy[i,j]))
for k in 2:9
ξ = ci[k,1] * velocityx[i,j] + ci[k,2] * velocityy[i,j]
feq[i, j, k] = wi[k] * height[i,j] * (1.5 * para.gravity * height[i,j] +
3.0 * ξ +
4.5 * ξ * ξ -
1.5 * (velocityx[i,j] * velocityx[i,j] + velocityy[i,j] * velocityy[i,j]))
end
end
end
return nothing
end
\ No newline at end of file
src/pressure.jl
View file @ 776fd029
......@@ -248,8 +248,8 @@ function capillarypressure!(pressure::Array{Float64,2}, height, θij::Array{Floa
pressure .-= (1 .- map.(cospi,θij)) .* spreadingfactor .* (power_broad.((stats.hstar./(hdisj .+ hᵪ)),stats.n) .- power_broad.((stats.hstar./(hdisj .+ hᵪ)), stats.m))
return nothing
end
#= TODO: Test the GPU implementation
function capillarypressure!(pressure::CuArray{Float64,2}, height, stats::pressurestats)
#= TODO: Test the GPU implementation
function capillarypressure!(pressure::CuArray, height, stats::pressurestats)
# The default compute size
wid, len = size(pressure)
linear = 0.0
......@@ -258,16 +258,16 @@ function capillarypressure!(pressure::CuArray{Float64,2}, height, stats::pressur
for i = 1:wid
for j = 1:len
# Compute the laplacian with periodic boundary conditions, height is padded
ii = i+1
jj = j+1
linear = 4 * (height[ii+1, jj] + height[ii-1, jj] + height[ii, jj+1] + height[ii, jj-1])
diagonal = height[ii+1, jj+1] + height[ii-1, jj+1] + height[ii-1, jj-1] + height[ii+1, jj-1]
@inbounds pressure[i,j] = -stats.γ*(1/6 * (linear + diagonal - 20 * height[ii, jj]))
@inbounds begin
linear = 4 * (height[i+1, j] + height[i-1, j] + height[i, j+1] + height[i, j-1])
diagonal = height[i+1, j+1] + height[i-1, j+1] + height[i-1, j-1] + height[i+1, j-1]
pressure[i,j] = -stats.γ*(1/6 * (linear + diagonal - 20 * height[i, j]))
end
end
end
# Next compute the disjoining term
spreadingfactor = stats.γ * (1 - CUDA.cospi(stats.θ)) * (stats.n-1)*(stats.m-1) / ((stats.n - stats.m)*stats.hstar)
hdisj = @view height[2:wid+1,2:len+1]
spreadingfactor = stats.γ * (1 - cospi(stats.θ)) * (stats.n-1)*(stats.m-1) / ((stats.n - stats.m)*stats.hstar)
hdisj = view(height, 2:wid+1,2:len+1)
pressure .-= spreadingfactor .* (power_broad.((stats.hstar./hdisj),stats.n) .- power_broad.((stats.hstar./hdisj), stats.m))
return nothing
end
......@@ -291,7 +291,55 @@ function capillarypressure!(pressure::CuArray{Float64,2}, θij::CuArray{Float64,
# Next compute the disjoining term
spreadingfactor = stats.γ * (stats.n-1)*(stats.m-1) / ((stats.n - stats.m)*stats.hstar)
hdisj = @view height[2:wid+1,2:len+1]
pressure .-= (1 .- CUDA.cospi.(θij)) .* spreadingfactor .* (power_broad.((stats.hstar./hdisj),stats.n) .- power_broad.((stats.hstar./hdisj), stats.m))
pressure .-= (1 .- cospi.(θij)) .* spreadingfactor .* (power_broad.((stats.hstar./hdisj),stats.n) .- power_broad.((stats.hstar./hdisj), stats.m))
return nothing
end
=#
function kernel_pressure!(p, height, stats::pressurestats)
# working good! :)
len, wid = size(height)
hₜ = 0.05
# Compute a constant disjoininh pressure based on the pressure stats struct
κ = stats.γ * (1 - map(CUDA.cospi, stats.θ)) * (stats.n-1) * (stats.m-1) / ((stats.n - stats.m)*stats.hstar)
@inbounds for j = 1:wid
for i = 1:len
# Find the nearest neighbors
ip = mod1(i+1+len,len)
im = mod1(i-1-len,len)
jp = mod1(j+1+wid,wid)
jm = mod1(j-1-wid,wid)
hstarbyh = stats.hstar / (height[i,j] + hₜ)
# Compute the laplacian with periodic boundary conditions
p[i,j] = -stats.γ*((2/3 * (height[ip, j] + height[im, j] + height[i, jp] + height[i, jm]) +
1/6 * (height[ip, jp] + height[im, jp] + height[im, jm] + height[ip, jm]) -
10/3 * height[i, j])) - κ * (power_broad(hstarbyh,stats.n) - power_broad(hstarbyh,stats.m))
end
end
return nothing
end
function kernel_h∇p!(fx, fy, pressure, height)
# working good! :)
len, wid = size(height)
@inbounds for j = 1:wid
for i = 1:len
# Find the nearest neighbors
ip = mod1(i+1+len,len)
im = mod1(i-1-len,len)
jp = mod1(j+1+wid,wid)
jm = mod1(j-1-wid,wid)
# Compute the diagonal term
diagonalx = pressure[ip,jp] - pressure[im,jp] - pressure[im,jm] + pressure[ip,jm]
diagonaly = pressure[ip,jp] + pressure[im,jp] - pressure[im,jm] - pressure[ip,jm]
# Compute the gradient
fx[i,j] = height[i,j] * (1/3 * (pressure[ip,j] - pressure[im,j]) + 1/12 * diagonalx)
fy[i,j] = height[i,j] * (1/3 * (pressure[i,jp] - pressure[i,jm]) + 1/12 * diagonaly)
end
end
return nothing
end
src/runsimulation.jl
View file @ 776fd029
......@@ -49,49 +49,49 @@ function runsimulation(parameters::JuThinFilm.params, colparams::JuThinFilm.coll
return mom
end
# function runsimulationGPU(parameters::JuThinFilm.params, colparams::JuThinFilm.collisionparams, pressureparams::JuThinFilm.pressurestats; kind="Droplet", h₀=30, θ₀=1/3, nx=1, ny=1, ϵ=0.1)
# # Some unuseful IO
# println("You choose to run on a simulation of a $kind")
# # Some parameter defining at the start
# T = eltype(parameters.gravity)
# JuThinFilm.culattice{T}()
function runsimulationGPU(parameters::JuThinFilm.params, colparams::JuThinFilm.collisionparams, pressureparams::JuThinFilm.pressurestats; kind="Droplet", h₀=30, θ₀=1/3, nx=1, ny=1, ϵ=0.1)
# Some unuseful IO
println("You choose to run on a simulation of a $kind")
# Some parameter defining at the start
T = eltype(parameters.gravity)
lat = JuThinFilm.culattice{T}()
# center = (parameters.lx÷2, parameters.ly÷2)
# if kind == "Droplet"
# mom = singledroplet(T, parameters.lx, parameters.ly, h₀, θ₀, center, device=true)
# elseif kind == "Sinewave"
# mom = sinewaves(T, parameters.lx, parameters.ly, nx=nx, ny=ny, ϵ=ϵ)
# end
center = (parameters.lx÷2, parameters.ly÷2)
if kind == "Droplet"
mom = singledroplet(T, parameters.lx, parameters.ly, h₀, θ₀, center, device=true)
elseif kind == "Sinewave"
mom = sinewaves(T, parameters.lx, parameters.ly, nx=nx, ny=ny, ϵ=ϵ)
end
# dists = default_cudistribution(T, parameters.lx, parameters.ly)
# Fpressure = default_cuforce(T, parameters.lx, parameters.ly)
# Fsum = default_cuforce(T, parameters.lx, parameters.ly)
# Fslip = default_cuforce(T, parameters.lx, parameters.ly)
# Fbody = default_cuforce(T, parameters.lx, parameters.ly)
dists = default_cudistribution(T, parameters.lx, parameters.ly)
Fpressure = default_cuforce(T, parameters.lx, parameters.ly)
Fsum = default_cuforce(T, parameters.lx, parameters.ly)
Fslip = default_cuforce(T, parameters.lx, parameters.ly)
Fbody = default_cuforce(T, parameters.lx, parameters.ly)
# # Get the first equilibrium here
# equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
# copy!(dists.ftemp,dists.feq)
# # Now move to the time loop
# for t in 1:parameters.mt
# if t % 1000 == 0
# mass = round(sum(mom.height),digits=3)
# maxh = maximum(mom.height)
# println("time step $t with mass $mass and maximum height $maxh")
# end
# Get the first equilibrium here
equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
copy!(dists.ftemp,dists.feq)
# Now move to the time loop
for t in 1:parameters.mt
if t % 10 == 0
mass = round(sum(mom.height),digits=3)
maxh = maximum(mom.height)
println("time step $t with mass $mass and maximum height $maxh")
end
# # Compute forces
# capillarypressure!(mom.pressure, padarray(mom.height, Pad(:circular, 1, 1)), pressureparams)
# h∇p!(Fpressure.x, Fpressure.y, mom.height, padarray(mom.pressure, Pad(:circular, 1, 1)))
# slippage!(Fslip.x, Fslip.y, mom.height, mom.velocity.x, mom.velocity.y, colparams)
# inclinedplane!(Fbody.x, Fbody.y, mom.height, t, 1e-4, true, false)
# sumforces!(Fsum, Fpressure, Fslip, Fbody)
# Compute forces
@cuda blocks=(8,8) threads=(8,8) JuThinFilm.kernel_pressure!(mom.pressure, mom.height, pressureparams)
# @cuda blocks=(4,4) threads=(4,4) JuThinFilm.kernel_h∇p!(Fpressure.x, Fpressure.y, mom.pressure, mom.height)
slippage!(Fslip.x, Fslip.y, mom.height, mom.velocity.x, mom.velocity.y, colparams)
inclinedplane!(Fbody.x, Fbody.y, mom.height, t, 1e-4, true, false)
sumforces!(Fsum, Fpressure, Fslip, Fbody)
# # Get the equilibrium and then performe a collision.
# equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
# BGKandStream!(dists, Fsum, colparams, kind="WFM")
# macroscopicmoments!(mom, dists)
# end
# return mom
# end
# Get the equilibrium and then performe a collision.
# @cuda blocks=(4,4) threads=(4,4) JuThinFilm.kernel_equilibrium!(dists.feq, mom.height, mom.velocity.x, mom.velocity.y, mom.pressure, lat.ci, lat.wi, parameters)
# @cuda blocks=(4,4) threads=(4,4) JuThinFilm.kernel_BGKandStream!(dists.fout, dists.feq, dists.ftemp, Fsum.x, Fsum.y, colparams)
# macroscopicmoments!(mom, dists)
end
return mom
end
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