Commit 7568aba3 authored by s.zitz's avatar s.zitz
Browse files

Many changes concerning the computation of forces

parent b28b9f3c
Pipeline #22793 failed with stage
in 2 minutes and 48 seconds
......@@ -135,6 +135,8 @@ notable changes.
- Performance tests for the computation of the capillary pressure, see [pressure.jl](test/pressure.jl).
- New collision kernels `BGKwithBuickGreatedStream!` and `BGKwithGuoStream!`, Buick and Guo use different force corrections.
- Performance tests for the collision kernels, `BGKwithGuoStream!` is the slowest option, further details see [collision.jl](test/collisions.jl).
- The gradient of the pressure is now calculated with a circular padded array and is called `h∇p!`.
### Changed
......@@ -142,6 +144,7 @@ notable changes.
- Collision kernels moved from looping to vectorization with broadcasting, for details see [collision.jl](src/collision.jl).
- The default collision operation was renamed to `BGKandStream!` and has now a flag to use the weigthing function correction.
- Vectorization of equilibrium calculation in function `equilibrium!`, see [equilibrium!.jl](src/equilibria.jl).
- structs are no longer exported, to initialize them use `JuThinFilm.name` with the name of the struct.
### Removed
......@@ -149,6 +152,7 @@ notable changes.
- Kernels have been removed such `kernel_laplacianperiodic!`, `kernel_gradientperiodic!` and `kernel_pressure!` are no longer existing.
- `kernel_equilibrium!` has been removed in favor of a `equilibrium!` GPU version.
- `velocitysquare` for cumoments has been removed.
- `filmpressure!` in favor of `h∇p!`, see [forcing.jl](src/forcing.jl).
## [0.0.1] - 2020-08-18
......
......@@ -12,8 +12,7 @@ include("forcing.jl")
include("runsimulation.jl")
include("measurements.jl")
# Export structs
export params, lattice, culattice, velocity, cuvelocity, force, cuforce, moments, cumoments, pressurestats, collisionparams, distributions, cudistributions
# Export generators for structs
export default_moments, default_cumoments, default_velocity, default_cuvelocity, default_force, default_cuforce
# Defined initial fluid states
export flatsurface, sinewaves, singledroplet
......@@ -24,7 +23,7 @@ export Π, Δ!, capillarypressure!
export equilibrium!
export macroscopicmoments!
export BGKandStream!, BGKwithBuickGreatedStream!, BGKwithGuoStream!
export slippage!, thermalflutuations!, filmpressure!, constantforce!, sumforces!
export slippage!, thermalflutuations!, h∇p!, inclinedplane!, sumforces!
# Put everything together to get a full simulation
export runsimulation, runrealsimulation
......
"""
slippage!(fx, fy, h, vx, vy, col, T)
Computes the force due to substrate friction based on the height `h` and slip length `δ`.
Computes the force due to a substrate friction based on the height `h`, the flow velocity `velocity` and slip length `δ`.
# Theory
Although the model does pretty good job simulating thin film dynamics, i.e. `` \\partial_t h = \\nabla(m(h)\\nabla p)``, it has a few shortcomings.
A rather serious one is the interaction with the substrate.
Our approach is one the hand to use a disjoining potential as used in `capillarypressure` but on the other hand we need to regularize velocity of the three phase contact line.
The later is taken care of with this forcing, which is nothing more than a parabolic interpolation of velocity profiles.
There are two natural boundary conditions, first being the stress free surface, here the velocity of the film is `v`.
However on the substrate depending on the strict choice of the mobility term `m(h)` we would like to observe a vanishing velocity, such a **no-slip** bc.
Now to interpolate between those two we use
`` F_{sl} = a h^2 + b h + c. ``
Which is simply a polynomial of second degree.
The force should be strong if `` h \\rightarrow 0 ``, therefore the full term is given by
`` F_{fric} = \\mathbf{u}\\frac{6h\\mu}{2h^2 + 6h\\delta + 3\\delta^2}. ``
This is calculated every time step and ensures stability of droplet or dewetting like simulations.
# Arguments
- `fx::Array{T,2}`: x-component of the computed force.
......@@ -56,33 +75,36 @@ julia> cuf.x
```
# References
Ideas of how to include slippage and normalize the contact line divergence are not new.
There is plenty of literature dealing with it, Fetzer and Jacobs are definitely explain it nicely.
- [Slippage of Newtonian Liquids: Influence on the Dynamics of Dewetting Thin Films](https://pubs.acs.org/doi/pdf/10.1021/la701746r)
- [Lattice Boltzmann method for thin-liquid-film hydrodynamics](https://journals.aps.org/pre/pdf/10.1103/PhysRevE.100.033313)
The second reference is the paper about the method.
See also: [`collidestreamBGK_periodic!`](@ref)
"""
function slippage!(Fx, Fy, height, velocityx, velocityy, col::collisionparams)
δₛ = 3 * col.δ * col.δ
Fx .= (6 .* col.μ .* height .* velocityx ./
(2 .* height .* height .+
6 .* col.δ .* height .+
δₛ))
Fy .= (6 .* col.μ .* height .* velocityy ./
(2 .* height .* height .+
6 .* col.δ .* height .+
δₛ))
Fx .= (6 .* col.μ .* height .* velocityx) ./ (2 .* height .* height .+
6 .* col.δ .* height .+
δₛ)
Fy .= (6 .* col.μ .* height .* velocityy) ./ (2 .* height .* height .+
6 .* col.δ .* height .+
δₛ)
return nothing
end
# With more structure
function slippage!(F::Force, mom::Moment, col::collisionparams)
δₛ = 3 * col.δ * col.δ
F.x = (6 .* col.μ .* mom.height .* mom.velocity.x ./
(2 .* mom.height .* mom.height .+
6 .* col.δ .* mom.height .+
δₛ))
F.y = (6 .* col.μ .* mom.height .* mom.velocity.y ./
(2 .* mom.height .* mom.height .+
6 .* col.δ .* mom.height .+
δₛ))
F.x .= (6 .* col.μ .* mom.height .* mom.velocity.x) ./ (2 .* mom.height .* mom.height .+
6 .* col.δ .* mom.height .+
δₛ)
F.y .= (6 .* col.μ .* mom.height .* mom.velocity.y) ./ (2 .* mom.height .* mom.height .+
6 .* col.δ .* mom.height .+
δₛ)
return nothing
end
......@@ -91,14 +113,32 @@ end
Computation of the force due to therocapillary waves with thermal energy `kᵦT`.
# Theory
Thermal fluctuations are in principle something that is neglected in lattice Boltzmann.
Actually from a historical point of view noise was a problem as the first lattice Boltzmann algorithms used boolean variables.
Later people spend a lot of effort to include those fluctuations again.
Nevertheless the idea is that due to a finite thermal energy something is wobbling all the time, similar to brownian motion.
In principle lattice Boltzmann can account for thermal fluctuation with the thermilation of non-conserved mode (ghost modes).
This however requires a collision kernel beyond **BGK**, a so called multi relaxation time kernel (MRT) where every mode can be tuned.
Here we choose a different approach, as the model is already rather phenomelogical we include the fluctuations with a fluctuating force.
This force is computed to match the fluctuating term of the fluctuating thin film equation
`` \\partial_t h = \\nabla(m(h)\\nabla p + \\sqrt{k_BT m(h)}\\mathcal{N}).``
Where the fluctuating term `` \\sqrt{k_BT m(h)}\\mathcal{N} `` is the ground assumption for our force matching.
# Arguments
- `F::Force`: force generated due to thermal flucs.
- `mom::Moment`: container for the macroscopic moments, e.g. `height`.
- `rands::Force`: dummy array to store random numbers.
- `Fx::Array`: x-component of the thermally induced fluctuating force.
- `Fy::Array`: y-component of the thermally induced fluctuating force.
- `height::Array`: free surface structure of the fluid.
- `col::collisionparams`: contains run time constants for viscosity `μ`, slippage `δ` and thermal energy `kᵦT`.
- `T::AbstractFloat`: retrun eltype of `F`.
# Examples
```jldoctest
julia> using JuThinFilm, Random, Statistics, Test
......@@ -122,36 +162,76 @@ Test Passed
```
# References
Literature containing the fluctuating thin film equation is growing these days (~2020).
A few good articles to get a general understanding
- [On thermal fluctuations in thin film flow](https://iopscience.iop.org/article/10.1088/0953-8984/17/45/042/pdf)
- [Thin-Film Flow Influenced by Thermal Noise](https://link.springer.com/content/pdf/10.1007/s10955-006-9028-8.pdf)
- [Wetting Phenomena in Nanofluidics](https://www.annualreviews.org/doi/pdf/10.1146/annurev.matsci.38.060407.132451)
See also: [`collidestreamBGK_periodic!`](@ref)
"""
function thermalflutuations!(F::Force, mom::Moment, rands::Force, col::collisionparams, T)
function thermalflutuations!(Fx, Fy, height, col::collisionparams)
# Some random number generation, have to switch between GPU and CPU
if isa(F, cuforce)
rands.x = CUDA.randn(T, size(mom.height))
rands.y = CUDA.randn(T, size(mom.height))
if isa(Fx, CuArray)
Fx .= CUDA.randn(Float64, size(height))
Fy .= CUDA.randn(Float64, size(height))
else
randn!(rands.x)
randn!(rands.y)
randn!(Fx)
randn!(Fy)
end
F.x .= sqrt.(T(2) .* T(col.kᵦT) .* T(col.μ) .* 6 .* mom.height ./
(T(2) .* mom.height .* mom.height .+
T(6) .* mom.height .* T(col.δ) .+
T(3) .* T(col.δ) .* T(col.δ))) .* rands.x
Fx .*= sqrt.(2 .* col.kᵦT .* col.μ .* 6 .* height ./
(2 .* height .* height .+
6 .* height .* col.δ .+
3 .* col.δ .* col.δ))
F.y .= sqrt.(T(2) .* T(col.kᵦT) .* T(col.μ) .* 6 .* mom.height ./
(T(2) .* mom.height .* mom.height .+
T(6) .* mom.height .* T(col.δ) .+
T(3) .* T(col.δ) .* T(col.δ))) .* rands.y
Fy .*= sqrt.(2 .* col.kᵦT .* col.μ .* 6 .* height ./
(2 .* height .* height .+
6 .* height .* col.δ .+
3 .* col.δ .* col.δ))
return nothing
end
function thermalflutuations!(F::Force, mom::Moment, col::collisionparams)
# Some random number generation, have to switch between GPU and CPU
randn!(F.x)
randn!(F.y)
F.x .*= sqrt.(2 .* col.kᵦT .* col.μ .* 6 .* mom.height ./
(2 .* mom.height .* mom.height .+
6 .* mom.height .* col.δ .+
3 .* col.δ .* col.δ))
F.y .*= sqrt.(2 .* col.kᵦT .* col.μ .* 6 .* mom.height ./
(2 .* mom.height .* mom.height .+
6 .* mom.height .* col.δ .+
3 .* col.δ .* col.δ))
return nothing
end
# pressure has to be circular padded for this function!
"""
filmpressure!(f, mom, T)
h∇p!(fx, fy, height, pressure)
Computation of h∇p.
Pressure gradient times the height, the driving force of the *thin film equation*.
Here we compute the gradient of the pressure field, which needs to be circular padded, and multiply the result with the height field.
The gradient is computed over a nine point stencil, that accounts for the *D2Q9* nature of the lattice Boltzmann.
# Theory
Starting with the thin film equation
`` \\partial_t h = \\nabla(m(h)\\nabla p), ``
one immediately sees that `h` can only change due to `` \\nabla p ``, assuming a constant mobility.
The pressure `p` here however is the sum of the laplacian of the height `` \\gamma\\Delta h `` and the wetting potential or disjoining pressure `` \\Pi(h) ``.
For details take a look at the `pressure.jl` file.
# Examples
```jldoctest
julia> using JuThinFilm, CUDA
......@@ -159,7 +239,7 @@ julia> mom = default_moments(Float64,5,5); f = force{Float64}(para = params(lx=5
julia> mom.pressure .= 1.0; mom.height .= 0.5; # Constant pressure yields no gradient!
julia> filmpressure!(f, mom, Float64)
julia> h∇p!(f, mom)
julia> f.x
5×5 Array{Float64,2}:
......@@ -177,7 +257,7 @@ julia> mom.pressure = reshape(collect(1.0:25.0),5,5)
4.0 9.0 14.0 19.0 24.0
5.0 10.0 15.0 20.0 25.0
julia> filmpressure!(f, mom, Float64)
julia> h∇p!(f, mom)
julia> f.x
5×5 Array{Float64,2}:
......@@ -185,7 +265,7 @@ julia> f.x
0.5 0.5 0.5 0.5 0.5
0.5 0.5 0.5 0.5 0.5
0.5 0.5 0.5 0.5 0.5
-0.75 -0.75 -0.75 -0.75 -0.75
-0.75, lx, ly -0.75 -0.75 -0.75 -0.75
julia> f.y
5×5 Array{Float64,2}:
......@@ -195,67 +275,55 @@ julia> f.y
-3.75 2.5 2.5 2.5 -3.75
-3.75 2.5 2.5 2.5 -3.75
```
The GPU case works a little different
```jldoctest
julia> using JuThinFilm, CUDA
julia> mom = default_cumoments(Float32,5,5); f = cuforce{Float32}(para = params{Float32}(lx=5,ly=5));
julia> mom.pressure .= 1.0; mom.height .= 0.5; # Constant pressure yields no gradient!
julia> @cuda blocks=(4,4) threads=(4,4) filmpressure!(f.x, f.y, mom.height, mom.pressure, Float32)
julia> f.x
5×5 CuArray{Float32,2}:
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
julia> mom.pressure = cu(reshape(collect(1.0:25.0),5,5));
julia> @cuda blocks=(4,4) threads=(4,4) filmpressure!(f.x, f.y, mom.height, mom.pressure, Float32)
# References
julia> f.x
5×5 CuArray{Float32,2}:
-0.75 -0.75 -0.75 -0.75 -0.75
0.5 0.5 0.5 0.5 0.5
0.5 0.5 0.5 0.5 0.5
0.5 0.5 0.5 0.5 0.5
-0.75 -0.75 -0.75 -0.75 -0.75
For the gradient of a *D2Q9* lattice Boltzmann field I suggested to use the paper by Junk and Klar
julia> f.y
5×5 CuArray{Float32,2}:
-3.75 2.5 2.5 2.5 -3.75
-3.75 2.5 2.5 2.5 -3.75
-3.75 2.5 2.5 2.5 -3.75
-3.75 2.5 2.5 2.5 -3.75
-3.75 2.5 2.5 2.5 -3.75
```
# References
- [Discretizations for the Incompressible Navier--Stokes Equations Based on the Lattice Boltzmann Method](https://epubs.siam.org/doi/pdf/10.1137/S1064827599357188)
- [Isotropic discrete Laplacian operators from lattice hydrodynamics](https://www.sciencedirect.com/science/article/pii/S0021999112004226?via%3Dihub)
See also: [`collidestreamBGK_periodic!`](@ref)
"""
function filmpressure!(f::Force, mom::Moment)
len, wid = size(mom.pressure)
Threads.@threads for j = 1:wid
@inbounds for i = 1:len
ip = mod1(i+1+len,len)
im = mod1(i-1-len,len)
jp = mod1(j+1+wid,wid)
jm = mod1(j-1-wid,wid)
f.x[i,j] = mom.height[i,j] * (1/3 * (mom.pressure[ip,j] - mom.pressure[im,j]) # up and down
+ 1/12 * (mom.pressure[ip,jp] - mom.pressure[im,jp] - mom.pressure[im,jm] + mom.pressure[ip,jm])) # diagonal
f.y[i,j] = mom.height[i,j] * (1/3 * (mom.pressure[i,jp] - mom.pressure[i,jm] ) # left & right
+ 1/12 * (mom.pressure[ip,jp] + mom.pressure[im,jp] - mom.pressure[im,jm] - mom.pressure[ip,jm])) # diagonal
function h∇p!(fx, fy, height::Array, pressure)
lx, ly = size(height)
diagonalx = 0.0
diagonaly = 0.0
for i = 1:lx
for j = 1:ly
# Compute the diagonal term
diagonalx = pressure[i+1,j+1] - pressure[i-1,j+1] - pressure[i-1,j-1] + pressure[i+1,j-1]
diagonaly = pressure[i+1,j+1] + pressure[i-1,j+1] - pressure[i-1,j-1] - pressure[i+1,j-1]
# Compute the gradient
@inbounds fx[i,j] = height[i,j] * (1/3 * (pressure[i+1,j] - pressure[i-1,j]) + 1/12 * diagonalx)
@inbounds fy[i,j] = height[i,j] * (1/3 * (pressure[i,j+1] - pressure[i,j-1]) + 1/12 * diagonaly)
end
end
return nothing
end
# pressure has to be circular padded for this function!
function h∇p!(fx, fy, height, pressure, lx, ly)
# The index information of padding is lost if cu(h) is used
function h∇p!(fx, fy, height::CuArray, pressure)
lx, ly = size(height)
diagonalx = 0.0
diagonaly = 0.0
for i = 1:lx
for j = 1:ly
ii = 1+1
jj = j+1
# Compute the diagonal term
diagonalx = pressure[ii+1,jj+1] - pressure[ii-1,jj+1] - pressure[ii-1,jj-1] + pressure[ii+1,jj-1]
diagonaly = pressure[ii+1,jj+1] + pressure[ii-1,jj+1] - pressure[ii-1,jj-1] - pressure[ii+1,jj-1]
# Compute the gradient
@inbounds fx[i,j] = height[i,j] * (1/3 * (pressure[ii+1,jj] - pressure[ii-1,jj]) + 1/12 * diagonalx)
@inbounds fy[i,j] = height[i,j] * (1/3 * (pressure[ii,jj+1] - pressure[ii,jj-1]) + 1/12 * diagonaly)
end
end
return nothing
end
function h∇p!(f::Force, mom::Moment, pressure)
lx, ly = size(mom.height)
diagonalx = 0.0
diagonaly = 0.0
for i = 1:lx
......@@ -264,19 +332,27 @@ function h∇p!(fx, fy, height, pressure, lx, ly)
diagonalx = pressure[i+1,j+1] - pressure[i-1,j+1] - pressure[i-1,j-1] + pressure[i+1,j-1]
diagonaly = pressure[i+1,j+1] + pressure[i-1,j+1] - pressure[i-1,j-1] - pressure[i+1,j-1]
# Compute the gradient
@inbounds fx[i,j] = height[i,j] * (1/3 * (pressure[i+1,j] - pressure[i-1,j]) + 1/12 * diagonalx)
@inbounds fy[i,j] = height[i,j] * (1/3 * (pressure[i,j+1] - pressure[i,j-1]) + 1/12 * diagonaly)
@inbounds f.x[i,j] = mom.height[i,j] * (1/3 * (pressure[i+1,j] - pressure[i-1,j]) + 1/12 * diagonalx)
@inbounds f.y[i,j] = mom.height[i,j] * (1/3 * (pressure[i,j+1] - pressure[i,j-1]) + 1/12 * diagonaly)
end
end
return nothing
end
"""
constantforce(F, mom, t, mag, x, y, T)
inclinedplane!(F, mom, t, mag, x, y, tonset=1000, tsmooth=10000)
Forcing as if the film was placed on an inclined plane with inclination direction `x` and `y`.
Applies a constant negative body force to the liquid, smoothed by a `tanh()`.
The inclination angle and strength of the gravitational pull can be tuned with `mag`.
With the boolean variables `x` and `y` one can choose the direction of the force.
# Theory
TODO: Write the theory part
# Arguments
- `F::Force`: output forcing due to the chosen forcing strength `mag` and the height profile `mom.height`.
- `mom::Moment`: container for the height profile, the applied force is roughly h*mag.
- `t::Int`: time step of the simulation, for `t` ≪ 1000 no force is applied.
......@@ -322,7 +398,7 @@ Test Passed
See also: [`collidestreamBGK_periodic!`](@ref), [`sumforce!`](@ref)
"""
function constantforce!(F::Force, mom::Moment, t::Int, mag, x::Bool, y::Bool; tonset=1000, tsmooth=10000)
function inclinedplane!(F::Force, mom::Moment, t::Int, mag, x::Bool, y::Bool; tonset=1000, tsmooth=10000)
smoother = 0.0
smoother = (0.5 .+ 0.5 .* tanh((t - tonset)/tsmooth))
F.x .= x .* (-1) .* mom.height .* mag .* smoother
......@@ -331,6 +407,15 @@ function constantforce!(F::Force, mom::Moment, t::Int, mag, x::Bool, y::Bool; to
return nothing
end
function inclinedplane!(Fx, Fy, height, t::Int, mag, x::Bool, y::Bool; tonset=1000, tsmooth=10000)
smoother = 0.0
smoother = (0.5 .+ 0.5 .* tanh((t - tonset)/tsmooth))
Fx .= x .* (-1) .* height .* mag .* smoother
Fy .= y .* (-1) .* height .* mag .* smoother
return nothing
end
"""
sumforces!(F, pressure, slip, thermal, body)
......
@testset "Different forcings" begin
params64 = JuThinFilm.params(lx=5,ly=5)
params32 = JuThinFilm.params{Float32}(lx=5, ly=5)
@testset "Slippage" begin
f = force{Float64}(para=params(lx=5,ly=5))
col = collisionparams{Float64}()
mom = default_moments(Float64, 5, 5)
lx = 5
ly = 5
f = default_force(Float64, lx, ly)
col = JuThinFilm.collisionparams{Float64}()
mom = default_moments(Float64, lx, ly)
mom.velocity.x .= 0.1
mom.velocity.y .= 0.2
slippage!(f.x, f.y, mom.height, mom.velocity.x, mom.velocity.y, col)
@test all(f.x .== 0.1/11)
@test all(f.y .== 0.2/11)
@test all(f.x . 1/110)
@test all(f.y . 2/110)
mom.height .= 0.1
slippage!(f.x, f.y, mom.height, mom.velocity.x, mom.velocity.y, col)
for i in eachindex(f.x)
@test f.x[i] 1/(2*181) atol=1e-9
@test f.y[i] 1/181 atol=1e-9
end
@test all(f.x . 1/(2*181))
@test all(f.y . 1/181)
f = force{Float64}(para=params(lx=5,ly=5))
f = default_force(Float64, lx, ly)
slippage!(f, mom, col)
for i in eachindex(f.x)
@test f.x[i] 1/(2*181) atol=1e-9
@test f.y[i] 1/181 atol=1e-9
end
@test all(f.x . 1/(2*181))
@test all(f.y . 1/181)
cuf = cuforce{Float64}(para=params{Float64}(lx=5,ly=5))
cuf = default_cuforce(Float64, lx, ly)
cumom = default_cumoments(Float64,5,5)
slippage!(cuf, cumom, col)
@test all(Array(cuf.x) .== 0.0)
......@@ -32,112 +33,96 @@
cumom.velocity.x .= 0.1
cumom.velocity.y .= 0.2
slippage!(cuf, cumom, col)
for i in eachindex(cuf.x)
@test Array(cuf.x)[i] . 0.1/11 atol = 1e-9
@test Array(cuf.y)[i] . 0.2/11 atol = 1e-9
end
slippage!(cuf.x, cuf.y, cumom.height, cumom.velocity.x, cumom.velocity.y, col)
@test all(Array(cuf.x) . 0.1/11)
@test all(Array(cuf.y) . 0.2/11)
end
@testset "Thermal fluctuations no kbt CPU" begin
@testset "Thermal fluctuations" begin
lx = 100
ly = 100
f = force{Float64}(para=params(lx=lx,ly=ly))
r = force{Float64}(para=params(lx=lx,ly=ly))
col = collisionparams{Float64}(kᵦT=0.0)
f = default_force(Float64,lx,ly)
mom = default_moments(Float64, lx, ly)
thermalflutuations!(f, mom, r, col, Float64)
# The mean should be somewhere near to zero
@test sum(r.x)/(lx*ly) 0.0 atol = 1e-1
@test sum(r.y)/(lx*ly) 0.0 atol = 1e-1
# Without kᵦT there are no fluctuations
@test all(f.x .== 0.0)
@test all(f.y .== 0.0)
@testset "CPU no k_BT" begin
col = JuThinFilm.collisionparams{Float64}(kᵦT=0.0)
thermalflutuations!(f, mom, col)
# Without kᵦT there are no fluctuations
@test all(f.x .== 0.0)
@test all(f.y .== 0.0)
end
@testset "CPU with k_BT" begin
col = JuThinFilm.collisionparams{Float64}(kᵦT=11/6, μ=2.0)
thermalflutuations!(f, mom, col)
# With kᵦT
@test sum(f.x)/(lx*ly) 0.0 atol = 1e-1
@test sum(f.y)/(lx*ly) 0.0 atol = 1e-1
@test var(f.x) 4.0 atol = 0.2
@test var(f.y) 4.0 atol = 0.2
end
@testset "GPU with k_BT" begin
col = JuThinFilm.collisionparams{Float32}(kᵦT=11.0f0/6.0f0, μ=2.0f0)
thermalflutuations!(f.x, f.y, mom.height, col)
# With kᵦT
@test sum(f.x)/(lx*ly) 0.0 atol = 1e-1
@test sum(f.x)/(lx*ly) 0.0 atol = 1e-1
@test var(Array(f.x)) 4.0 atol = 2e-1
@test var(Array(f.y)) 4.0 atol = 2e-1
end
end
@testset "Thermal fluctuations with kbt CPU" begin
lx = 100
ly = 100
f = force{Float64}(para=params(lx=lx,ly=ly))
r = force{Float64}(para=params(lx=lx,ly=ly))
col = collisionparams{Float64}(kᵦT=11/6, μ=2.0)
mom = default_moments(Float64, lx, ly)
@testset "pressure gradient h∇p" begin
@testset "gradient of collect CPU" begin
lx = 5
ly = 5
mom = default_moments(Float64, lx, ly)
f = default_force(Float64,lx,ly)
mom.pressure = reshape(collect(1.0:25.0),5,5)
pres = padarray(mom.pressure, Pad(:circular,1,1))
h∇p!(f, mom, pres)
resultx = ones(5,5)
resultx[1,1] = -1.5; resultx[1,2] = -1.5; resultx[1,3] = -1.5; resultx[1,4] = -1.5; resultx[1,5] = -1.5;
resultx[5,1] = -1.5; resultx[5,2] = -1.5; resultx[5,3] = -1.5; resultx[5,4] = -1.5; resultx[5,5] = -1.5;
thermalflutuations!(f, mom, r, col, Float64)
# The mean should be somewhere near to zero
@test sum(r.x)/(lx*ly) 0.0 atol = 1e-1
@test sum(r.y)/(lx*ly) 0.0 atol = 1e-1
# With kᵦT
@test sum(f.x)/(lx*ly) 0.0 atol = 1e-1
@test sum(f.y)/(lx*ly) 0.0 atol = 1e-1
@test var(f.x) 4.0 atol = 0.2
@test var(f.y) 4.0 atol = 0.2
end
@testset "Thermal fluctuations with kbt GPU" begin
lx = 100
ly = 100
f = cuforce{Float32}(para=params{Float32}(lx=lx,ly=ly))
r = cuforce{Float32}(para=params{Float32}(lx=lx,ly=ly))
col = collisionparams{Float32}(kᵦT=11.0f0/6.0f0, μ=2.0f0)
mom = default_cumoments(Float32, lx, ly)
thermalflutuations!(f, mom, r, col, Float64)
# The mean should be somewhere near to zero
@test sum(r.x)/(lx*ly) 0.0 atol = 1e-1
@test sum(r.y)/(lx*ly) 0.0 atol = 1e-1
# With kᵦT
@test sum(f.x)/(lx*ly) 0.0 atol = 1e-1
@test sum(f.x)/(lx*ly) 0.0 atol = 1e-1
@test var(Array(f.x)) 4.0 atol = 2e-1
@test var(Array(f.y)) 4.0 atol = 2e-1
end
@testset "Filmpressure CPU" begin
mom = default_moments(Float64, 5, 5)
f = force{Float64}(para=params(lx=5,ly=